Regular Article
MUSCL scheme for Single Well Chemical Tracer Test simulation, design and interpretation
IFP Energies nouvelles, 1 et 4 avenue de BoisPréau, 92852 RueilMalmaison Cedex, France
^{*} Corresponding author: benjamin.braconnier@ifpen.fr
Received:
23
February
2018
Accepted:
15
November
2018
Our paper presents an improved numerical scheme to simulate Single Well Chemical Tracer Test (SWCTT) method. SWCTT is mainly applied to determine the residual oil saturation of reservoirs. It consists in injecting an aqueous slug of a primary tracer into the reservoir formation and displacing it at a certain distance from the well. This tracer is partly miscible with oil on the one hand, and generates in situ a secondary tracer on the other hand. As a consequence, a shift is observed between the primary and the secondary tracers arrival times when production is resumed. This time shift is used to evaluate the residual oil saturation. In our paper, we propose a numerical scheme based on a fractional time stepping technique to decouple the resolution of the phases mass conservation equations and the chemical tracers mole conservation equations. For the phases resolution, we use an implicit scheme to ensure stability and robustness. For the chemical tracers, we propose an explicit secondorder scheme in time and in space via MUSCL technique to improve the tracers timeshift calculation. The proposed numerical method is implemented on a realistic simulation model consisting of a vertical well crossing a reservoir consisting of a stack of homogeneous layers. By reducing the numerical dispersion, the proposed scheme improves the accuracy of predicted concentration profiles, without significantly increasing the computation time. Finally, the advantages of using a secondorder scheme for residual oil saturation assessment are discussed on the basis of a radial 1D mesh convergence study.
© B. Braconnier et al., published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
List of symbols
Symbols concerning the rockSymbols concerning the phase ψ
kr _{ ψ } : Relative permeability
: Maximum relative permeability
v _{ ψ } : Hydrodynamic gradient
q _{ ψ } : Well source terms per volume
Symbols concerning the chemical tracer α
C _{ α,ψ } : Mole concentration in the phase ψ
C _{ α,R } : Adsorbed mole concentration
: Maximum adsorbed mole concentration
b _{ α } : Langmuir coefficient
: Molecular diffusion in phase ψ
F _{ α } : Molar advection flux
J _{ α } : Molar diffusion/dispersion fluxes
: Degradation kinetic constant
K _{ α,ψ } − ψ _{ α } : Partition coefficient for phase ψ
Symbols concerning the reaction
γ _{ α,β} : Stoechiometric coefficient
: Chemical reaction half life time
: Chemical reaction kinetic constant
Symbols concerning the mesh
Π_{prod} : Set of perforations of production wells
Π_{inj} : Set of perforations of injection wells
Symbols concerning the mesh cell K
KL : Face between cell and cell
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
: Finite volume flux of hydrodynamic gradient at face KL
: Discrete adsorption operator
: Numerical scheme operator in concise form
: Limited gradient of at face KL
: Secondorder approximation of at face KL
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
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
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 firstorder upwind scheme [5]. Regarding the flow model, we adopt a classical blackoil 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], DattaGupta et al. [9] where other chemicals such as polymer, surfactant or foam were involved. Also, many sophisticated highorder 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 firstorder upwind finite volume methods [12]. For this, we select among wellknown 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 blackoil 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 highorder scheme for the resolution of these variables and we consider a firstorder 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 secondorder scheme in time and space. For the secondorder in time, we propose to use a predictorcorrector scheme [14, 15] (or also referred as the Heun’s method) which is quite simple to implement. For the secondorder 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 secondorder 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 timestepping 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 secondorder 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 1D mesh convergence study, the benefit of a secondorder scheme is discussed regarding the assessment of residual oil saturation.
2 Physical model
We consider a model for a threephase 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 Blackoil model
To model the water and hydrocarbon components flow, we consider a standard blackoil model [6, 12] which is defined in mass per unit volume of saturated porous medium as (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 modeling production wells and the inflow source terms devoted to injection wells. In the following, we suppose that we only inject water or gas in the reservoir: .
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 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 (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 and the hydrodynamic gradient such that .
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 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 (3) where n _{ α } is the tracer mole number given by (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 (5) where is the maximum adsorbed mole concentration of chemical tracer α on the rock. It is deduced from the adsorption maximum mass fraction through the relation where ρ _{ R } is the rock density. The coefficient is deduced from Langmuir coefficient b _{ α } (dimensionless) with the relation .
The advection flux of the chemical tracer is given by (6)
The diffusion and dispersion contributions are written as follows in a homogeneous porous medium (7) where is the molecular diffusion constant of the chemical tracer α in the phase ψ, τ is the tortuosity and 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 (8) where we recall that is the injected reference phase flow rate per volume and is the injected mole concentration of chemical tracer α.
The production terms are given by (9) where we recall that is the produced flow rate of phase .
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 (10) where is the kinetic constant of tracer . In practice, the half life time is commonly used to characterize this kinetics and is related to the former constant with the relation .
A classical SWCTT generally involves a single chemical reaction. Let us then suppose that α and are the two reactants of this chemical reaction, with α the primary tracer to be followed. Reaction of α with generates chemical tracer β and other products denoted . This type of reaction reads (11) where is the stoechiometric coefficient of tracer β with respect to tracer α. We suppose also that the reactant is present in very large quantity compared to the reaction consumption quantities (usually is the reference phase molecule, namely H_{2}O for an aqueous carrier phase) such that there is no need to model the reactant consumption in our model. In addition, the products 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 and , this type of reaction simplifies as (12)
This chemical reaction has its own kinetics characterized by a half life time that we will denote or the kinetic constant . 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 reactionrateimpacting 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 and writes as follows: (13)
This chemical reaction couples the mole conservation equation of chemical tracers and via and 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 as follows (14)
Using this relation, the chemical tracer mole conservation equation can be formulated with its mole concentration in the reference phase as unknown (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 firstorder explicit discretization of chemical tracers mole conservation equations, we describe its extension to a secondorder scheme in time and space.
3.1 Domain discretization
Let be an admissible finite volume mesh of the reservoir given by a family of control volumes or cells noted K. For any K of , K is its measure, 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 , that is .
To construct our secondorder in space scheme, we will need to define the center of each cell K. For each neighbor cell , we denote d _{ K,KL } the distance of the point 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 . For instance, K is a perforated cell if it exists a perforation or such that M(j) = K. In the next subsections, we note and for concision.
Finally, we introduce an increasing sequence of discrete times such that t ^{0} = 0 and t ^{ N } = T. Then the time interval (0, T) is subdivided into N variable time steps Δt = t ^{ n+1} − t ^{ n }, 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 highorder scheme for the water and hydrocarbon components mass conservation equations. In fact, the main requirement is only to dispose of a robust and unconditionallystable resolution scheme. Therefore, we propose to use an implicit scheme in time and an upwind finite volume method for the spatial derivatives (16) where is a finite volume discretization of the flux at the time t ^{ n+1}. For each phase, we define the operator 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 secondorder numerical scheme is implemented because it is robust and more accurate than a firstorder 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 nearwellbore 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 diffusiondispersion terms (17)
In the next subsection, we firstly describe the chemical tracer basic scheme that is firstorder in space and time. This scheme is extended in Section 3.3.2 to a secondorder scheme in time by using a predictorcorrector 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 firstorder scheme
We propose an explicit scheme in time and upwind finitevolume approximations for the spatial derivatives, as follows: (18) where we denote and the total chemical tracer mole per unit volume respectively at times t ^{ n } and t ^{ n+1}. 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 }, 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 (19) where denotes the adsorption operator and . Then, the molar concentration is computed from the total mole number of chemical tracer in the cell by inverting the previous relation with l = n + 1. This inversion leads to the resolution of a secondorder polynomial with the unknown (with the previous notation) which admits only one positive root (20) with (21)
The values and denote respectively the mole concentration in the cell and if a firstorder scheme in space is considered.
In order to simplify the notation when describing how to extend our numerical scheme to secondorder in time and space, we propose to rewrite the system (18) in a more concise form to obtain from as follows (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 because they are taken at time t ^{ n+1} from implicit resolution of (16). In addition, as only a firstorder scheme is considered in this subsection, we set and .
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 (23) with (24) and (25)
Application of the above firstorder scheme in Section 4 will be performed with an approximate CFL criterion to ensure stability.
3.3.2 Secondorder scheme in time for the chemical tracers
The numerical method used in this paper to build a secondorder method is classical [14, 15]. Considering formula (22), secondorder scheme is written as (26)
This scheme involves the computation of an intermediate solution . The coefficient a _{1}, a _{2}, h _{1} and h _{2} 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.
Coefficient values to recover classical secondorder scheme in time.
Even if switching from one scheme to another is simple, only the predictorcorrector scheme will be used in the simulation Section 4 for the sake of briefness of this paper.
In the framework of a firstorder scheme in space, the slopes of the chemical tracer mole concentration are not computed such that , , and . The next subsection describes a method to further improve the resolution by including also a secondorder scheme in space, based on the MUSCL method.
3.3.3 Secondorder 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 and are reconstucted by using the states with . First, let us first determine the gradient of the chemical tracer molar concentration across each interface that separates the cell and the cell (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 such that is strictly negative and minimum. If one such cell exists, we compute the gradient (28) and then compute a limited slope (29) where the function 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 . 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 secondorder approximate molar concentration at the interface between the cells K and L as (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 and to compute .
By using these reconstructed molar concentrations, schemes (22) or (26) become secondorder in space.
As slope limiter, we can choose the Minmod limiter defined as (31) or the Superbee limiter (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 firstorder or a secondorder 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 .
Layers thickness, horizontal and vertical permeabilities.
The considered well test is supposed to be an isothermal twophase 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 and , 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, and 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 m^{3}/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 , and an ester, denoted e, are injected for day with the same flow rate and aqueous mass fractions of ppm, followed by a 1.5 days water drive and a shutin 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 m^{3}/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 shutin period at t = 15 days, and during the recovery period at t = 15.5 days.
Well test schedule and fluids composition.
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] (33) where t _{ a } and t _{ e } denote the alcohol and ester arrival times, respectively; t _{0} = 15 days denotes the backflow starting time and is the ester partition coefficient, defined as the ratio of the ester molar concentration in the oil phase (expressed in mol/m^{3}) 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:

Instantaneously partition between the mobile water phase and the residual immobile oil phase according to the constant partition coefficient K _{ e,W−O }.

Undergo a slow hydrolysis that yields alcohol according to the firstorder kinetics where is the kinetic constant of ester hydrolysis producing alcohol.
The corresponding half life time has been set to 3 days so that the alcohol is mainly generated during the shutin period, whose duration (3 days) is however not too long so that a nonnegligible 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 pistontype, and its recovery that is delayed compared to the alcohol one, one expects its displacement to be caterpillarlike 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 firstorder scheme and with a secondorder scheme both in time and space. For the sake of simplicity, we will denote hereafter the firstorder scheme O(1), the secondorder scheme with a minmod slope limiter O(2^{m}), and the secondorder scheme with a superbee slope limiter O(2^{s}).
4.2 Simulations on a threedimensional regular grid
We consider a threedimensional 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 xy 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 secondorder resolution of the chemical tracer equation (17) and its relevance to achieve threedimensional 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 shutin period at t = 15 days,

and finally during the recovery period at t = 15.5 days.
Fig. 1 Spatial distribution of ester, alcohol and tracer concentrations, after the ester placement at t = 12 days, obtained with the first and secondorder schemes. 
Fig. 2 Spatial distribution of ester, alcohol and tracer concentrations, after the shutin period at t = 15 days, obtained with the first and secondorder schemes. 
Fig. 3 Spatial distribution of ester, alcohol and tracer concentrations, during the recovery period at t = 15.5 days, obtained with the first and secondorder schemes. 
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 shutin 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 shutin 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 shutin 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).
Fig. 4 Comparison of the simulated tracer, ester and alcohol outlet concentration profiles using first and secondorder schemes with the coarsest gridblock size Δr _{0} = 0.56 m (see the text). Symbols denote the maximum concentration profiles peak amplitudes. 
Fig. 5 Comparison of the simulated tracer, ester and alcohol concentration profiles using first and secondorder schemes when gradually refining the gridblock size to , with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size (see the text). Symbols denote the maximum concentration profiles peak amplitudes. Normalized profiles are reported in Figures 6– 8. 
Fig. 6 Normalized aqueous tracer outlet concentration profiles simulated using first and secondorder schemes when refining the gridblock size to with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size. Evolutions of the maximum peak amplitude and of the non normalized concentration profiles mid height width as a function of are reported in Figure 9. 
Fig. 7 Normalized aqueous ester outlet concentration profiles simulated using first and secondorder schemes when refining the gridblock size to with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size. 
Fig. 8 Normalized aqueous alcohol outlet concentration profiles simulated using first and secondorder schemes when refining the gridblock size to with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size. 
Fig. 9 Evolution (a) of the tracer concentration profiles maximum peak amplitudes and (b) of the concentration profiles width at mid height as a function of using and minmod and superbee schemes (see the text; the corresponding numerical values are given in Tab. 5). 
Fig. 10 Comparison of the simulated (a) tracer, (b) ester and (c) alcohol outlet concentration profiles obtained with the firstorder scheme and the highest spatial resolution n = 16 with the ones obtained with the secondorder schemes O(2^{m}) and O(2^{s}), with the spatial resolutions n = 4 and n = 2 that approximately match the profiles, respectively. This figure shows that the and profiles approximately frame the profile: therefore, the O(2^{m}) and O(2^{s}) schemes reproduce the O(1) profiles with approximately 4 to 8 times less cells, respectively. 
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 % 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 secondorder 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 secondorder scheme for chemical tracer modeling and its ability to simulate a realistic threedimensional 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 onedimensional radial geometry.
Comparison of the tracer, ester and alcohol arrival times and of the resulting estimate for each scheme (see the text).
4.3 Mesh convergence study and sensitivity analysis using a onedimensional radial geometry
In order to quantify more precisely the benefit of a secondorder 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 verticaltohorizontal 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 onedimensional 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 with a refinement level n equal to 1, 2, 4, 8 and 16, where Δr _{0} = 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 shutin 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).
Fig. 11 Spatial distribution of ester, alcohol and tracer concentrations (a) after the ester placement at t = 12 days, (b) after the shutin period at t = 15 days, (c) during the recovery period at t = 15.5 days. These profiles are obtained with the firstorder scheme with a grid refinement level n = 16 (see the text). 
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 secondorder schemes convergence over the firstorder 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 , using first and secondorder schemes when performing a grid refinement.
Comparison of several convergence indicators such as the tracer, ester and alcohol outlet concentration profiles maximum peak concentrations max C_{(i,W)} (i = t, e, a) and the tracer width at mid height , when using first and secondorder schemes and performing a grid refinement.
Figure 10 compares in a detailed manner the simulated tracer, ester and alcohol outlet concentration profiles obtained with the firstorder scheme and the highest spatial resolution n = 16 with the ones obtained with the secondorder 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 and profiles approximately frame the 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 firstorder 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 shutin 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 , which would integrate that plateau.
Figure 12 and Table 6 compare the residual oil saturation computed from the ester and alcohol arrival times when using first and secondorder schemes and performing a grid refinement, as well as arrival times t _{ t } − t _{0}, t _{ e } − t _{0} and t _{ a } − t _{0}. Even if the secondorder schemes yield again more accurate results compared to the firstorder 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 case, as shown in Figure 12d, whereas that error is about 35% in the 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.
Fig. 12 Comparison of the residual oil saturation S _{ orw } computed from the ester and alcohol arrival times, according to equation (33), when using first and secondorder schemes and performing a grid refinement (see also Tab. 6). Arrival times t _{ t } − t_{0}, t _{ e } − t_{0} and t _{ a } − t _{0} are computed from the tracer, ester and alcohol outlet concentration profiles maximum peak amplitudes. 
Comparison of the residual oil saturation computed from the ester and alcohol arrival times, according to equation (33), when using first and secondorder schemes and performing a grid refinement (see also Fig. 12). Expected residual oil saturation is . Arrival times , and are computed from the tracer, ester and alcohol outlet concentration profiles maximum peak amplitudes.
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 secondorder schemes, fluctuations of S _{ orw } with n stem from a large peak width, which makes the determination of maximum peak times inaccurate. So, although secondorder 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 firstorder 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 nonnegligible time interval. Therefore, the benefit of a secondorder 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 nonnegligible size, a less accurate firstorder scheme may be sufficient because numerical dispersion in space leads to a bellshaped 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 and 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 and , a systematic late shift in the alcohol and ester arrival times and , an early shift in the alcohol arrival time and a late shift in the ester one, or a late shift in the alcohol arrival time and an early shift in the ester one. The corresponding estimated residual saturations, which are denoted in Table 7 for each grid refinement and scheme order, lead to a maximum residual oil saturation estimate uncertainty range with and , and a residual oil saturation estimate uncertainty . 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.
Comparison of the residual oil saturation estimates uncertainty ranges resulting from a 25 or 50 ppm error in the alcohol and ester outlet concentration profiles maximum peak amplitudes measurement, when using first and minmod secondorder schemes and performing a grid refinement (see the text). Negative values are reported as zero values.
Figure 13 reports the resulting residual oil saturation estimates lower and upper uncertainty bound 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 secondorder scheme residual oil saturation uncertainty range is narrower than the firstorder scheme one.
Fig. 13 Comparison of the residual oil saturation estimates uncertainty ranges resulting from a 50 (left column) or 25 ppm (right column) error in the alcohol and ester outlet concentration profiles maximum peak amplitudes measurement, when using first and minmod secondorder schemes and performing a grid refinement (see the text). 
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 firstorder 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 firstorder scheme which provides unconditional stability and robustness. The chemical tracer resolution is based on an upwind explicit scheme which has been upgraded to the secondorder in time and to the secondorder in space of MUSCL type. In order to prove our method viability for realistic application cases, we simulated a SWCTT process on a threedimensional 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 secondorder schemes for the chemical tracer, we performed a mesh convergence analysis on a onedimensional 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 secondorder scheme. As a consequence, the secondorder scheme gives good approximation with coarser grids compared to the firstorder scheme.
In a forthcoming work, we expect to upgrade the chemical tracer scheme by considering a subcycling technique for its resolution to reduce significantly the simulation times. In addition, we plan to consider more complex physics involving dispersion, adsorption, temperaturedependent partition coefficient and noninstantaneous 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
 Deans H. (1971) Method of determining fluid saturations in reservoirs. U.S. Patent 3623842. [Google Scholar]
 Deans H., Carlisle C. (2007) The SingleWell Chemical Tracer Test – A method for measuring reservoir fluid saturations in situ, Chapter 5, in: Holstein E.D. (ed), Petroleum Engineering Handbook, Vol. 5, SPE, Richardson, Texas, pp. 615–648. [Google Scholar]
 Khaledialidusti R., Kleppe J. (2017) A comprehensive framework for the theoretical assessment of the singlewellchemicaltracer tests, J. Pet. Sci. Eng. 159, 164–181, . doi: 10.1016/j.petrol.2017.09.027 . [Google Scholar]
 Braconnier B., Flauraud É., Nguyen Q.L. (2014) Efficient scheme for chemical flooding simulation, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 69, 4, 585–601. doi: 10.2516/ogst/2013189 . [CrossRef] [Google Scholar]
 Braconnier B., Preux C., Flauraud É., Tran Q.H., Berthon C. (2017) An analysis of physical models and numerical schemes for polymer flooding simulations, Comput. Geosci. 21, 5, 1267–1279. doi: 10.1007/s1059601796370 . [Google Scholar]
 Peaceman D.W. (1977) Fundamentals of Numerical Reservoir Simulation, volume 6 of Developments in Petroleum Science, Elsevier Science, Amsterdam, The Netherlands. [Google Scholar]
 Delshad M., Pope G.A., Sepehrnoori K. (1996) A compositional simulator for modeling surfactant enhanced aquifer remediation: 1. Formulation, J. Contam. Hydrol. 23, 4, 303–327. doi: 10.1016/01697722(95)001069. [Google Scholar]
 Pope G.A., Nelson R.C. (1978) A chemical flooding compositional simulator, SPE J. 18, 5, 339–354. doi: 10.2118/6725PA. [Google Scholar]
 DattaGupta A., Pope G.A., Seperhrnoori K., Thrasher R.L. (1986) A symmetric, positive definite formulation of a threedimensional micellar/polymer simulator, SPE Reserv. Eng. 1, 6, 622–632. doi: 10.2118/13504PA . [CrossRef] [Google Scholar]
 Saad N., Pope G.A., Sepehrnoori K. (1990) Application of higherorder methods in compositional simulation, SPE Reserv. Eng. 5, 4, 623–630. doi: 10.2118/18585PA . [CrossRef] [Google Scholar]
 Liu J., Delshad M., Pope G.A., Sepehrnoori K. (1994) Application of higherorder fluxlimited methods in compositional simulation, Transp. Porous Media 16, 41, 1–29. doi: 10.1007/BF01059774 . [Google Scholar]
 Trangenstein J.A., Bell J.B. (1989) Mathematical structure of the blackoil model for petroleum reservoir simulation, SIAM J. Appl. Math. 49, 3, 749–783. doi: 10.1137/0149044 . [Google Scholar]
 Khaledialidusti R., Kleppe J. (2018) Significance of geochemistry in singlewell chemicaltracer tests by coupling a multiphaseflow simulator to the geochemical package, SPE J. 23, 1126–1144. doi: 10.2118/189971PA. https://www.onepetro.org/journalpaper/SPE189971PA . [CrossRef] [Google Scholar]
 Butcher J.C. (2003) Numerical methods for ordinary differential equations, John Wiley & Sons, New York. ISBN 9780471967583. [CrossRef] [Google Scholar]
 Press W.H., + Teukolsky S.A., Vetterling W.T., Flannery B.P. (2007) Numerical recipes: the art of scientific computing, Cambridge University Press, New York. ISBN 9780521880688. [Google Scholar]
 van Leer B. (1979) Towards the ultimate conservative difference scheme. V. A secondorder sequel to Godunov’s method, J. Comput. Phys. 32, 1, 101–136. doi: 10.1016/00219991C79)901451. ISSN 00219991. http://www.sciencedirect.com/science/article/pii/0021999179901451 . [Google Scholar]
 Toro E.F. (2009) Riemann solvers and numerical methods for fluid dynamics, Springer, Heidelberg, Berlin. doi: 10.1007/b79761 . [CrossRef] [Google Scholar]
 Perthame B., Qiu Y. (1994) A variant of van leer’s method for multidimensional systems of conservation laws, J. Comput. Phys. 112, 2, 370–381. doi: 10.1006/jcph.1994.1107. ISSN 00219991. http://www.sciencedirect.com/science/article/pii/S0021999184711077 . [Google Scholar]
 Berthon C. (2006) Why the musclhancock scheme is L1stable, Numer. Math. 104, 27–46. [Google Scholar]
 Buffard T., Clain S. (2010) Monoslope and multislope muscl methods for unstructured meshes, J. Comptut. Phys. 229, 3745–3776. doi: 10.1016/j.jcp.2010.01.026. http://hdl.handle.net/1822/11976 . [CrossRef] [Google Scholar]
 Jiang J., Younis R. (2016) An efficient fullyimplicit highresolution MFDMUSCL method for twophase flow with gravity in discrete fractured media, ECMOR XV – 15th European Conference on the Mathematics of Oil Recovery, 29 August to 1 September 2016, Amsterdam, The Netherlands. doi: 10.3997/22144609.201601753 . [Google Scholar]
 Jiang J., Younis R. (2016) C1PPU schemes for efficient simulation of coupled flow and transport with gravity, ECMOR XV – 15th European Conference on the Mathematics of Oil Recovery, 08, 29 August to 1 September 2016, Amsterdam, The Netherlands. [Google Scholar]
 Jiang J., Younis R. (2017) C1PPU schemes for efficient simulation of coupled flow and transport with gravity, SPE Reservoir Simulation Conference, 20–22 February, Montgomery, Texas, USA. doi: 10.2118/182695MS. SPE182695MS. [Google Scholar]
 Berthon C. (2006) Robustness of MUSCL schemes for 2D unstructured meshes, J. Comput. Phys. 218, 495–509. [Google Scholar]
 Berthon C. (2005) Stability of the MUSCL schemes for the Euler equations, Commun. Math. Sci. 3, 133–158. [Google Scholar]
 Cheng H., Shook M.S., Malik T., Varadarajan D., Smith B.R. (2012) Interwell tracer tests to optimize operating conditions for a surfactant field trial: design, evaluation, and implications, SPE Reserv. Eval. Eng. 899, 229–242. doi: 10.2118/144899PA . [CrossRef] [Google Scholar]
 Huseby O.K., Sagen J., Dugstad O. (2012) Single well chemical tracer tests  fast and accurate simulations, SPE EOR Conference at Oil and Gas West Asia, April 2012, Muscat, Oman. doi: 10.2118/155608MS. Paper SPE 155 [Google Scholar]
 De Zwart A.H., Van Batenburg D.W., Stoll M., AlHarthi S. (2011) Numerical interpretation of single chemical tracer tests for ASP injection, 16th European Symposium on Improved Oil Recovery, April 2011, Cambridge, UK. Paper B18. [Google Scholar]
All Tables
Comparison of the tracer, ester and alcohol arrival times and of the resulting estimate for each scheme (see the text).
Comparison of several convergence indicators such as the tracer, ester and alcohol outlet concentration profiles maximum peak concentrations max C_{(i,W)} (i = t, e, a) and the tracer width at mid height , when using first and secondorder schemes and performing a grid refinement.
Comparison of the residual oil saturation computed from the ester and alcohol arrival times, according to equation (33), when using first and secondorder schemes and performing a grid refinement (see also Fig. 12). Expected residual oil saturation is . Arrival times , and are computed from the tracer, ester and alcohol outlet concentration profiles maximum peak amplitudes.
Comparison of the residual oil saturation estimates uncertainty ranges resulting from a 25 or 50 ppm error in the alcohol and ester outlet concentration profiles maximum peak amplitudes measurement, when using first and minmod secondorder schemes and performing a grid refinement (see the text). Negative values are reported as zero values.
All Figures
Fig. 1 Spatial distribution of ester, alcohol and tracer concentrations, after the ester placement at t = 12 days, obtained with the first and secondorder schemes. 

In the text 
Fig. 2 Spatial distribution of ester, alcohol and tracer concentrations, after the shutin period at t = 15 days, obtained with the first and secondorder schemes. 

In the text 
Fig. 3 Spatial distribution of ester, alcohol and tracer concentrations, during the recovery period at t = 15.5 days, obtained with the first and secondorder schemes. 

In the text 
Fig. 4 Comparison of the simulated tracer, ester and alcohol outlet concentration profiles using first and secondorder schemes with the coarsest gridblock size Δr _{0} = 0.56 m (see the text). Symbols denote the maximum concentration profiles peak amplitudes. 

In the text 
Fig. 5 Comparison of the simulated tracer, ester and alcohol concentration profiles using first and secondorder schemes when gradually refining the gridblock size to , with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size (see the text). Symbols denote the maximum concentration profiles peak amplitudes. Normalized profiles are reported in Figures 6– 8. 

In the text 
Fig. 6 Normalized aqueous tracer outlet concentration profiles simulated using first and secondorder schemes when refining the gridblock size to with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size. Evolutions of the maximum peak amplitude and of the non normalized concentration profiles mid height width as a function of are reported in Figure 9. 

In the text 
Fig. 7 Normalized aqueous ester outlet concentration profiles simulated using first and secondorder schemes when refining the gridblock size to with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size. 

In the text 
Fig. 8 Normalized aqueous alcohol outlet concentration profiles simulated using first and secondorder schemes when refining the gridblock size to with n = 1, 2, 4, 8, 16 and Δr _{0} = 0.56 m the coarsest gridblock size. 

In the text 
Fig. 9 Evolution (a) of the tracer concentration profiles maximum peak amplitudes and (b) of the concentration profiles width at mid height as a function of using and minmod and superbee schemes (see the text; the corresponding numerical values are given in Tab. 5). 

In the text 
Fig. 10 Comparison of the simulated (a) tracer, (b) ester and (c) alcohol outlet concentration profiles obtained with the firstorder scheme and the highest spatial resolution n = 16 with the ones obtained with the secondorder schemes O(2^{m}) and O(2^{s}), with the spatial resolutions n = 4 and n = 2 that approximately match the profiles, respectively. This figure shows that the and profiles approximately frame the profile: therefore, the O(2^{m}) and O(2^{s}) schemes reproduce the O(1) profiles with approximately 4 to 8 times less cells, respectively. 

In the text 
Fig. 11 Spatial distribution of ester, alcohol and tracer concentrations (a) after the ester placement at t = 12 days, (b) after the shutin period at t = 15 days, (c) during the recovery period at t = 15.5 days. These profiles are obtained with the firstorder scheme with a grid refinement level n = 16 (see the text). 

In the text 
Fig. 12 Comparison of the residual oil saturation S _{ orw } computed from the ester and alcohol arrival times, according to equation (33), when using first and secondorder schemes and performing a grid refinement (see also Tab. 6). Arrival times t _{ t } − t_{0}, t _{ e } − t_{0} and t _{ a } − t _{0} are computed from the tracer, ester and alcohol outlet concentration profiles maximum peak amplitudes. 

In the text 
Fig. 13 Comparison of the residual oil saturation estimates uncertainty ranges resulting from a 50 (left column) or 25 ppm (right column) error in the alcohol and ester outlet concentration profiles maximum peak amplitudes measurement, when using first and minmod secondorder schemes and performing a grid refinement (see the text). 

In the text 