MUSCL scheme for Single Well Chemical Tracer Test simulation, design and interpretation

. 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 second-order scheme in time and in space via MUSCL technique to improve the tracers time-shift 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 proﬁles, without signiﬁcantly increasing the computation time. Finally, the advantages of using a second-order scheme for residual oil saturation assessment are discussed on the basis of a radial 1D mesh convergence study.

MUSCL scheme for Single Well Chemical Tracer Test simulation, design and interpretation Abstract. 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 second-order scheme in time and in space via MUSCL technique to improve the tracers time-shift 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 second-order scheme for residual oil saturation assessment are discussed on the basis of a radial 1D mesh convergence study. 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

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 secondorder 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][17][18][19][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][22][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.

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.

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 where U is the rock porosity. For each phase denoted w = W, O, G, S w is the saturation, q w is the mass density and q w denotes the source term per volume. It will be useful to distinguish the outflow source terms q þ w ! 0 modeling production wells and the inflow source terms q À w 0 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 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 where K is the rock permeability, kr w is the relative permeability for the phase w, l w is the pure phase viscosity, P w is the pressure of the phase w and g is the gravity acceleration. In order to simplify notations, we introduce the phase mobility k w ¼ kr w =l w and the hydrodynamic gradient v w ¼ ÀK rP w À q w g À Á ; such that u w ¼ k w v w .

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 a 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 w a 2 fW ; O; Gg for tracer a. This phase has to be present wherever the chemical tracer is present.
The mass conservation equation of the chemical tracer a is converted to a mole conservation equation dividing it by the molar mass M a . The mole conservation equation is given by where n a is the tracer mole number given by where C a,w is the mole concentration of tracer a in the phase w and C a,R the mole concentration of tracer adsorbed on the rock. C a,R is governed by a Langmuir isotherm which has been extended to the case of a phase partitioning tracer where C max a;R is the maximum adsorbed mole concentration of chemical tracer a on the rock. It is deduced from the adsorption maximum mass fraction X max a;R through the relation C max a;R ¼ q R X max a;R =M a where q R is the rock density. The coefficientb a is deduced from Langmuir coefficient ba The advection flux of the chemical tracer is given by The diffusion and dispersion contributions are written as follows in a homogeneous porous medium where D df a;w is the molecular diffusion constant of the chemical tracer a in the phase w, s is the tortuosity and D ds w 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 where we recall that q À w a is the injected reference phase flow rate per volume and C inj a;w a is the injected mole concentration of chemical tracer a.
The production terms are given by where we recall that q þ w is the produced flow rate of phase w.
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 where k Deg a is the kinetic constant of tracer a. In practice, the half life time t Deg a;1=2 is commonly used to characterize this kinetics and is related to the former constant with the relation t Deg a;1=2 ¼ ln 2=k Deg a . A classical SWCTT generally involves a single chemical reaction. Let us then suppose that a and R are the two reactants of this chemical reaction, with a the primary tracer to be followed. Reaction of a with R generates chemical tracer b and other products denoted P. This type of reaction reads where c a;b is the stoechiometric coefficient of tracer b with respect to tracer a. We suppose also that the reactant R is present in very large quantity compared to the reaction consumption quantities (usually R 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 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 R and P, this type of reaction simplifies as This chemical reaction has its own kinetics characterized by a half life time that we will denote t Rea ab;1=2 or the kinetic constant k Rea ab ¼ ln 2=t Rea ab;1=2 . 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-rateimpacting 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 a and b writes as follows: 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 a, 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 a;wÀw a as follows C a;w ¼ K a;wÀw a C a;w a 8w 2 fW ; O; Gg: Using this relation, the chemical tracer mole conservation equation can be formulated with its mole concentration in the reference phase as unknown 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.

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.

Domain discretization
Let 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 M h , |K| is its measure, KL ¼ oK \ oL 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 N ðKÞ, that is N ðKÞ ¼ fL 2 M h ; oK \ oL 6 ¼ ;g. To construct our second-order in space scheme, we will need to define the center M K of each cell K. For each neighbor cell L 2 N ðKÞ, we denote d K,KL the distance of the point M K to the interface KL and d L,KL the distance of the point M L to the interface KL.
Let P prod be the set of perforations of the production wells and P inj be the set of perforations of the injection wells and M the set of all the perforated cells in M h . For instance, K is a perforated cell if it exists a perforation j 2 P prod or j 2 P inj such that M(j) = K. In the next subsections, we note P p ðKÞ ¼ fj 2 P prod jMðjÞ ¼ Kg and P i ðKÞ ¼ fj 2 P inj jMðjÞ ¼ Kg for concision.
Finally, we introduce an increasing sequence of discrete times ft n g 0 n N such that t 0 = 0 and t N = T. Then the time interval (0, T) is subdivided into N variable time steps Dt = tn+1 À tn, 0 n N À 1.

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 where v nþ1 w;KL is a finite volume discretization of the flux R KL v w Á n KL dr 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].

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 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][17][18][19][20] is also considered because monotonous and stable [24,25].

Chemical tracer first-order scheme
We propose an explicit scheme in time and upwind finitevolume approximations for the spatial derivatives, as follows: where we denote n n K and n nþ1 K 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 , n l K 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 where A denotes the adsorption operator and ðb a Þ Then, the molar concentration ðC a;w a Þ nþ1 K is computed from the total mole number of chemical tracer n nþ1 K 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 a;w a Þ l K (with the previous notation) which admits only one positive root with The values ðC a;w a Þ n KL and ðC a;w a Þ n LK denote respectively the mole concentration in the cell ðC a;w a Þ n K and ðC a;w a Þ n L 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 nþ1 K from n n K as follows 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 F because they are taken at time t n+1 from implicit resolution of (16). In addition, as only a first-order scheme is considered in this subsection, we set ðC a;w a Þ Application of the above first-order scheme in Section 4 will be performed with an approximate CFL criterion to ensure stability.

Second-order 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), second-order scheme is written as This scheme involves the computation of an intermediate solution n ½nþ1 K . 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.
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 ðC a;w a Þ . The next subsection describes a method to further improve the resolution by including also a secondorder scheme in space, based on the MUSCL method.

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 ðC a;w a Þ l KL and ðC a;w a Þ l KL are reconstucted by using the states ðC a;w a Þ l K 0 with K 0 2 M h . First, let us first determine the gradient of the chemical tracer molar concentration across each interface KL that separates the cell K and the cell L 2 N ðKÞ 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 dw 2 N ðKÞ such that n KL Á n KL dw is strictly negative and minimum. If one such cell exists, we compute the gradient and then compute a limited slopẽ r KL C n a;w a Á n KL ¼ Lðr KL C n a;w a Á n KL ; r KL dw C n a;w a Á n KL dw Þ; ð29Þ where the function 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 setr KL C n a;w a ¼ 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 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 slopẽ r LK C n a;w a Á n KL and to compute ðC a;w a Þ n LK . 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 Lðx; yÞ ¼ 1 2 x jxj þ y jyj minðjxj; jyjÞ; ð31Þ or the Superbee limiter Lðx; yÞ ¼ 1 2 x jxj þ y jyj max minð2jxj; jyjÞ; minðjxj; 2jyjÞ ð Þ : 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.

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 kr W ðSÞ ¼ kr 0 W S n W and kr O ðSÞ ¼ kr 0 O ð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, kr 0 W ¼ 0:2 and kr 0 O ¼ 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 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 t, and an ester, denoted e, are injected for 0:5 day with the same flow rate and aqueous mass fractions of C inj t;W ¼ C inj e;W ¼ 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  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 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 orw ¼ 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 0 = 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 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: (i) Instantaneously partition between the mobile water phase and the residual immobile oil phase according to the constant partition coefficient K e;W ÀO . (ii) Undergo a slow hydrolysis that yields alcohol according to the first-order kinetics C e;W ðtÞ ¼ C inj e;W expðÀk Rea ea tÞ where k Rea ea is the kinetic constant of ester hydrolysis producing alcohol.
The corresponding half life time t Rea ea;1=2 ð¼ lnð2Þ=k Rea ea Þ 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 firstorder 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 ).

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 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  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.

Mesh convergence study and sensitivity analysis using a one-dimensional 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 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].

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 Ár n ¼ Ár 0 =n with a refinement level n equal to 1, 2, 4, 8 and 16, where Dr 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 Dt 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 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 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 supp ðC t;W ! max C t;W 2 Þ, using first-and second-order schemes when performing a grid refinement. 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  6. Normalized aqueous tracer outlet concentration profiles simulated using first-and second-order schemes when refining the gridblock size to Ár n ¼ Ár 0 =n with n = 1, 2, 4, 8, 16 and Dr 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 n are reported in Figure 9.

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 5 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 1 Át R Át 0 C a;W ðt À t 0 Þ dt, which would integrate that plateau. Figure 12 and Table 6 compare the residual oil saturation S 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 0 , t e À t 0 and t a À t 0 . 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ð2 m Þj n¼4 case, as shown in Figure 12d, whereas that error   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 second-order 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. is about 35% in the Oð2 m Þj 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 Table 5. 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 supp À C t;W ! max C t;W 2 Á , when using first-and second-order schemes and performing a grid refinement.  Table 6. Comparison of the residual oil saturation S orw computed from the ester and alcohol arrival times, according to equation (33), when using first-and second-order schemes and performing a grid refinement (see also Fig. 12). Expected residual oil saturation is S orw ¼ 0:2. 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.
Scheme n S orw (À) t t À t 0 (days) t e À t 0 (days) t a À t 0 (days) t e À t a (days) 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 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 second-order schemes and performing a grid refinement (see the text).
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 À DC and maxC e,W À DC, with DC set to 50 or 25 ppm. Doing so, ester and alcohol arrival times 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 ½AEAE orw in Table 7 for each grid refinement and scheme order, lead to a maximum residual oil saturation estimate uncertainty range ½S À orw ; S þ orw with S À orw ¼ min S ½AEAE orw and S þ orw ¼ max S ½AEAE orw , and a residual oil saturation estimate uncertainty ÁS orw ¼ S þ orw À S À 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 AE orw as a function of the grid refinement level and scheme order. It is worth noting that the shaded DS 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.

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 Table 7. 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 second-order schemes and performing a grid refinement (see the text 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.