© B. Braconnier et al., published by IFP Energies nouvelles, 2019

List of symbols

Symbols concerning the rock

Φ: Rock porosity

K : Rock permeability

ρ _R : Rock density

τ : Rock tortuosity

Symbols concerning the phase ψ

ρ_ψ: Mass density

S _ψ : Saturation

u _ψ : Darcy velocity

P _ψ : Pressure

μ _ψ : Viscosity

kr _ψ : Relative permeability

$k{r}_{\psi }^0$ : Maximum relative permeability

n _ψ : Corey exponent

λ _ψ : Mobility

v _ψ : Hydrodynamic gradient

q _ψ : Well source terms per volume

$D_{\psi }^{\mathrm{ds}}$ : Dispersion coefficient

Symbols concerning the chemical tracer α

ψ _α : Reference phase

M _α : Molar mass

n _α : Mole number

C _α,ψ : Mole concentration in the phase ψ

C _α,R : Adsorbed mole concentration

$C_{\alpha ,R}^{\max }$ : Maximum adsorbed mole concentration

b _α : Langmuir coefficient

$D_{\alpha ,\psi }^{\mathrm{df}}$ : Molecular diffusion in phase ψ

F _α : Molar advection flux

J _α : Molar diffusion/dispersion fluxes

$S_{\alpha }^I$ : Injection source term

$C_{\alpha ,{\psi }_{\alpha }}^{\mathrm{inj}}$ : Injected mole concentration

$S_{\alpha }^{\mathrm{P}}$ : Production source term

$S_{\alpha }^{\mathrm{Deg}}$ : Degradation source term

${\lambda }_{\alpha }^{\mathrm{Deg}}$ : Degradation kinetic constant

$t_{\alpha ,1/2}^{\mathrm{Deg}}$ : Degradation half life time

$S_{\alpha }^{\mathrm{Rea}}$ : Reaction source term

K _α,ψ  − ψ _α : Partition coefficient for phase ψ

Symbols concerning the reaction $\mathrm{Rea} (\mathcal{R}+\alpha \to \mathcal{P}+{\gamma }_{\alpha ,\beta }\beta )\mathrm{ :}$

$\mathcal{R}$ : Chemical reaction reactant

$\mathcal{P}$ : Chemical reaction product

γ _α,β : Stoechiometric coefficient

$t_{\alpha \beta ,1/2}^{\mathrm{Rea}}$ : Chemical reaction half life time

${\lambda }_{\alpha \beta }^{\mathrm{Rea}}$ : Chemical reaction kinetic constant

Symbols concerning the mesh ${\mathcal{M}}_h$

Π_prod : Set of perforations of production wells

Π_inj : Set of perforations of injection wells

Symbols concerning the mesh cell K

M _K : Cell center

|K|: Cell measure

$\mathcal{N}\left(K\right)$ : Set of neighbors cells

KL : Face between cell $K$ and cell $L$

n _KL : Unit normal vector to face KL

d _K,KL : Distance between M _K and face KL

Π _p (K): Set of production perforations for cell K

Π _i (K): Set of injection perforations for cell K

Symbols concerning the numerical scheme

t ⁿ : n ^th time step

Δt : Time step increment

${\mathbf{v}}_{\psi ,\mathrm{KL}}^n$ : Finite volume flux of hydrodynamic gradient at face KL

$\mathcal{A}$ : Discrete adsorption operator

$\mathcal{F}$ : Numerical scheme operator in concise form

${\nabla }_{\mathrm{KL}}{C}_{\alpha ,\psi }^n$ : Gradient of $(C_{\alpha ,\psi }{)}^n$ at face KL

${\tilde{\nabla }}_{\mathrm{KL}}{C}_{\alpha ,\psi }^n$ : Limited gradient of $(C_{\alpha ,\psi }{)}^n$ at face KL

$({\tilde{C}}_{\alpha ,\psi }^n{)}_{\mathrm{KL}}$ : Second-order approximation of $(C_{\alpha ,\psi }{)}_K^n$ at face KL

$\mathcal{L}$ : Slope limiter

O(1): First order scheme for chemical tracer

O(2^m): Second order scheme with Minmod limiter for chemical tracer

O(2^s): Second order scheme with Superbee limiter for chemical tracer

Misc symbols

g : Gravity acceleration

C _h : Heavy components mass fraction in oil phase

C _v : Volatile components mass fraction in oil phase

S : Normalized (mobile) water saturation

S _wi : (Immobile) irreducible water saturation

S _orw : Residual oil saturation to waterflooding

t _a : Alcohol arrival time

t _e : Ester arrival time

1 Introduction

The Single Well Chemical Tracer Test (SWCTT) method has been the subject of longstanding patent of Harry Deans [1] on behalf of Esso and is detailed in [2]. Main application of the SWCTT method is the determination of the residual oil saturation of a reservoir formation after a secondary or tertiary waterflooding process. This measurement is essential for assessing the process effectiveness under conditions as representative as possible, i.e. on a large reservoir rock volume of 10–20 feet around the test wells. The SWCTT method thus consists of injecting an aqueous slug of a reactive tracer and then pushing it by water at a certain distance, currently a few meters, from the well. This tracer is called the primary tracer. By reacting slowly with water, the primary tracer generates in situ another tracer called the secondary tracer. These primary and secondary tracers do not have the same behavior in the presence of oil. The primary tracer is called miscible, i.e. miscible in oil, because it undergoes partitioning between the aqueous phase and the residual oil. The secondary tracer is called immiscible because it is generally assumed to remain entirely in the aqueous phase. However, that assumption is not valid for all tracers. Production is resumed at the end of an optimal duration in order to obtain a sufficient amount of both tracers in place. The analysis of the produced water reveals a time shift or chromatographic delay of the miscible tracer with respect to the immiscible tracer because, unlike the latter, the miscible tracer is gradually restored by the immobile residual oil by succession of equilibria with the water devoid of this tracer and coming from the distant parts of the well. In recent years, the SWCTT method has gained a renewed interest along with EOR methods because it gives the possibility to measure in situ the efficiency of EOR process through the comparison of SWCT tests performed before and after EOR implementation. One may refer to the comprehensive framework proposed by Khaledialidusti and Kleppe [3] to get an insight of SWCTT stages and design parameters.

This paper deals with the design of a numercial scheme devoted to the SWCTT method simulation to be integrated in industrial simulators. This scheme must rely on finite volume method because widely used by industrial simulators [4, 5], must be almost easy to implement and must improve significantly the results accuracy in comparison with a first-order upwind scheme [5]. Regarding the flow model, we adopt a classical black-oil system [6] upgraded with the chemical tracers mole conservation, phase partitioning, adsorption, degradation and chemical reaction. The chemical tracer has no signicant effect on the fluid phases present in the reservoir so that the modeling of these phases does not need to be modified. The resulting system is quite simple and has already been addressed in more complex context, in particular in the works of Delshad et al. [7], Pope and Nelson [8], Datta-Gupta et al. [9] where other chemicals such as polymer, surfactant or foam were involved. Also, many sophisticated high-order numerical schemes have been proposed in the literature, such as [10, 11], but cannot be implemented quickly in our numerical simulator framework. The goal of this work is to propose an easy and fast to implement numerical scheme dedicated to the SWCTT method simulation only, with an improved accuracy compared to classical first-order upwind finite volume methods [12]. For this, we select among well-known and widely used numerical methods the most suitable one on the basis of studying carefully the SWCTT method and its underlying physics.

The chemical tracers have no signicant effect on the fluids present in the reservoir such that the classical black-oil system [6] is unchanged. Thus, the resolution of chemical tracers concentrations is decoupled from the resolution of pressure and saturations. In the context of a SWCTT process, the oil and gas quantities remaining in the reservoir are supposed to be low when the tracers are injected: the oil saturation is closed to the residual oil saturation and the gas has already been produced. Thus, during the SWCTT method, the saturations changes are negligible and the pressure profiles are smooth. In this context, there is no need to consider a high-order scheme for the resolution of these variables and we consider a first-order implicit in time and upwind scheme [12] because robust and unconditionally stable.

An important physical feature of the SWCTT method is that physical phenomena occur sequentially with two different characteristic specific times. During the injection and production of the chemical tracer slugs, the main physical process is advection. The duration of these advection steps is short so that the chemical reaction may be ignored in first approximation during those steps. This is indeed an approximation because the chemical reaction is initiated at the beginning of tracer injection and ongoing all during the test. In addition, peculiar reservoir conditions such as a high temperature and a low buffer capacity may emphasize the role played by reaction during the advection steps of SWCTT, as shown by Khaledialidusti and Kleppe in [13] for the hydrolysis of an ester. During the rest time, the chemical tracers are immobile and the primary tracer generates a secondary tracer according to certain kinetics. In such a context, it is not necessary to consider a complex scheme that couples these two phenomena, namely the advection and the chemical reactions.

As using the SWCTT method requires to determine the delay between several chemical tracer breakthrough, it is valuable to use a precise numerical scheme. As in practice, the SWCTT method takes place on a localized part of a reservoir, namely 10–20 feet around a testing well, we assume here that the simulation model would entail a reasonable number of mesh cells such that an explicit scheme for the chemical tracers is appropriate. In fact, an explicit scheme has two main advantages. It induces less numerical diffusion and it is easier to implement than an implicit scheme. However a CFL stability criterion has to be satisfied. In order to improve the accuracy of the SWCTT method simulation, in particular the calculation of sharp concentration profiles, we propose the use of a second-order scheme in time and space. For the second-order in time, we propose to use a predictor-corrector scheme [14, 15] (or also referred as the Heun’s method) which is quite simple to implement. For the second-order in space, we suggest to use MUSCL technique [16– 20] because it is simple to implement and robust. This well known technique has been recently applied in the field of chemical enhanced oil recovery in the work [21– 23] by considering fluxes reconstruction and limitation. In our context, the second-order is obtained by computing slopes to reconstruct chemical tracer mole concentration at each cell interface to compute the numerical tracer fluxes more accurately. In the framework of the MUSCL scheme, the slopes are limited with slope limiter functions to avoid spurious oscillations of the resulting solutions and to preserve the numerical scheme stability [19, 24, 25].

The paper is organized as follows. First, we present the physical model governing the flow. It is composed of a black oil subsystem upgraded with the chemical tracer mole conservation equation. Then, we describe our numerical scheme that is based on a fractional time-stepping approach with an implicit resolution of the black oil subsystem and an explicit resolution of the chemical tracers mole conservation equation. Afterwards, this basic explicit numerical scheme for the chemical tracer mole conservation equation is extended to second-order in time and space by using the MUSCL technique. Finally, we propose a realistic example of SWCTT simulation model to illustrate the benefit and robustness of our easy and fast to implement numerical scheme. Finally, on the basis on a radial 1-D mesh convergence study, the benefit of a second-order scheme is discussed regarding the assessment of residual oil saturation.

2 Physical model

We consider a model for a three-phase flow in a porous medium that involves several chemical tracers. We distinguish three phases: a liquid phase W consisting of water; a second liquid phase O consisting of oil and a gas phase G. This flow involves chemical tracers in view of performing a SWCTT. The chemical tracers are either mobile or fixed and adsorbed on the rock. When they are mobile, they can be transported, depending on their physical properties, in the three phases, as a result of advection, diffusion and dispersion effects. These chemical tracers can degrade and react with each other. In addition, as required in the SWCTT process, the chemical tracers have no significant effect on any of the three phases.

In this section we recall briefly the physical model for the three phases and then we describe the chemical physical model in detail.

2.1 Black-oil model

To model the water and hydrocarbon components flow, we consider a standard black-oil model [6, 12] which is defined in mass per unit volume of saturated porous medium as $$\{\begin{array}{l}{\mathrm{\partial }}_t\left(\mathrm{\Phi }{C}_h{\rho }_O{S}_O\right)\mathrm{ }+\nabla \cdot \left({\rho }_O{C}_h{\mathbf{u}}_O\right)+C_h{\rho }_O{q}_O=0,\\ \\ {\mathrm{\partial }}_t(\mathrm{\Phi }{\rho }_G{S}_G+\mathrm{\Phi }{C}_v{\rho }_O{S}_O)+\nabla \cdot \left({\rho }_O{C}_v{\mathbf{u}}_O+{\rho }_G{\mathbf{u}}_G\right)\\ +{\rho }_G{q}_G+C_v{\rho }_O{q}_O=0,\\ \\ {\mathrm{\partial }}_t\left(\mathrm{\Phi }{\rho }_W{S}_W\right)+\nabla \cdot \left({\rho }_W{\mathbf{u}}_W\right)+{\rho }_W{q}_W=0,\end{array}$$(1) where Φ is the rock porosity. For each phase denoted ψ = W, O, G, S _ψ is the saturation, ρ _ψ is the mass density and q _ψ denotes the source term per volume. It will be useful to distinguish the outflow source terms $q_{\psi }^{+}\ge 0$ modeling production wells and the inflow source terms $q_{\psi }^{-}\le 0$ devoted to injection wells. In the following, we suppose that we only inject water or gas in the reservoir: $q_O^{-}=0$.

C _h and C _v  = 1 − C _h denote respectively the mass fraction of the heavy and volatile components in the oil phase. The volatile component mass transfer between the oil and gas phases is governed by the equilibrium $X_v=1/K_v$ where K _v is an equilibrium constant depending on pressure. X _v is the molar fraction of the volatile component in the oil phase defined from C _h and C _v and components molar masses.

Under laminar flow conditions, the pure phase velocities in permeable porous media are governed by generalized Darcy laws $${\mathbf{u}}_{\psi }=-\frac{\mathrm{Kk}{r}_{\psi }}{{\mu }_{\psi }}\left(\nabla P_{\psi }-{\rho }_{\psi }\mathbf{g}\right),$$(2) where K is the rock permeability, kr _ψ is the relative permeability for the phase ψ, μ _ψ is the pure phase viscosity, P _ψ is the pressure of the phase ψ and g is the gravity acceleration. In order to simplify notations, we introduce the phase mobility ${\lambda }_{\psi }=k{r}_{\psi }/{\mu }_{\psi }$ and the hydrodynamic gradient $${\mathbf{v}}_{\psi }=-K\left(\nabla P_{\psi }-{\rho }_{\psi }\mathbf{g}\right),$$ such that ${\mathbf{u}}_{\psi }={\lambda }_{\psi }{\mathbf{v}}_{\psi }$.

2.2 Chemical tracer model

We consider a generic chemical tracer model without presuming of the nature of the chemical tracer. As several chemical tracers are involved in the SWCTT method, we consider here a chemical tracer denoted α evolving in a reservoir containing a fluid composed of at least one of the water, oil and gas phases as described by the system (1). When the three phases are present at a given bulk volume, the chemical tracer molecules are possibly partitioned between these three phases. We suppose that this partitioning between phases is instantaneous and that we can distinguish one phase for which the chemical tracer concentration is the highest. This particular phase is the reference phase and is denoted ${\psi }_{\alpha }\in \{W,O,G\}$ for tracer α. This phase has to be present wherever the chemical tracer is present.

The mass conservation equation of the chemical tracer α is converted to a mole conservation equation dividing it by the molar mass M _α . The mole conservation equation is given by $${\mathrm{\partial }}_t\left(n_{\alpha }\right)+\nabla \cdot {\mathbf{F}}_{\alpha }+\nabla \cdot {\mathbf{J}}_{\alpha }+S_{\alpha }^I+S_{\alpha }^P+S_{\alpha }^{\mathrm{Deg}}+S_{\alpha }^{\mathrm{Rea}}=0,$$(3) where n _α is the tracer mole number given by $$n_{\alpha }=\mathrm{\Phi }\sum _{\psi \in \{W,O,G\}}\mathrm{}{S}_{\psi }{C}_{\alpha ,\psi }+(1-\mathrm{\Phi })C_{\alpha ,R},$$(4) where C _α,ψ is the mole concentration of tracer α in the phase ψ and C _α,R the mole concentration of tracer adsorbed on the rock. C _α,R is governed by a Langmuir isotherm which has been extended to the case of a phase partitioning tracer $$C_{\alpha ,R}=C_{\alpha ,R}^{\max }\frac{{\tilde{b}}_{\alpha }\sum _{\psi \in \{W,O,G\}}\mathrm{}{S}_{\psi }{C}_{\alpha ,\psi }}{1+{\tilde{b}}_{\alpha }\sum _{\psi \in \{W,O,G\}}\mathrm{}{S}_{\psi }{C}_{\alpha ,\psi }},$$(5) where $C_{\alpha ,R}^{\max }$ is the maximum adsorbed mole concentration of chemical tracer α on the rock. It is deduced from the adsorption maximum mass fraction $X_{\alpha ,R}^{\max }$ through the relation $C_{\alpha ,R}^{\max }={\rho }_R{X}_{\alpha ,R}^{\max }/M_{\alpha }$ where ρ _R is the rock density. The coefficient ${\tilde{b}}_{\alpha }$ is deduced from Langmuir coefficient b _α (dimensionless) with the relation ${\tilde{b}}_{\alpha }=b_{\alpha }{M}_{\alpha }/{\rho }_{{\psi }_{\alpha }}{S}_{{\psi }_{\alpha }}$.

The advection flux of the chemical tracer is given by $${\mathbf{F}}_{\alpha }=\sum _{\psi \in \{W,O,G\}}\mathrm{}{\mathbf{u}}_{\psi }{C}_{\alpha ,\psi }.$$(6)

The diffusion and dispersion contributions are written as follows in a homogeneous porous medium $${\mathbf{J}}_{\alpha }=-\sum _{\psi \in \{W,O,G\}}\mathrm{}\left(\mathrm{\Phi }{S}_{\psi }\frac{D_{\alpha ,\psi }^{\mathrm{df}}}{\tau }+D_{\psi }^{\mathrm{ds}}\left|{\mathbf{u}}_{\psi }\right|\right)\nabla C_{\alpha ,\psi },$$(7) where $D_{\alpha ,\psi }^{\mathrm{df}}$ is the molecular diffusion constant of the chemical tracer α in the phase ψ, τ is the tortuosity and $D_{\psi }^{\mathrm{ds}}$ is the dispersion coefficient.

In this paper, we suppose that every chemical tracer is injected in the reservoir only through its reference phase such that the injection source term reads $$S_{\alpha }^I=q_{{\psi }_{\alpha }}^{-}{C}_{\alpha ,{\psi }_{\alpha }}^{\mathrm{inj}},$$(8) where we recall that $q_{{\psi }_{\alpha }}^{-}$ is the injected reference phase flow rate per volume and $C_{\alpha ,{\psi }_{\alpha }}^{\mathrm{inj}}$ is the injected mole concentration of chemical tracer α.

The production terms are given by $$S_{\alpha }^P=\sum _{\psi \in \{W,O,G\}}\mathrm{}{q}_{\psi }^{+}{C}_{\alpha ,\psi },$$(9) where we recall that $q_{\psi }^{+}$ is the produced flow rate of phase $\psi$.

The chemical tracer molecules may degrade with time. As a main assumption, we suppose that the degradation of a tracer only occurs in its reference phase such that the degradation effects are written $$S_{\alpha }^{\mathrm{Deg}}=\mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha }^{\mathrm{Deg}}{C}_{\alpha ,{\psi }_{\alpha }},$$(10) where ${\lambda }_{\alpha }^{\mathrm{Deg}}$ is the kinetic constant of tracer $\alpha$. In practice, the half life time $t_{\alpha ,1/2}^{\mathrm{Deg}}$ is commonly used to characterize this kinetics and is related to the former constant with the relation $t_{\alpha ,1/2}^{\mathrm{Deg}}=\ln 2/{\lambda }_{\alpha }^{\mathrm{Deg}}$.

A classical SWCTT generally involves a single chemical reaction. Let us then suppose that α and $\mathcal{R}$ are the two reactants of this chemical reaction, with α the primary tracer to be followed. Reaction of α with $\mathcal{R}$ generates chemical tracer β and other products denoted $\mathcal{P}$. This type of reaction reads $$\left(\mathrm{Rea}\right) \mathcal{R}+\alpha \to \mathcal{P}+{\gamma }_{\alpha ,\beta }\beta ,$$(11) where ${\gamma }_{\alpha ,\beta }$ is the stoechiometric coefficient of tracer β with respect to tracer α. We suppose also that the reactant $\mathcal{R}$ is present in very large quantity compared to the reaction consumption quantities (usually $\mathcal{R}$ is the reference phase molecule, namely H₂O for an aqueous carrier phase) such that there is no need to model the reactant consumption in our model. In addition, the products $\mathcal{P}$ do not need to be modeled because they have no effect on the flow of phases. In that context where it is useless to model $\mathcal{R}$ and $\mathcal{P}$, this type of reaction simplifies as $$\left(\mathrm{Rea}\right) \alpha \to {\gamma }_{\alpha ,\beta }\beta .$$(12)

This chemical reaction has its own kinetics characterized by a half life time that we will denote $t_{\alpha \beta ,1/2}^{\mathrm{Rea}}$ or the kinetic constant ${\lambda }_{\alpha \beta }^{\mathrm{Rea}}=\ln 2/t_{\alpha \beta ,1/2}^{\mathrm{Rea}}$. Herein, the reaction rate is assumed invariant, although that may not be valid in some cases. Indeed, the rate of reaction may change if temperature, salinity or pH conditions, that are reaction-rate-impacting variables, change significantly in the volume investigated by the test and/or in the course of the test [13]. Thus, the physical model corresponding to this kinetic chemical reaction for the chemical tracers $\alpha$ and $\beta$ writes as follows: $$\{\begin{array}{l}{S}_{\alpha }^{\mathrm{Rea}}=\mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha \beta }^{\mathrm{Rea}}{C}_{\alpha ,{\psi }_{\alpha }},\\ S_{\beta }^{\mathrm{Rea}}=-\mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha \beta }^{\mathrm{Rea}}{\gamma }_{\alpha ,\beta }{C}_{\alpha ,{\psi }_{\alpha }}.\end{array}$$(13)

This chemical reaction couples the mole conservation equation of chemical tracers $\alpha$ and $\beta$ via $S_{\alpha }^{\mathrm{Rea}}$ and $S_{\beta }^{\mathrm{Rea}}$ terms. As a consequence, it will not be possible to perform their numerical resolutions separately.

Considering the phase transfer mechanisms of the reactive chemical tracer α, its mole concentrations in the water, oil and gas phases are linked to the mole concentration in the reference phase by the partition coefficients $K_{\alpha ,\psi -{\psi }_{\alpha }}$ as follows $$C_{\alpha ,\psi }=K_{\alpha ,\psi -{\psi }_{\alpha }}{C}_{\alpha ,{\psi }_{\alpha }}\hspace{1em}\forall \psi \in \{W,O,G\}.$$(14)

Using this relation, the chemical tracer mole conservation equation can be formulated with its mole concentration in the reference phase as unknown $$\{\begin{array}{lll}& & {\mathrm{\partial }}_t\left(\mathrm{\Phi }\sum _{\psi =W,O,G}\mathrm{}{S}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\right)\\ & +& {\mathrm{\partial }}_t\left((1-\mathrm{\Phi })C_{\alpha ,R}^{\max }\frac{{\tilde{b}}_{\alpha }\sum _{\psi =W,O,G}\mathrm{}{S}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}}{1+{\tilde{b}}_{\alpha }\sum _{\psi =W,O,G}\mathrm{}{S}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}}\right)\\ & +& \nabla \cdot \left(\sum _{\psi =W,O,G}\mathrm{}{\mathbf{u}}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\right)\\ & -& \nabla \cdot \left(\sum _{\psi =W,O,G}\mathrm{}\left(\mathrm{\Phi }{S}_{\psi }\frac{D_{\alpha ,\psi }^{\mathrm{df}}}{\tau }+D_{\psi }^{\mathrm{ds}}\left|{\mathbf{u}}_{\psi }\right|\right)\nabla \left(K_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\right)\right)\\ & +& q_{{\psi }_{\alpha }}^{-}{C}_{\alpha ,{\psi }_{\alpha }}^{\mathrm{inj}}+\sum _{\psi =W,G}\mathrm{}{q}_{\psi }^{+}{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\\ & +& \mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha }^{\mathrm{Deg}}{C}_{\alpha ,{\psi }_{\alpha }}+\mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha \beta }^{\mathrm{Rea}}{C}_{\alpha ,{\psi }_{\alpha }}=0.\\ & & \end{array}$$(15)

The system under this form does not exhibit any singularity. The good symmetry of the system with respect to the present phases simplifies the implementation of the proposed numerical schemes.

3 Numerical scheme

The chemical tracer has no significant effect on the fluid present in the reservoir, namely water, oil and gas. Thus the resolution of chemical tracers mole conservation equations is decoupled from the resolution of phases mass conservation equation. The benefit is that different types of scheme can be easily considered for the tracers without underlying mathematical difficulties. In particular, we consider an implicit scheme for the phases and an explicit scheme for the chemical tracers. In the first subsection, we detail the notations used for the mesh discretization. In the second subsection, we quickly describe the water oil and gas phases mass conservation discretization which is classical. In the third subsection, after describing the first-order explicit discretization of chemical tracers mole conservation equations, we describe its extension to a second-order scheme in time and space.

3.1 Domain discretization

Let ${\mathcal{M}}_h$ be an admissible finite volume mesh of the reservoir given by a family of control volumes or cells noted K. For any K of ${\mathcal{M}}_h$, |K| is its measure, $\mathrm{KL}=\mathrm{\partial }K\cap \mathrm{\partial }L$ is the common interface between K and a neighbouring cell L, and n _KL denotes the unit normal vector to the interface KL oriented from K to L. The set of neighbors of the cell K is denoted $\mathcal{N}\left(K\right)$, that is $\mathcal{N}\left(K\right)=\{L\in {\mathcal{M}}_h;\mathrm{\partial }K\cap \mathrm{\partial }L\ne \mathrm{\varnothing }\}$.

To construct our second-order in space scheme, we will need to define the center ${\mathcal{M}}_K$ of each cell K. For each neighbor cell $L\in \mathcal{N}\left(K\right)$, we denote d _K,KL the distance of the point ${\mathcal{M}}_K$ to the interface KL and d _L,KL the distance of the point M _L to the interface KL.

Let Π_prod be the set of perforations of the production wells and Π_inj be the set of perforations of the injection wells and M the set of all the perforated cells in ${\mathcal{M}}_h$. For instance, K is a perforated cell if it exists a perforation $j\in {\mathrm{\Pi }}_{\mathrm{prod}}$ or $j\in {\mathrm{\Pi }}_{\mathrm{inj}}$ such that M(j) = K. In the next subsections, we note ${\mathrm{\Pi }}_p\left(K\right)=\{j\in {\mathrm{\Pi }}_{\mathrm{prod}}|M\left(j\right)=K\}$ and ${\mathrm{\Pi }}_i\left(K\right)=\{j\in {\mathrm{\Pi }}_{\mathrm{inj}}|M\left(j\right)=K\}$ for concision.

Finally, we introduce an increasing sequence of discrete times $\{t^n{\}}_{0\le n\le N}$ such that t ⁰ = 0 and t ^N  = T. Then the time interval (0, T) is subdivided into N variable time steps Δt = t ⁿ⁺¹ − t ⁿ , 0 ≤ n ≤ N − 1.

3.2 Water and hydrocarbons scheme

In the context of a SWCTT, the oil and gas quantities remaining in the reservoir are supposed to be low as soon as the chemical tracers are injected. This is especially the case for the oil as its saturation has to be close to the residual oil saturation to avoid its movement during the SWCTT. Thus, the injection and production of the aqueous carrier phase of the tracer do not change the oil and water saturations within the reservoir. In this context, there is no need to consider a high-order scheme for the water and hydrocarbon components mass conservation equations. In fact, the main requirement is only to dispose of a robust and unconditionally-stable resolution scheme. Therefore, we propose to use an implicit scheme in time and an upwind finite volume method for the spatial derivatives $$\{\begin{array}{l}\left|K\right|\frac{{\left(\mathrm{\Phi }{\rho }_O{C}_h{S}_O\right)}_K^{n+1}-{\left(\mathrm{\Phi }{\rho }_O{C}_h{S}_O\right)}_K^n}{\Delta t}\\ +\sum _{L\in \mathcal{N}\left(K\right)}\mathrm{}\left({\rho }_O{\lambda }_O{C}_h{)}_K^{n+1}\right({\mathbf{v}}_{O,\mathrm{KL}}^{n+1}{)}^{+}+\left({\rho }_O{\lambda }_O{C}_h{)}_L^{n+1}\right({\mathbf{v}}_{O,\mathrm{KL}}^{n+1}{)}^{-}\\ +\sum _{j\in {\mathrm{\Pi }}_p\left(K\right)}\mathrm{}\left(C_h{)}_K^{n+1}\right({\rho }_O{)}_K^{n+1}(q_O^{+}{)}_j^{n+1}=0,\\ \\ \left|K\right|\frac{{\left(\mathrm{\Phi }{\rho }_O{C}_v{S}_O+\mathrm{\Phi }{\rho }_G{S}_G\right)}_K^{n+1}-{\left(\mathrm{\Phi }{\rho }_O{C}_v{S}_O+\mathrm{\Phi }{\rho }_G{S}_G\right)}_K^n}{\Delta t}\\ +\sum _{L\in \mathcal{N}\left(K\right)}\mathrm{}\left({\rho }_O{\lambda }_O{C}_v{)}_K^{n+1}\right({\mathbf{v}}_{O,\mathrm{KL}}^{n+1}{)}^{+}+\left({\rho }_O{\lambda }_O{C}_v{)}_L^{n+1}\right({\mathbf{v}}_{O,\mathrm{KL}}^{n+1}{)}^{-}\\ +\sum _{L\in \mathcal{N}\left(K\right)}\mathrm{}\left({\rho }_G{\lambda }_G{)}_K^{n+1}\right({\mathbf{v}}_{G,\mathrm{KL}}^{n+1}{)}^{+}+\left({\rho }_G{\lambda }_G{C}_v{)}_L^{n+1}\right({\mathbf{v}}_{G,\mathrm{KL}}^{n+1}{)}^{-}\\ +\sum _{j\in {\mathrm{\Pi }}_p\left(K\right)}\mathrm{}\left({\rho }_G{)}_K^{n+1}\right(q_G^{+}{)}_j^{n+1}+\left(C_v{\rho }_O{)}_K^{n+1}\right(q_O^{+}{)}_j^{n+1}\\ +\sum _{j\in {\mathrm{\Pi }}_i\left(K\right)}\mathrm{}\left({\rho }_G{)}_K^{n+1}\right(q_G^{-}{)}_j^{n+1}=0,\\ \\ \left|K\right|\frac{{\left(\mathrm{\Phi }{\rho }_W{S}_W\right)}_K^{n+1}-{\left(\mathrm{\Phi }{\rho }_W{S}_W\right)}_K^n}{\Delta t}\\ +\sum _{L\in \mathcal{N}\left(K\right)}\mathrm{}\left({\rho }_W{\lambda }_W{)}_K^{n+1}\right({\mathbf{v}}_{W,\mathrm{KL}}^{n+1}{)}^{+}+\left({\rho }_W{\lambda }_W{)}_L^{n+1}\right({\mathbf{v}}_{W,\mathrm{KL}}^{n+1}{)}^{-}\\ +\sum _{j\in {\mathrm{\Pi }}_p\left(K\right)}\mathrm{}\left({\rho }_W{)}_K^{n+1}\right(q_W^{+}{)}_j^{n+1}\\ +\sum _{j\in {\mathrm{\Pi }}_i\left(K\right)}\mathrm{}\left({\rho }_W{)}_K^{n+1}\right(q_W^{-}{)}_j^{n+1}=0,\end{array}$$(16) where ${\mathbf{v}}_{\psi ,\mathrm{KL}}^{n+1}$ is a finite volume discretization of the flux ${\int }_{\mathrm{KL}}\mathrm{}{\mathbf{v}}_{\psi }\cdot {\mathbf{n}}_{\mathrm{KL}}\mathrm{d}\sigma$ at the time t ⁿ⁺¹. For each phase, we define the operator $$({\mathbf{v}}_{\psi ,\mathrm{KL}}{)}^{+}=\max (0,{\mathbf{v}}_{\psi ,\mathrm{KL}}) \mathrm{and} ({\mathbf{v}}_{\psi ,\mathrm{KL}}{)}^{-}=\min (0,{\mathbf{v}}_{\psi ,\mathrm{KL}}),$$ in order to have the upwind scheme. This numerical scheme is classical and addressed in the literature [12].

3.3 Chemical tracer scheme

In this work, a second-order numerical scheme is implemented because it is robust and more accurate than a first-order scheme to solve the chemical model of a SWCTT. Furthermore, an explicit scheme is adopted to limit numerical diffusion. Explicit scheme remains computationally acceptable for models comprising a limited number of cells, such as the near-wellbore models used to simulate SWCTT. The CFL criterion of explicit scheme was not a major issue in the numerical experiments of this paper.

In this section, for the sake of briefness, the system (15) is simplified by omitting the diffusion-dispersion terms $$\{\begin{array}{lll}& & {\mathrm{\partial }}_t\left(\mathrm{\Phi }\sum _{\psi =W,O,G}\mathrm{}{S}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\right)\\ & +& {\mathrm{\partial }}_t\left((1-\mathrm{\Phi })C_{\alpha ,R}^{\max }\frac{{\tilde{b}}_{\alpha }\sum _{\psi =W,O,G}\mathrm{}{S}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}}{1+{\tilde{b}}_{\alpha }\sum _{\psi =W,O,G}\mathrm{}{S}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}}\right)\\ & +& \nabla \cdot \left(\sum _{\psi =W,O,G}\mathrm{}{\mathbf{u}}_{\psi }{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\right)\\ & +& q_{{\psi }_{\alpha }}^{-}{C}_{\alpha ,{\psi }_{\alpha }}^{\mathrm{inj}}+\sum _{\psi =W,G}\mathrm{}{q}_{\psi }^{+}{K}_{\alpha ,{\psi }_{\alpha }-\psi }{C}_{\alpha ,{\psi }_{\alpha }}\\ & +& \mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha \beta }^{\mathrm{Deg}}{C}_{\alpha ,{\psi }_{\alpha }}+\mathrm{\Phi }{S}_{{\psi }_{\alpha }}{\lambda }_{\alpha \beta }^{\mathrm{Rea}}{C}_{\alpha ,{\psi }_{\alpha }}=0.\\ & & \end{array}$$(17)

In the next subsection, we firstly describe the chemical tracer basic scheme that is first-order in space and time. This scheme is extended in Section 3.3.2 to a second-order scheme in time by using a predictor-corrector scheme. Then, in Section 3.3.3, the second order scheme in space via a MUSCL technique [16– 20] is also considered because monotonous and stable [24, 25].

3.3.1 Chemical tracer first-order scheme

We propose an explicit scheme in time and upwind finite-volume approximations for the spatial derivatives, as follows: $$\{\begin{array}{lll}{n}_K^{n+1}& =& n_K^n\\ & -& \frac{\Delta t}{\left|K\right|}\sum _{L\in \mathcal{N}\left(K\right)}\mathrm{}\sum _{\psi =W,O,G}\mathrm{}\left({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right({\lambda }_{\psi }{)}_K^{n+1}({\mathbf{v}}_{\psi ,\mathrm{KL}}^{n+1}{)}^{+}\\ & & +\left({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right({\lambda }_{\psi }{)}_L^{n+1}({\mathbf{v}}_{\psi ,\mathrm{KL}}^{n+1}{)}^{-}\\ & -& \frac{\Delta t}{\left|K\right|}\sum _{j\in {\mathrm{\Pi }}_p\left(K\right)}\mathrm{}\sum _{\psi =W,O,G}\mathrm{}\left(q_{\psi }^{+}{)}_j^{n+1}{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n\\ & -& \frac{\Delta t}{\left|K\right|}\sum _{j\in {\mathrm{\Pi }}_i\left(K\right)}\mathrm{}\left(q_{{\psi }_{\alpha }}^{-}{)}_j^{n+1}\right(C_{\alpha ,{\psi }_{\alpha }}{)}_j^{n+1}\\ & -& \left(\mathrm{\Phi }{S}_{{\psi }_{\alpha }}{)}_K^{n+1}\right(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n\left(e^{-{\lambda }_{\alpha }^{\mathrm{Deg}}\Delta t}+e^{-{\lambda }_{\alpha \beta }^{\mathrm{Rea}}\Delta t}\right),\\ & & \end{array}$$(18) where we denote $n_K^n$ and $n_K^{n+1}$ the total chemical tracer mole per unit volume respectively at times t ⁿ and t ⁿ⁺¹. They are defined as the sum of the mobile moles present in the three phases and moles adsorbed on the rock in the cell. For each timestep l with 0 ≤ l ≤ N − 1 corresponding to time t ^l , $n_K^l$ has to be shared between the phases and the rock surface by solving the mass conservation equation of cell tracer with adsorption isotherm as follows $$\begin{array}{ll}{n}_K^l& =\mathcal{A}\left((C_{\alpha ,{\psi }_{\alpha }}{)}_K^l;{\mathrm{\Phi }}_K^0,{\mathrm{\Phi }}_K^l,(S_{\psi }{)}_K^l\right)\\ & ={\mathrm{\Phi }}_K^l\left(\sum _{\psi =W,O,G}\mathrm{}(S_{\psi }{)}_K^l{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right)(C_{\alpha ,{\psi }_{\alpha }}{)}_K^l\\ & +(1-{\mathrm{\Phi }}_K^0)C_{\alpha ,R}^{\max }\frac{\left({\tilde{b}}_{\alpha }{)}_K^l\left(\sum _{\psi =W,O,G}\mathrm{}(S_{\psi }{)}_K^l{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right)\right(C_{\alpha ,{\psi }_{\alpha }}{)}_K^l}{1+\left({\tilde{b}}_{\alpha }{)}_K^l\left(\sum _{\psi =W,O,G}\mathrm{}(S_{\psi }{)}_K^l{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right)\right(C_{\alpha ,{\psi }_{\alpha }}{)}_K^l},\\ & \end{array}$$(19) where $\mathcal{A}$ denotes the adsorption operator and $({\tilde{b}}_{\alpha }{)}_K^l=b_{\alpha }{M}_{\alpha }/({\rho }_{{\psi }_{\alpha }}{S}_{{\psi }_{\alpha }}{)}_K^l$. Then, the molar concentration $(C_{\alpha ,{\psi }_{\alpha }}{)}_K^{n+1}$ is computed from the total mole number of chemical tracer $n_K^{n+1}$ in the cell by inverting the previous relation with l = n + 1. This inversion leads to the resolution of a second-order polynomial with the unknown $(C_{\alpha ,{\psi }_{\alpha }}{)}_K^l$ (with the previous notation) which admits only one positive root $$\begin{array}{ll}(C_{\alpha ,{\psi }_{\alpha }}{)}_K^l& ={\mathcal{A}}^{-1}\left((n_{\alpha }{)}_K^l;{\mathrm{\Phi }}_K^0,{\mathrm{\Phi }}_K^l,(S_{\psi }{)}_K^l\right)\\ & =\frac{\sqrt{4(\mathrm{\Phi }{n}_{\alpha }{\tilde{b}}_{\alpha }{)}_K^l+{\theta }^2}-\theta }{2(\mathrm{\Phi }{\tilde{b}}_{\alpha }{)}_K^l\left(\sum _{\psi =W,O,G}\mathrm{}(S_{\psi }{)}_K^l{K}_{\alpha ,{\psi }_{\alpha }-\psi }\right)},\end{array}$$(20) with $$\theta ={\mathrm{\Phi }}_K^l+({\tilde{b}}_{\alpha }{)}_K^l\left((1-{\mathrm{\Phi }}_K^0)C_{\alpha ,R}^{\max }-(n_{\alpha }{)}_K^l\right).$$(21)

The values $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n$ and $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n$ denote respectively the mole concentration in the cell $(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n$ and $(C_{\alpha ,{\psi }_{\alpha }}{)}_L^n$ if a first-order scheme in space is considered.

In order to simplify the notation when describing how to extend our numerical scheme to second-order in time and space, we propose to rewrite the system (18) in a more concise form to obtain $n_K^{n+1}$ from $n_K^n$ as follows $$n_K^{n+1}=n_K^n+\mathcal{F}\left(\Delta t,({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n,({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n,(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n\right).$$(22)

This concise form does not explicitly detail the porosity, the phases densities, the saturations, the velocities, the injection and production source terms needed in the function $\mathcal{F}$ because they are taken at time t ⁿ⁺¹ from implicit resolution of (16). In addition, as only a first-order scheme is considered in this subsection, we set $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n=(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n$ and $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n=(C_{\alpha ,{\psi }_{\alpha }}{)}_L^n$.

For the advection, injection and production terms of this numerical scheme, a CFL criterion can de derived. For this, we extend the work of Braconnier et al. [4] based on a von Neumann analysis to our more complex scheme to obtain $$\Delta t\le \underset{\alpha }{\min }\left(\underset{K\in {\mathcal{M}}_h}{\min }\left(\frac{1}{({\lambda }_{\alpha }^{\mathrm{adv}}{)}_K^n+({\lambda }_{\alpha }^{\mathrm{well}}{)}_K^n}\right)\right),$$(23) with $$({\lambda }_{\alpha }^{\mathrm{adv}}{)}_K^n=\sum _{L\in \mathcal{N}\left(K\right)}\mathrm{}\sum _{\psi =W,O,G}\mathrm{}{K}_{\alpha ,{\psi }_{\alpha }-\psi }({\lambda }_{\psi }{)}_K^{n+1}(v_{\psi ,\mathrm{KL}}^{n+1}{)}^{+},$$(24) and $$({\lambda }_{\alpha }^{\mathrm{well}}{)}_K^n=\sum _{j\in {\mathrm{\Pi }}_p\left(K\right)}\mathrm{}\sum _{\psi =W,O,G}\mathrm{}(q_{\psi }^{+}{)}_j^{n+1}{K}_{\alpha ,{\psi }_{\alpha }-\psi }.$$(25)

Application of the above first-order scheme in Section 4 will be performed with an approximate CFL criterion to ensure stability.

3.3.2 Second-order scheme in time for the chemical tracers

The numerical method used in this paper to build a second-order method is classical [14, 15]. Considering formula (22), second-order scheme is written as $$\{\begin{array}{ll}{n}_K^{[n+1]}& =n_K^n\\ & +\mathcal{F}\left(h_1\Delta t,({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n,({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n,(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n\right),\\ & \\ n_K^{n+1}& =a_1{n}_K^n+a_2{n}_K^{[n+1]}\\ & +\mathcal{F}\left(h_2\Delta t,({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^{[n+1]},({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^{[n+1]},(C_{\alpha ,{\psi }_{\alpha }}{)}_K^{[n+1]}\right).\end{array}$$(26)

This scheme involves the computation of an intermediate solution $n_K^{[n+1]}$. The coefficient a ₁, a ₂, h ₁ and h ₂ are set to different values, depending on the integration method choosen. Some classical values, corresponding to classical schemes [14, 15], are given in Table 1.

Even if switching from one scheme to another is simple, only the predictor-corrector scheme will be used in the simulation Section 4 for the sake of briefness of this paper.

In the framework of a first-order scheme in space, the slopes of the chemical tracer mole concentration are not computed such that $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n=(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n$, $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n=(C_{\alpha ,{\psi }_{\alpha }}{)}_L^n$, $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^{[n+1]}=(C_{\alpha ,{\psi }_{\alpha }}{)}_K^{[n+1]}$ and $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^{[n+1]}=(C_{\alpha ,{\psi }_{\alpha }}{)}_L^{[n+1]}$. The next subsection describes a method to further improve the resolution by including also a second-order scheme in space, based on the MUSCL method.

3.3.3 Second-order scheme in space for the chemical tracers

As said in the former section of this paper, we focus on a method to improve the S _orw assessment from numerical simulations. As we use a simple upwind and explicit scheme for the chemical tracers resolution (22) and (26), one way to improve easily the solution accuracy is achieved via the MUSCL method. One important feature of the MUSCL method is that the positivity of the reconstructed unknowns can easily be preserved with slope limiters such as Minmod or Superbee. Thus, using this method, it is possible to improve the accuracy of chemical tracers concentration peaks.

Our MUSCL technique is based on slope reconstruction to improve the molar concentration needed in the fluxes of the schemes (22) and (26). The states $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^l$ and $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^l$ are reconstucted by using the states $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{K\mathrm{\prime }}^l$ with $K\mathrm{\prime }\in {\mathcal{M}}_h$. First, let us first determine the gradient of the chemical tracer molar concentration across each interface $\mathrm{KL}$ that separates the cell $K$ and the cell $L\in \mathcal{N}\left(K\right)$ $${\nabla }_{\mathrm{KL}}{C}_{\alpha ,{\psi }_{\alpha }}^n=\frac{(C_{\alpha ,{\psi }_{\alpha }}{)}_L^n-(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n}{d_{K,\mathrm{KL}}+d_{L,\mathrm{KL}}}{\mathbf{n}}_{\mathrm{KL}}.$$(27)

In the context of the MUSCL technique, we have to search the cells L ^dw which position is downward to the cells K and L in the direction defined by the cell interface normal n _KL . From a mathematical point of view, we search $L^{\mathrm{dw}}\in \mathcal{N}\left(K\right)$ such that ${\mathbf{n}}_{\mathrm{KL}}\cdot {\mathbf{n}}_{K{L}^{\mathrm{dw}}}$ is strictly negative and minimum. If one such cell exists, we compute the gradient $${\nabla }_{K{L}^{\mathrm{dw}}}{C}_{\alpha ,{\psi }_{\alpha }}^n=\frac{(C_{\alpha ,{\psi }_{\alpha }}{)}_{L^{\mathrm{dw}}}^n-(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n}{d_{K,K{L}^{\mathrm{dw}}}+d_{L^{\mathrm{dw}},K{L}^{\mathrm{dw}}}}{\mathbf{n}}_{K{L}^{\mathrm{dw}}},$$(28) and then compute a limited slope $${\tilde{\nabla }}_{\mathrm{KL}}{C}_{\alpha ,{\psi }_{\alpha }}^n\cdot {\mathbf{n}}_{\mathrm{KL}}=\mathcal{L}({\nabla }_{\mathrm{KL}}{C}_{\alpha ,{\psi }_{\alpha }}^n\cdot {\mathbf{n}}_{\mathrm{KL}},{\nabla }_{K{L}^{\mathrm{dw}}}{C}_{\alpha ,{\psi }_{\alpha }}^n\cdot {\mathbf{n}}_{K{L}^{\mathrm{dw}}}),$$(29) where the function $\mathcal{L}$ is a slope limiter for which exemples are given below. If there is no downwind cell L ^dw to the cell K, we simply set ${\tilde{\nabla }}_{\mathrm{KL}}{C}_{\alpha ,{\psi }_{\alpha }}^n=0$. The case where several downwind cells L ^dw to the cell K is not considered here because we limited the simulation tests to a regular mesh in the frame of this paper.

Then we reconstruct a second-order approximate molar concentration at the interface between the cells K and L as $$({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{KL}}^n=(C_{\alpha ,{\psi }_{\alpha }}{)}_K^n+d_{K,\mathrm{KL}}{\tilde{\nabla }}_{\mathrm{KL}}{C}_{\alpha ,{\psi }_{\alpha }}^n\cdot {\mathbf{n}}_{\mathrm{KL}},$$(30) to feed the relation (22) or (26). A similar approach is applied to cell L to find its upwind cell relatively to the interface KL, to determine the limited slope ${\tilde{\nabla }}_{\mathrm{LK}}{C}_{\alpha ,{\psi }_{\alpha }}^n\cdot {\mathbf{n}}_{\mathrm{KL}}$ and to compute $({\tilde{C}}_{\alpha ,{\psi }_{\alpha }}{)}_{\mathrm{LK}}^n$.

By using these reconstructed molar concentrations, schemes (22) or (26) become second-order in space.

As slope limiter, we can choose the Minmod limiter defined as $$\mathcal{L}(x,y)=\frac{1}{2}\left(\frac{x}{\left|x\right|}+\frac{y}{\left|y\right|}\right)\min \left(\right|x|,|y\left|\right),$$(31) or the Superbee limiter $$\mathcal{L}(x,y)=\frac{1}{2}\left(\frac{x}{\left|x\right|}+\frac{y}{\left|y\right|}\right)\max \left(\min \left(2\right|x|,|y\left|\right),\min \left(\right|x|,2|y\left|\right)\right).$$(32)

As mentioned in the work of Berthon [19, 24, 25], we can demonstrate the stability and robustness of such a space scheme, whether a first-order or a second-order time scheme is used.

4 Numerical experiments

4.1 Case study description

We consider a part of a real reservoir simulation model which is centered around a vertical testing well. Although the chemical tracer slugs propagate only 2–10 m from wellbore, the lateral extension of that model grid was chosen large enough to allow the use of closed (no flow) boundary condition, without any impact on SWCTT simulation results. Such a model remained computationally affordable. Thus, we consider a square bounded domain of the reservoir that is centered around the well with a lateral extension of about 200 m. This reservoir is vertically layered into 18 layers whose porosity, horizontal and vertical permeabilities are constant. Porosity is set to 10% in every layer whereas the permeability differs from one layer to another as reported in Table 2. As a rule of thumb, the vertical permeability was set to one tenth of the horizontal permeability in each layer, that is $K_V/K_H=0.1$.

The considered well test is supposed to be an isothermal two-phase flow process involving water and oil phases, whose transport is described by the generalized Darcy equations, as explained in Section 2. Capillary pressures are neglected and relative permeabilities to water and oil are assumed to be power laws of the water saturation and read $k{r}_W\left(S\right)=k{r}_W^0{S}^{n_W}$ and $k{r}_O\left(S\right)=k{r}_O^0(1-S{)}^{n_O}$, where S = (S _W  − S _wi )/(1 − S _orw  − S _wi ) denotes the normalized (mobile) water saturation, S _wi  = 0.1 and S _orw  = 0.2 are the (immobile) irreducible water saturation and residual oil saturation to waterflooding, $k{r}_W^0=0.2$ and $k{r}_O^0=0.9$ are the maximum relative permeabilites to water and oil, and n _W  = n _O  = 2 their exponents.

The well test schedule and fluids composition are given in Table 3: prior to the injection steps reported in Table 3, the reservoir is initially saturated with oil such that S _O (r, t = 0) = 1 − S _wi . Starting from this initial condition, the reservoir is waterflooded for 10 days with a constant flow rate $q_W^{-}=150$ m³/day at bottom conditions. This waterflood leads to a homogeneous residual oil saturation state S _O (r,t) = S _orw in the near wellbore region for r ≲ 10–20 m. Then, an inert tracer, denoted $t$, and an ester, denoted e, are injected for $0.5$ day with the same flow rate and aqueous mass fractions of $C_{t,W}^{\mathrm{inj}}=C_{e,W}^{\mathrm{inj}}=10^3$ ppm, followed by a 1.5 days water drive and a shut-in period of 3 days, that aims at generating a significant amount of alcohol, denoted a, from the ester which is at rest about 5 m away from the injection well. The flow is then reverted for 5 days with a constant liquid flow rate $q_W^{+}+q_O^{+}=150$ m³/day so that the alcohol, inert tracer and ester are recovered. As an illustration, Figure 11 reports the spatial distribution of ester, alcohol and tracer concentrations, after the ester placement at t = 12 days, after the shut-in period at t = 15 days, and during the recovery period at t = 15.5 days.

It has been shown by many authors, under specific assumptions, that the residual oil saturation can be computed from the recovered ester and alcohol concentration profiles as [1, 2, 26, 27] $$S_{\mathrm{orw}}=\frac{t_e-t_a}{t_e-t_a+K_{e,W-O}(t_a-t_0)},$$(33) where t _a and t _e denote the alcohol and ester arrival times, respectively; t ₀ = 15 days denotes the backflow starting time and $K_{e,W-O}=5$ is the ester partition coefficient, defined as the ratio of the ester molar concentration in the oil phase (expressed in mol/m³) over the ester molar concentration in the water phase, as explained in Section 2.

It is worth noting that the ester, that is initially transported in the water phase when being injected in the porous medium, is assumed to:

1. Instantaneously partition between the mobile water phase and the residual immobile oil phase according to the constant partition coefficient K _e,W−O .

2. Undergo a slow hydrolysis that yields alcohol according to the first-order kinetics $C_{e,W}\left(t\right)=C_{e,W}^{\mathrm{inj}}\exp (-{\lambda }_{\mathrm{ea}}^{\mathrm{Rea}}t)$ where ${\lambda }_{\mathrm{ea}}^{\mathrm{Rea}}$ is the kinetic constant of ester hydrolysis producing alcohol.

The corresponding half life time $t_{\mathrm{ea},1/2}^{\mathrm{Rea}}(=\ln (2)/{\lambda }_{\mathrm{ea}}^{\mathrm{Rea}})$ has been set to 3 days so that the alcohol is mainly generated during the shut-in period, whose duration (3 days) is however not too long so that a non-negligible amount of ester remains at the end of the sequence and can be recovered at wellbore. The duration of tracer placement (0.5 + 1.5 days) was also short to minimize ester hydrolysis during these injection steps; however, that duration might have been reduced further to allow us neglecting ester reaction during placement, as will be shown later on. Needless to say, such a well test can be carefully designed and optimized in several ways depending on the radius of investigation of the ester that is aimed at during the injection process, and according to the time that one can afford in the field. To finish with, point (ii) highlights a key feature of the ester transport and of the SWCTT process: since the ester undergoes instantaneous partition between the mobile water phase and the residual immobile oil phase, or contrasted behavior between its placement that is piston-type, and its recovery that is delayed compared to the alcohol one, one expects its displacement to be caterpillar-like and its arrival time to be late compared to the alcohol one, because the latter is only transported in water.

In the following, we perform simulations of the SWCTT scenario described above, with a first-order scheme and with a second-order scheme both in time and space. For the sake of simplicity, we will denote hereafter the first-order scheme O(1), the second-order scheme with a minmod slope limiter O(2^m), and the second-order scheme with a superbee slope limiter O(2^s).

4.2 Simulations on a three-dimensional regular grid

We consider a three-dimensional regular grid composed of (N _x  = 189) × (N _y  = 189) × (N _z  = 18) = 642 978 cells. The testing well is perforated along all the reservoir height and is located in the center of each x-y layer. We consider the simulation model defined in the previous Section 4.1 and the process schedule reported in Table 3. The goal is to prove the benefit of a second-order resolution of the chemical tracer equation (17) and its relevance to achieve three-dimensional realistic simulations. Simulations were performed in order to compare the O(1), O(2^m) and O(2^s) schemes.

Figures 1– 3 report the simulated ester, alcohol and tracer in situ concentration profiles for three key dates according to the well test schedule of Table 3:

• after the ester placement at t = 12 days,

• after the shut-in period at t = 15 days,

• and finally during the recovery period at t = 15.5 days.

Specifically, Figure 1 compares the inert tracer (simply denoted tracer), ester and alcohol in situ concentration profiles obtained with the O(1), O(2^m) and O(2^s) schemes after the ester placement. These three components are correctly displaced at about 5 m away from the well. The inert tracer goes farther than the ester because the ester undergoes partition in the residual immobile oil phase, as previously explained. A small amount of alcohol is generated from the ester during the forward flow process whereas a significant amount of alcohol is generated during the shut-in period that lasts longer. O(2^s) and O(2^m) schemes yield sharper profiles than the O(1) scheme, the O(2^s) scheme profiles being a bit sharper than the O(2^m) scheme ones.

Figure 2 compares the simulated tracer, ester and alcohol in situ concentration profiles obtained with the O(1), O(2^m) and O(2^s) schemes after the shut-in period. Whereas the inert tracer profile remains the same as in Figure 1, a very significant amount of alcohol has been generated from the ester during the shut-in period. Figure 3 compares the simulated tracer, ester and alcohol in situ concentrations obtained with the O(1), O(2^m) and O(2^s) schemes during the recovery period, and shows in a more pronounced fashion than in Figure 1 at early times that O(2^s) and O(2^m) schemes yield sharper profiles than the O(1) scheme, the O(2^s) scheme profiles being a bit sharper than the O(2^m) scheme ones.

That diffuse behavior of schemes is clearly illustrated by Figure 4, that compares the simulated tracer, ester and alcohol outlet concentration profiles obtained with the O(1), O(2^m) and O(2^s) schemes. Specifically, the O(2^s) scheme reproduces a plateau in the alcohol profile tail which is not visible with the other schemes unless using finer meshes as it will be shown in the next section ( Fig. 5).

Table 4 reports the tracer, ester and alcohol arrival times that are taken at maximum concentration of the peaks for each scheme (peaks are shown with symbols in Fig. 4), and the residual oil saturation S _orw that is computed according to equation (33). It is worth noting that this S _orw assessment is better with the O(2^s) scheme that yields less than $3$% relative error, whereas the O(2^m) and O(1) schemes yield a 12% and 18% error, respectively. These results show that the numerical diffusion clearly degrades the SWCTT simulation results and thus the S _orw assessment. The second-order scheme both in time and space with the superbee limiter gives satisfying results as its numerical diffusion does not spoil the simulation and its results. In addition, the simulation times are of the same order; compared to the O(1) scheme, the O(2^m) and O(2^s) schemes simulation times are about 6% higher which is little considering the results improvement. It proves the benefits of a second-order scheme for chemical tracer modeling and its ability to simulate a realistic three-dimensional SWCTT process. Nevertheless, the fact that the O(2^s) scheme clearly does not preserve the concentration profiles symmetry compared to the O(1) and O(2^m) schemes raises questions. One may wonder if that is linked to the actual physical behavior of tracers (i.e. the coupling between the radial flow configuration and/or transport and reaction), or to numerical artifacts linked to the numerical scheme for solving tracer concentration. For this reason, we were careful before drawing conclusions. Thus, we performed a mesh convergence study, which is reported in the next section, with a one-dimensional radial geometry.

4.3 Mesh convergence study and sensitivity analysis using a one-dimensional radial geometry

In order to quantify more precisely the benefit of a second-order scheme for the chemical tracers equation (17) resolution, the reservoir model has been simplified to a 1D radial model. We did not consider a 2D model because previous 3D results showed similar and independent responses from each layer, as a consequence of the vertical-to-horizontal permeability anisotopy and of the absence of any capillary interaction between layers. In this framework, it has been possible to simulate SWCTT on very refined grids and study the convergence of solutions with mesh size. The simulations have been performed with the O(1), O(2^m) and O(2^s) schemes. For each case, the residual oil saturation has been estimated using equation (33) [1, 2].

4.3.1 Mesh convergence study

In order to carry out a mesh convergence study that aims at evaluating the benefits of the O(2^s) and O(2^m) schemes over the O(1) scheme, the near wellbore reservoir model presented in Section 4.1 has been turned into a one-dimensional radial model by merging all the layers into a single one of 15 m height, 100 mD average absolute permeability and 10% porosity. Because of its static and dynamic properties which are homogeneous and isotropic, this domain is only meshed in the radial coordinate r. A grid refinement has been performed by reducing the radial gridblock size to $\Delta r_n=\Delta r_0/n$ with a refinement level n equal to 1, 2, 4, 8 and 16, where Δr ₀ = 0.56 m is the coarsest gridblock size. Hence, the radial gridblock size ranges from 3.5 cm to 0.56 m. In each case, the time step Δt is computed with the CFL relationship (23) as explained in Section 3.3.

As an illustration, Figure 11 reports the spatial distribution of ester, alcohol and tracer concentrations, obtained with the O(1) scheme with a mesh refinement level n = 16, after the ester placement at t = 12 days, after the shut-in period at t = 15 days, and during the recovery period at t = 15.5 days. Figure 5 compares the simulated tracer, ester and alcohol outlet aqueous concentration profiles using the O(1), O(2^m) and O(2^s) schemes when performing a grid refinement. Symbols refer to the maximum concentration profiles peak amplitudes from which the alcohol and ester arrival times t _a and t _e are estimated, and the residual oil saturation S _orw is computed according to equation (33). We will come back to that S _orw interpretation later on. Regarding convergence of predicted concentration profiles with grid refinement, Figure 5 shows a general trend which is expected: the smaller the gridblock size, that is the higher n, the less diffuse the concentration profiles and the higher their maximum peak amplitudes are. It can also be observed that the tracer, ester and alcohol arrival times are very sensitive to the scheme that is used and to the grid refinement level; this is discussed later on. Figures 6– 8 report the normalized aqueous tracer, ester and alcohol outlet concentration profiles simulated using the O(1), O(2^m) and O(2^s) schemes when performing a grid refinement. These figures show that (i) the smaller the gridblock size is, the less diffuse the concentration profiles are, (ii) a convergence trend can be observed in the peak width at mid height as n increases, and (iii) the O(2^s) scheme is less diffuse than the O(2^m) scheme, which is much less diffuse than the O(1) scheme, when n is fixed. Evolution of the tracer concentration profile width at mid height as a function of n using O(1), O(2^m) and O(2^s) schemes is reported in Figure 9b. The tracer concentration profiles width at mid height approximately converges as a power law of n, whose exponent is close to −0.4 for O(1) and O(2^m) schemes and −0.6 for the O(2^s) scheme, before converging towards the 0.5 day tracer injection duration (see Tab. 3).

Figure 9a reports the evolution of the tracer concentration profiles maximum peak amplitudes as a function of the grid refinement level n when using O(1), O(2^m) and O(2^s) schemes. The tracer concentration profiles maximum peak amplitude approximately converges as a power law of n whose exponent is −0.41 for the O(1) scheme, −2.15 for the O(2^s) scheme, approximately, and varies between −0.94 and −3.6 for the O(2^m) scheme. Figure 9 clearly demonstrates the superiority of the second-order schemes convergence over the first-order one; the O(2^s) scheme convergence rate appears to be higher than the O(2^m) one. Table 5 summarizes and compares several convergence indicators such as the tracer, ester and alcohol maximum peak concentrations max C _i,W (i = t, e, a) and the tracer width at mid height $\mathrm{s}\mathrm{upp} (C_{t,W}\ge \frac{\max C_{t,W}}{2})$, using first- and second-order schemes when performing a grid refinement.

Figure 10 compares in a detailed manner the simulated tracer, ester and alcohol outlet concentration profiles obtained with the first-order scheme and the highest spatial resolution n = 16 with the ones obtained with the second-order schemes O(2^m) and O(2^s), with the spatial resolutions n = 4 and n = 2 that approximately match the reference O(1) profiles, respectively. This figure shows that the $O\left(2^{\mathrm{m}}\right){|}_{n=4}$ and $O\left(2^{\mathrm{s}}\right){|}_{n=2}$ profiles approximately frame the $O\left(1\right){|}_{n=16}$ profiles: therefore, the O(2^m) and O(2^s) schemes reproduce the O(1) profiles with approximately 4–8 times less cells, respectively.

To finish with, the alcohol outlet concentration profile tail displayed in Figures 8 and 10 deserves some comments. As a matter of fact, the alcohol concentration profile exhibits a plateau in its tail, which is sharply captured for n = 16 only when using the first-order scheme, for n ≥ 4 when using the O(2^m) scheme, and for n ≥ 2 when using the O(2^s) scheme. This plateau is due to the ester degradation into alcohol, that happens most significantly during the shut-in period but also to a lesser extent during ester placement, where a small amount of alcohol is generated from the ester and pushed downstream the ester, as clearly demonstrated in Figure 11.

4.3.2 Residual oil saturation interpretation from alcohol and ester concentration profiles

We now come back to the residual oil saturation to waterflooding S _orw that can be determined from the alcohol and ester outlet concentration profiles reported in Figure 5, according to equation (33). As explained above, ester is late compared to the alcohol whenever S _orw  ≠ 0 because it undergoes partition between water and oil, oil being immobile, whereas the alcohol is always transported in mobile water. Alcohol and ester arrival times t _a and t _e are estimated from their maximum concentration profiles peak amplitudes, shown with symbols in Figure 5. Because of the plateau that exhibits the alcohol production tail and which is due to the forward flow, it is worth noting that it would be meaningless to derive the alcohol arrival time from the averaged alcohol concentration profile $\frac{1}{\Delta t}{\int }_0^{\Delta t}\mathrm{}{C}_{a,W}(t-t_0) \mathrm{d}t$, which would integrate that plateau.

Figure 12 and Table 6 compare the residual oil saturation $S_{\mathrm{orw}}$ computed from the ester and alcohol arrival times when using first- and second-order schemes and performing a grid refinement, as well as arrival times t _t  − t ₀, t _e  − t ₀ and t _a  − t ₀. Even if the second-order schemes yield again more accurate results compared to the first-order one – which we will focus on hereafter –, it is worth noting that the exact imposed residual oil saturation of S _orw  = 0.2 is at best predicted with a relative error of about 5% for the $O\left(2^{\mathrm{m}}\right){|}_{n=4}$ case, as shown in Figure 12d, whereas that error is about 35% in the $O\left(2^{\mathrm{m}}\right){|}_{n=1}$ case. Overall, the residual oil saturation that can be computed from the alcohol and ester arrival times shows a 10–20% relative error for n = 2 – 16, whatever the scheme that is used. To finish with, the O(2^s) scheme obviously leads to erratic S _orw estimates based on arrival times when refining the mesh because of the very sharp outlet concentration profiles it yields.

For a given numerical scheme, the evolution of S _orw with grid refinement does not show any clear trend, and seems to show unexplainable fluctuations. For the second-order schemes, fluctuations of S _orw with n stem from a large peak width, which makes the determination of maximum peak times inaccurate. So, although second-order schemes yield accurate concentration profiles, S _orw determination remains inaccurate but this inaccuracy is not due to numerical reasons but has physical origin, that is the ester slug size. The latter should have been smaller to obtain a narrow sharp peak. Regarding the first-order scheme, the determination of the maximum peak time is made easier but that scheme induces some numerical dispersion. Indeed, numerical dispersion leads to a smoothing of physical (real) concentration variations in space, which better defines the maximum concentration and corresponding arrival time for a slug of finite size recovered over a non-negligible time interval. Therefore, the benefit of a second-order numerical scheme is maximized for fairly small ester slug sizes, because minimal numerical dispersion preserves the peak sharpness and minimizes peak amplitude attenuation. For slugs of finite non-negligible size, a less accurate first-order scheme may be sufficient because numerical dispersion in space leads to a bell-shaped peak with the possibility to identify a maximum concentration and a unique arrival time. All above interpretation stands for ideal flow conditions obtained in homogeneous layers where physical dispersion can be neglected as assumed herein.

We now highlight the sensitivity of such S _orw estimates based on the simulated alcohol and ester arrival times, as a function of grid refinement and numerical scheme. More precisely, we consider that two systematic errors can be done in reading the outlet concentration profiles maximum peak amplitudes: instead of reading the actual maximum values, maxC _a,W and maxC _e,W , one would instead read maxC _a,W  − ΔC and maxC _e,W  − ΔC, with ΔC set to 50 or 25 ppm. Doing so, ester and alcohol arrival times t _i shift to $t_i^{+}=t_i+\Delta t_i^{+}$ and $t_i^{-}=t_i-\Delta t_i^{-}$ with i = a, e. Table 7 reports the corresponding estimated residual oil saturations when considering a systematic early shift in the alcohol and ester arrival times $t_a^{-}$ and $t_e^{-}$, a systematic late shift in the alcohol and ester arrival times $t_a^{+}$ and $t_e^{+}$, an early shift in the alcohol arrival time $t_a^{-}$ and a late shift $t_e^{+}$ in the ester one, or a late shift in the alcohol arrival time $t_a^{+}$ and an early shift $t_e^{-}$ in the ester one. The corresponding estimated residual saturations, which are denoted $S_{\mathrm{orw}}^{[\pm \pm ]}$ in Table 7 for each grid refinement and scheme order, lead to a maximum residual oil saturation estimate uncertainty range $[S_{\mathrm{orw}}^{-},S_{\mathrm{orw}}^{+}]$ with $S_{\mathrm{orw}}^{-}=\min S_{\mathrm{orw}}^{[\pm \pm ]}$ and $S_{\mathrm{orw}}^{+}=\max S_{\mathrm{orw}}^{[\pm \pm ]}$, and a residual oil saturation estimate uncertainty $\Delta S_{\mathrm{orw}}=S_{\mathrm{orw}}^{+}-S_{\mathrm{orw}}^{-}$. This uncertainty analysis has not been carried out with the O(2^s) scheme that obviously leads to erratic S _orw estimates based on arrival times when refining the mesh because of the very sharp outlet concentration profiles it yields, as explained above.

Figure 13 reports the resulting residual oil saturation estimates lower and upper uncertainty bound $S_{\mathrm{orw}}^{\pm }$ as a function of the grid refinement level and scheme order. It is worth noting that the shaded ΔS _orw estimate uncertainty area frames the S _orw estimated from the alcohol and ester outlet concentration profiles maximum peak amplitudes as well as the exact S _orw of 0.2 for all grid refinements and scheme orders. That uncertainty range decreases quickly as n increases, as shown in Figures 13e and 13f, and evolves as a power law of the grid refinement level n with an exponent close to −0.9 for the first and second order schemes, as shown in Figures 13g and 13h. Figures 13a– d also show that the second-order scheme residual oil saturation uncertainty range is narrower than the first-order scheme one.

5 Conclusion and perspectives

In this work, we have tested an improved numerical scheme for simulating accurately the SWCTT process. This scheme was designed in order to be fast to implement with a classic first-order upwind scheme as starting point in order to be easily integrated in industrial reservoir simulators. Thus, the resolution of the pressure and phases mass conservation has been decoupled from the chemical tracers resolution. The resolution of the pressure and phases mass conservation is based on standard implicit first-order scheme which provides unconditional stability and robustness. The chemical tracer resolution is based on an upwind explicit scheme which has been upgraded to the second-order in time and to the second-order in space of MUSCL type. In order to prove our method viability for realistic application cases, we simulated a SWCTT process on a three-dimensional model composed of 18 layers. These simulations show the robustness of our upgraded numerical scheme and the improved accuracy of the chemical tracers profiles by reducing numerical dispersion, with only a 6% increase in the computation time. In order to further compare first- and second-order schemes for the chemical tracer, we performed a mesh convergence analysis on a one-dimensional radial simulation model. The obtained results have been analysed in terms of concentration profiles maximum peak amplitude and width at mid height. Peak smoothing due to the numerical dispersion is strongly reduced with a second-order scheme. As a consequence, the second-order scheme gives good approximation with coarser grids compared to the first-order scheme.

In a forthcoming work, we expect to upgrade the chemical tracer scheme by considering a sub-cycling technique for its resolution to reduce significantly the simulation times. In addition, we plan to consider more complex physics involving dispersion, adsorption, temperature-dependent partition coefficient and non-instantaneous partitioning. Due to this complex physics, a dedicated SWCTT history matching methodology, not relying only on the maximum peak amplitudes of ester and alcohol outlet concentration profiles, is recommended for the S_orw assessment. Actually, all these effects combined with the numerical resolution specificities may have a significant impact on the S_orw uncertainty range, which is the main output of a SWCTT. Interpretation of such a SWCTT may become even trickier when the near wellbore region is made of different layers that hydrodynamically communicate as shown in [28]. That complexity justifies to dispose of:

• A precise numerical scheme, as the one proposed in this work.

• A dedicated SWCTT uncertainty analysis and history matching methodology to design the test and interpret results.

References

All Tables

All Figures