G M Porta

G Ceriotti

J.-F Thovert

Comparative assessment of continuum-scale models of bimolecular reactive transport in porous media under pre-asymptotic conditions

Keywords: reactive transport, continuum-scale models, upscaling, parameter estimation, porescale modeling

Comparative assessment of continuum-scale models of bimolecular reactive transport in porous media under pre-asymptotic conditions

INTRODUCTION

Recent developments in pore-scale numerical modeling and imaging technologies are providing remarkable added value to our knowledge of dynamics of reactive transport processes in complex pore spaces (e.g., [START_REF] Molins | Pore-scale controls on calcite dissolution rates from flow-through laboratory and numerical experiments[END_REF][START_REF] Menke | Dynamic three-dimensional pore-scale imaging of reaction in a carbonate at reservoir conditions[END_REF]. These techniques and approaches allow characterizing and comparing the relative importance of the dynamics associated with advective transport, diffusive mixing and reactive processes at the (micron) scale of individual pores. Reactive transport applications in hydrogeology typically involve large spatial scales. It would then be desirable to explore the extent at which our knowledge on porescale dynamics could be transferred into models that can be applied at a continuum-(or Darcy-) scale. In this work we consider solute transport in porous media in the presence of a homogeneous bimolecular irreversible reaction of the kind A + B → C, which is typically classified as a fluid-fluid reaction. Upscaling this reactive transport problem from the pore-to a continuum-scale is especially challenging under conditions of advection dominated transport and fast reaction, respectively associated with large values of the Péclet, Pe, and Damköhler, Da, numbers (e.g. [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF]Tartakovsky et al., 2009;Porta et al., 2012a). Available column experiments at the laboratory scale document the space-time dynamics of the reaction product C [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF][START_REF] Raje | Experimental study of bimolecular reaction kinetics in porous media[END_REF]. These works show that the assumption of complete mixing of the reactants at the pore-scale can lead to considerable discrepancies between measured concentration of the reaction product C and predictions based on continuum-scale approximations. This result is consistent with theoretical upscaling by [START_REF] Battiato | Applicability regimes for macroscopic models of reactive transport in porous media[END_REF]. Experimental observations documenting the evolution of the pore-scale reaction rate in a Hele-Shaw cell are illustrated by de [START_REF] De Anna | Mixing and reaction kinetics in porous media: An experimental pore-scale quantification[END_REF][START_REF] De Anna | The filamentary structure of mixing fronts and its control on reaction kinetics in porous media flows[END_REF]. These authors show that the spatial distribution of the reaction rates within the pore space is highly non uniform in space and is characterized by a filamentary structure. The pore-scale dynamics of the process of mixing of the reactants, which is inherently linked to the reaction product, has also been investigated through detailed numerical simulations in various settings (e.g., Tartakovsky et al., 2009;[START_REF] Willingham | Using dispersivity values to quantify the effects of pore-scale flow focusing on enhanced reaction along a transverse mixing zone[END_REF]Porta et al. 2012b[START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF][START_REF] Hochstetler | The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime[END_REF][START_REF] Rolle | Effects of compound-specific dilution on transient transport and solute breakthrough: A pore-scale analysis[END_REF][START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF]. These studies emphasize the critical importance of a proper model of local mixing for the interpretation of reactive transport processes. Similar conclusions are drawn from column scale experiments where pH-controlled reactions take place (e.g., [START_REF] Jose | Measurement of mixing-controlled reactive transport in homogeneous porous media and its prediction from conservative tracer test data[END_REF][START_REF] Loyaux-Lawniczak | Acid/base front propagation in saturated porous media: 2D laboratory experiments and modeling[END_REF][START_REF] Edery | Anomalous reactive transport in porous media: Experiments and modeling[END_REF].

Several alternative approaches have been proposed with the aim of embedding effects of pore-scale incomplete mixing within Darcy-scale effective reactive transport models (see also the review by [START_REF] Edery | Reactive transport in disordered media: Role of fluctuations in interpretation of laboratory experiments[END_REF]. These include particle-based Lagrangian methodologies (e.g., [START_REF] Edery | Modeling bimolecular reactions and transport in porous media[END_REF][START_REF] Edery | Particle tracking model of bimolecular reactive transport in porous media[END_REF][START_REF] Ding | Modeling bimolecular reactions and transport in porous media via particle tracking[END_REF] as well as Eulerian approaches (e.g., [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF][START_REF] Chiogna | Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume[END_REF][START_REF] Hochstetler | The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime[END_REF]. Here we focus on the latter set of approaches as applied to bimolecular fluid-fluid reactions. Continuum-scale approaches typically embed effective formulations, which encode incomplete mixing effects into transport and reaction parameters. These effective formulations typically require model calibration against a set of observations, and can be used to interpret available data. In this context, a series of works (e.g., [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF][START_REF] Chiogna | Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume[END_REF][START_REF] Rubio | Numerical solution of the advection-reaction-diffusion equation at different scales[END_REF] are focused on the interpretation of the experimental dataset of [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF].

These authors document the spatial distribution of the reaction product concentration within a fully saturated column where the host porous medium is a bead pack. Longitudinal (i.e., along the mean flow direction of flow) profiles of solute concentration are reported for diverse time levels and transport conditions. A comparison of the results obtained in the above referenced works shows that very different modeling approaches allow interpreting this particular experimental dataset with comparable levels of accuracy. This suggests that the information content associated with the distribution of the reaction product in the system is insufficient for a unique characterization of incomplete mixing through effective continuum-scale approaches. [START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF] perform a three-dimensional pore-scale simulation in an attempt to replicate the reactive experiment scenario of [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF]. They show that spreading of the reactants across the medium can be quantified in terms of a time dependent dispersion coefficient. The temporal evolution of longitudinal spreading is different for the two reactants A and B. However, to the best of our knowledge, the ability of effective models to characterize the concentration of the two reactants (A and B) within the regions where they mix and the reactive process takes place has not been tested to date.

A mathematical formulation relating pore-scale characteristics and effective reactive transport parameters at a continuum-scale can be derived through formal upscaling methodologies, such as the method of volume averaging (e.g., Whitaker, 1999;[START_REF] Orgogozo | Upscaling of transport processes in porous media with biofilms in non-equilibrium conditions[END_REF]Valdes-Parada et al., 2011;[START_REF] Guo | Dispersion in Porous Media with Heterogeneous Nonlinear Reactions[END_REF]. The latter is applied to bimolecular reactive transport by Porta et al. (2012a). These authors show that it is possible to obtain a closed nonlocal formulation of a reaction term associated with incomplete mixing in the presence of of a fast reaction (i.e., large Da values). The assumptions underlying the theoretical analysis of Porta et al. (2012a) have been verified by comparison against direct numerical upscaling of pore-scale simulations performed in simple two-dimensional geometries, i.e., a plane channel (Porta et al., 2012b) and an ordered array of cylinders [START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF]. The continuum-scale model proposed by Porta et al. (2012a) is presented in terms of a nonlocal integro-differential formulation, consistent with a recent formulation obtained by [START_REF] Hansen | First-principles derivation of reactive transport modeling parameters for particle tracking and PDE approaches[END_REF] starting from the generalized master equation. Note that the formulation proposed by Porta et al. (2012) accounts for the temporal dynamics of the spreading process discussed by [START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF] by including a nonlocal in time dispersive transport term.

In this work we start from pore-scale numerical simulations of a bimolecular homogeneous and fast reaction. We consider the latter to take place within a two-dimensional pore space, which is generated by randomly placing in space a set of circular grains of uniform diameter. The key target of the work is to test our ability to interpret the continuum-scale results obtained from direct upscaling of pore-scale simulations in such a disordered system through continuum-scale models. The reactive transport setting we analyze is similar to the one considered, e.g., in [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF], and we model it in terms of an effective onedimensional unsteady transport process at the continuum-scale. We distinguish two different approaches, depending on the type of information employed to compute the associated model parameters: (1) an approach based on the use of effective models which require calibration against reactive transport features that are typically observed at the continuum scale, e.g., spatial distribution or breakthrough curves of reactants; and (2) an approach based on upscaled models which only require information on pore-space geometry and on the velocity field to predict reactive transport at the continuum-scale.

A first goal of this work is to provide the numerical approximation of the nonlocal upscaled model of bimolecular reactive transport derived by Porta et al. (2012a). We then consider three additional effective reactive transport models, i.e., the typical advection dispersion reaction equation, which implicitly assumes complete mixing at pore-scale, and the models proposed by [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF] and [START_REF] Hochstetler | The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime[END_REF]. These models differ in terms of the formulation employed to characterize the effective reaction parameters.

Consistent with typical applications, we estimate these parameters through results obtained from numerical upscaling of pore-scale calculations.

Available experimental datasets typically document the evolution of the reaction product in space and time [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF][START_REF] Raje | Experimental study of bimolecular reaction kinetics in porous media[END_REF]. Therefore, the majority of previous works dealing with this reactive transport setting are solely based on the interpretation of the space-time evolution of the reaction product. In this work we explicitly consider the model performance also in terms of its ability to capture the dynamics of the concentrations of the two reactants A and B. This aspect is key in a pre-asymptotic transport regime, where the characterization of transport can be significantly affected by the initial distribution of the reactants in the domain (e.g., [START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF]. We compare the ability of all these model formulations to interpret the considered reactive transport process and analyze the impact of incorporating diverse types of information in the estimation of model parameters. By doing so, we also identify the critical relevance of specific formats assumed by continuum-scale reactive and transport terms. The results of this work allow quantifying the relevance of nonlocal terms in the reactive transport setting we analyze.

The work is organized as follows. Section 2 provides a description of the pore-scale problem setting. Section 3 is devoted to a presentation of the continuum-scale modeling options we implement. In Section 4 we assess the capabilities of the selected continuum-scale models to interpret the pore-scale simulation results. Concluding remarks end the work.

PROBLEM SETTING

We consider the two-dimensional domain  filled by a fully saturated porous medium, consisting of a liquid ( l  ) and a solid fraction ( s  ). The porous domain is constituted by a collection of periodic unit cells,   . The geometry of the unit cell is generated by the disordered superposition of circular grains of uniform diameter ˆ0.08 w  mm and is characterized by a porosity 0.6  

(see Figure 1). Note that here and in the following all primed variables are defined within a unit cell and hat-signed variables are dimensional. The liquid and the solid phases share an impermeable boundary surface, ls  . Fluid flow through the pore space is associated with a steady laminar two-dimensional velocity field, ˆ() ux , which satisfies the Stokes problem. The flow field is driven by a unit pressure gradient along the x-direction and a zero pressure gradient along the y-direction (Figure 1). The numerical approximation of the velocity field is obtained by the methodology introduced by [START_REF] Bekri | Dissolution of porous media[END_REF] and [START_REF] Coelho | Geometrical and transport properties of random packings of spheres and aspherical particles[END_REF].

A desired value of average fluid velocity, Û , along the x-direction is imposed by multiplying the computed velocity field by an appropriate constant. Figure 1 ) throughout most of the domain, values of u > 3 being mostly concentrated only within a few pore throats.

One can clearly identify the presence of regions which are characterized by small velocity values and are poorly connected to the flow preferential path (e.g., regions highlighted in Figure 1). The dimensionless velocities varies in space over a range of values which is comparable to that obtained by flow simulation performed within millimeter-scale imaged rock samples (see e.g., [START_REF] Bijeljic | Predictions of non-Fickian solute transport in different classes of porous media using direct simulation on pore-scale images[END_REF]. → C takes place in the liquid phase and it is assumed not to affect the velocity field and the pore space geometry. Molecular diffusion is modeled through the standard Fick's law. We employ the following dimensionless formulation of the reactive transport problem [START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF] 2 We consider a replacement setting: the system is initially saturated by a constant concentration 0 ĉ of species A for x > 0; species B is introduced continuously in the system for t > 0 at a constant concentration 0 ĉ through the boundary in  . Equations ( 1)-( 2) are then completed by the following boundary and initial conditions

The external boundary

1 i ii c cc t Pe        u , i D E  (1) 2 1 () B B B B B D c Da c c c c c t Pe Pe          u (2) where 0 ˆîi c c c  , Û  uu , ˆtt U  , ˆ  , m U Pe D  , 2 0 ˆm kc Da D  (3) ˆi c is concentration, 0 ĉ is a reference concentration; k [m 3 /(mol s)] is the reaction kinetic constant; ˆm D is molecular diffusion ˆ is
0 i c    n ,, i B D E  ;
,,

ls imp out     x (4) 1 ( , ) ( , ) e i e i i c t c t N Pe      u n x n x , ,, i B D E  in  x (5) ( ,0) 0 B c  x , ( 0) 1 D c,  x , ( ,0) 0 E c  x ,0 l x   x (6) ( ,0) 1 B c  x , ( 0) 1 D c,  x , ( ,0) 1 E c  x ,0 l x   x (7)
A flux boundary conditions is assumed at in  so that 1

B D E c c c    .
Zero flux is imposed on the liquid-solid interface ls  , and imp  . Note that the boundary in  is located at 60 x  so that the pore space is saturated with a constant concentration of B at 0 t  within the first unit cell adjacent to the inlet section (i.e., for 60 < x < 0). This choice allows avoiding any influence of the inlet boundary condition on the reactive process at early times.

The reactive transport problem is solved within the pore space by means of the particle tracking methodology described by Porta et al. (2012b[START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF] and adapted from routines developed by [START_REF] Bekri | Dissolution of porous media[END_REF] and [START_REF] Debenest | Smouldering in fixed beds of oil shale grains. A threedimensional microscale numerical model[END_REF] for the pore scale description of reactive transport. We select ˆˆ0.08 w  mm and set Pe = 24.88 (i.e., average velocity ˆ0.62 U  mm/s, Reynolds number Re = 0.05) and Da = 1038. This value of Pe is computed by relying on the average fluid velocity and suggests that the process takes place under moderately advection-dominated conditions. We note that the local velocity distribution displays large variations within the considered disordered cell (see Figure 1). The relative strength of the time scales related to transport by advection and diffusion can be evaluated through the local Péclet number [START_REF] Porta | Continuum-scale characterization of solute transport based on pore-scale velocity distributions[END_REF]. One can expect transport to be strongly advection-dominated along the fast channels detected in the system, where the largest values of velocity are about one order of magnitude larger than the average velocity (i.e.,

  ˆ/ loc m Pe D  xu

 

10, 250

loc Pe  ux ), while diffusion drives solute transport within the cavities, where the fluid is basically immobile, i.e.

 

0 loc Pe  x .
The total longitudinal length of the system is L = 1200 ( ˆ9.6 L  cm), i.e. we consider a collection of 20 unit cells   along the x-direction. Our computations rely on a discretization of the pore space into elementary squares K (pixels) of size 0.02 mm. As in [START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF], we consider a unit concentration to be equivalent to a number of particles 100

K NP 
per elementary square. Note that, as a result of these modeling choices, simulating a uniform unit concentration of a single chemical species along the whole computational domain requires a number of approximately 1.7×10 7 particles. This yields reliable local (pixel-scale) concentration values for c > 5  10 -3 [START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF](Porta et al., , 2012b)).

Figure 2 shows a snapshot of the spatial distribution of cA, cB and cC within the pore space for t = 448. Note that reactant A is initially residing in the system. Positive values of cA are observed for x < 500 solely within pores which are scarcely accessible to flow, these areas corresponding to regions (cavities) highlighted in Figure 1. The invading reactant B tends to follow the principal pathways, characterized by large velocity values. As a consequence, local concentrations of reactants exhibit significant and sharp variations across relatively small distances. We also observe that reactant concentrations are anti-correlated, i.e., B c is small in pores where concentration of A is significant and vice versa, due to the occurrence of the fast reaction. This is consistent with previous computational analyses in regular two-dimensional geometries [START_REF] Porta | Numerical investigation of pore and continuum-scale formulations of bimolecular reactive transport in porous media[END_REF]. The concentration of the reaction product C in Figure 2 displays a peak value at x ≈ 450. Contributions to C c in the forward tail (x > 500) are mainly associated with regions characterized by large velocity, stagnant regions chiefly contributing to the backward tail (x < 400).

Pore-scale simulation results are used in the following as calibration and validation data against which outputs of continuum-scale models are compared. For this purpose we consider section-averaged concentrations in the fluid phase. These are defined as

    1 , ii H c x c x y dy H   ;   1 , H G x y dy H    (8)
where H is the total length of the porous medium along y, i c is section-averaged concentration of local concentrations i c of chemical species i = A, B, or C and G(x,y) is an indicator function, i.e. G(x,y) = 0, 1 respectively for solid grains and fluid.

CONTINUUM-SCALE MODELS

A number of approaches have been proposed to model this reactive problem at the Darcyscale. In this work we consider three existing effective models together with the upscaled formulation stemming from volume averaging of the system (1)-( 2), as provided by Porta et al. (2012a). In Section 3.1 we summarize the salient features of the considered effective models. Section 3.2 recalls the definition of the continuum-scale system resulting from volume averaging of the pore-scale system (1)-(2).

Effective models

Effective continuum models of reactive transport are typically formulated through an advection dispersion reaction equation (ADRE)

2 2 1 ( 1) i i i C C C DU t Pe x x           , i D E  (9)     2 * 2 1 1 ˆB B B eff B B D C C C D U K C C C t Pe x x            (10) where ˆ/ m D D D   is the dimensionless longitudinal dispersion coefficient, eff K is an effective reaction constant, i C (i = B, D, E) are continuum-scale concentrations.
We consider in the following three models for Keff proposed in the literature:

 Model 1 (complete mixing): eff Da K Pe  (11)  Model 2 (Sanchez Vila et al., 2010): eff Da K Pe   ; 0 ˆˆm t k     ( 12
)
where 0  and m are parameters to be estimated through model calibration.

 Model 3 [START_REF] Hochstetler | The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime[END_REF]:

eff F Da KE Pe  ; F E Da     ( 13
)
where EF is the reaction effectiveness factor,  and  being model parameters.

Model 1 implies that the same reaction constant can be employed at the pore and continuumscales, i.e. reactants are assumed to be completely mixed at the pore-scale. Models 2 and 3 account for the effect of incomplete mixing of the reactants through the effective reaction parameters ( 12)-( 13). The definitions in ( 12)-( 13 Porta et al. (2012a) perform an upscaling of ( 1)-( 2) through the volume averaging method (Whitaker, 1999). Due to the nature of the scenario we consider, we recall here the corresponding one-dimensional formulation, which is labeled as model 4 in the following. A key target of the volume averaging analysis is to obtain a closure formulation in terms of the intrinsic volume averaged concentration

Upscaled transport model

  1 () l ii V c c dV V    x x ( 14
)
where V(x) is a volume of porous medium within . The volume averaged formulation for the conservative species reads

  1 1 * , l l l l i i i U c c c u D i D E t x Pe x t x                       (15) where   0 * ll t ii UU cc D D t d t x x               ; 1 ls l Ux V D bn dA ub V   (16)
is the dispersion coefficient, l u u u  is a spatial fluctuation (or deviation) of the velocity in the x-direction within the pore space, and b is a closure variable. Equations ( 15)-( 16) are based on the following unsteady closure relationship [START_REF] Moyne | Two-equation model for a diffusive process in porous media using the volume averaging method with an unsteady-state closure[END_REF][START_REF] Chastanet | Mass transfer process in a two-region medium[END_REF]

)   , l l i i i i c c c c b t tx         x (
  2 1 b b u b t Pe        u ,0 l t     x (18) x bn    n ,0 ls t      x (19)   ,0 0 b   x ,0 l t      x (20)
Here, we have considered the unit cell reference system  x , with the corresponding liquid

    12 1 1* l l l l l B B B D U c c c c Da u M D R R t x Pe x x t x Pe                          (21) where   1 l l l B B D R c c c  (22)   2 2 1 l D l c R M M b tx           (
      ll l l l B B D B B D B B D c c c c c c c c c      (24) 2 ll BD B D D ll BD cc c c Mc Mb tx cc         (25)
Substitution of ( 25) in ( 24) leads to

  2 l B B D R c c c  .
We refer to the formulation ( 15)-( 21) as model 4 in Section 4 for simplicity. The parameters appearing in ( 15)-( 21) can be fully characterized from pore-scale geometry and velocity distribution upon solving the system ( 18)-( 20). We do so by implementing a finite element formulation of ( 18)-( 20) in the FreeFEM++ environment [START_REF] Hecht | New development in Freefem++[END_REF]. We discretize the solution in space through a stabilized streamline diffusion finite element formulation while employing a first order implicit backward Euler finite difference scheme in time. We compute the solution by employing a structured grid formed by approximately 10 5 triangular elements.

Figure 3 depicts the time evolution of the coefficients U D (16) and of the quantity

    2 2 , l B t b t   x (26)
We select B2 as a measure of the nonlocal reaction term (23), i.e., the localized counterpart of (23) reads

    2 2, 2 1 l D LOC c R M M B t x          (27)
We observe that both DU and B2 vary over a wide range of values and are always increasing within the considered time window. This result evidences that the reactive transport process we observe takes place within a pre-asymptotic regime. Therefore, the nonlocal terms embedded in ( 15)-( 21) can be expected to play a relevant role in the reactive system characterization.

The numerical solution of the upscaled system requires the discretization of the integrodifferential terms appearing in ( 15)-( 21). We employ a trapezoid integration rule to discretize the integral terms ( 16) and ( 23) in the time domain. System (15)-( 21) is approximated through a fixed space and time discretization. We select a time step t = 2 and a spatial discretization step x = 1. We numerically verified through a grid convergence analysis that the results do not depend on the spatial discretization and time step size (details not shown). We also verify that the results of the nonlocal solution reproduce those of the corresponding local formulation when parameters are constant in time. The computational time associated with a run of the nonlocal model 4 is approximately 10 5 s (28 hours) on a processor Intel(R) Pentium(R) III Xeon processor 2.83 GHz. For comparison, note that the numerical solution of the effective models ( 9)-( 10) requires approximately 10 3 s on the same processor.

COMPARATIVE ASSESSEMENT OF CONTINUUM-SCALE MODELS

In this section we discuss the ability of the continuum-scale models introduced in Section 3 to reproduce the results obtained by pore-scale simulation of the reactive transport process.

As anticipated in the Introduction, we compare two classes of models which are distinguished on the basis of the information sources employed to characterize the model parameters and are associated with diverse levels of complexity. Models 1-3 have a simple structure, but their effective parameters are estimated through model calibration, i.e. they incorporate information on the output concentration profiles. On the other hand, model 4 is characterized by a complex structure, and the related parameters are estimated relying solely on pore-scale information, i.e. no model calibration is required. Our model comparison aims at assessing the impact of the critical differences amongst the considered continuum-scale models on the ability of the models to reproduce our pore-scale reference solution. These key differences can be classified as illustrated in the following. With reference to transport, model 4 is nonlocal in time and allows accounting for non-Fickian dispersion effects, as opposed to models 1-3 which are all based on an assumed Fickian dispersion process. With reference to the way reaction is embedded in the model, model 4 entails a space and time dependent nonlocal formulation of the reaction rate constant; model 2 embeds a time dependent effective reaction parameter Keff through (12); while models 1 and 3 employ a constant effective reaction coefficient, Keff .

In the following, we start by discussing the characterization of the parameters of models 1-3, which is performed through calibration against concentration profiles extracted from the reference pore-scale simulation. We then compare the results of the set of the four continuumscale models illustrated in Section 3 against (section-averaged) pore-scale simulation results.

Calibration of the effective models

The characterization of the effective models (9)-( 13) (models 1-3) is here performed through the estimation of the embedded effective parameters, performed via model calibration.

The model calibration parameters we consider are listed in Table 1 for each of the three models.

The longitudinal dispersion coefficient is considered as a calibration parameter for all three models. This choice is motivated from results of previous studies focusing on the interpretation of the dataset of [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF]. In this context, [START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF] demonstrate through numerical pore-scale simulation that transport in the system is still in a preasymptotic stage. Otherwise, results by [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF] show that a Fickian model can be employed to interpret the same reactive transport dataset up to a reasonable accuracy by including the dispersion coefficient as a parameter to be estimated through model calibration against the reactive transport data. Models 2 and 3 both embed two additional calibration parameters which are associated with the respective reaction models, i.e. 0  and m in (12) for model 2,  and  in (13) for model 3. Given the structure of the reactive term in (13), we note that infinite combinations of ,  can lead to the same effective reaction constant for a given Da. This implies that joint estimation of  and  is possible only when calibration data associated with diverse values of Da are available. Since we consider here a unique value of Da, we follow the indication of Hochstetler and Kitanidis (2013) and set 1   while estimating  .

Our aim is here to study the performance of these models in the presence of different sets of calibration data. We consider two sets of such data, i.e., concentration profiles of the invading reactant B (dataset Bi c  ) and of the reaction product C (dataset Ci c  ) at a fixed dimensionless time,

( 448 t 
). We select this time level because it is the longest time considered in the pore-scale simulation. We then calibrate all three models by employing separately the two datasets. Entries Entries of dataset Bi c  are selected in a similar way, i.e., as the spatial average of the concentration in the liquid phase

  1 1 Li Bi B Li c c dx L       1, ,8 i  (29)
Note that positive values of 31) in a maximum likelihood framework [START_REF] Carrera | Estimation of aquifer parameters under transient and steady state conditions: 1. Maximum likelihood method incorporating prior information[END_REF], where we assume the covariance matrix associated with measurements error to be equal to a diagonal matrix with constant entries. This enables us to compute the covariance matrix of the estimation error 1 2 , , , , This result is consistent with previous analysis [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF][START_REF] Hochstetler | The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime[END_REF] and supports the idea that effective models 2-3 can interpret the distribution of the reaction product in the reactive transport setting we analyze in the presence of pore-scale incomplete mixing of reactants. We also observe that the values of JBmin are similar for models 1 and 2 and JBmin is largest for model 3. We note that values of standard deviation p  associated with the estimates of the dispersion parameters are about 5-10 % of the corresponding ML estimate for all three models and for the two considered datasets, thus resulting in relatively small coefficients of variation. These results suggest that the data convey appropriate information to characterize the longitudinal dispersion coefficients. However, note that the 4b). Previous studies (e.g., [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF]Chiogna and Bellin, 2014) show that the evolution of the peak concentration is typically linked to incomplete mixing of the reactant at pore-scale. We further note from Table 1 that the standard deviations related to the effective reaction parameters embedded in models 2 (i.e, 0  and m) and 3 (i.e, ) are much larger when calibration is performed through dataset Bi c  than by relying on Ci c  . Therefore, results in Figure 4b and4d and Table 1 suggest that an accurate estimation of the incomplete mixing parameters embedded in models 2 and 3 requires information on the reaction product concentration. This result can be explained upon observing that the output

  2 1 1 DC N C Ci C i i DC J c C x N       (30)   2 1 1 DB N B Bi B i i DB J c C x N       (
j j j j T k C k C k C k C    Q J J k = 1, 2, 3; j = B, C (32)

Model comparison

In this Section we aim at assessing the performance of the diverse effective and upscaled models in reproducing the pore-scale simulation results, which represent our reference solution.

We compare in our discussion the results obtained through calibration of the effective models 1-3 against the dataset Ci c  (see Section 4.1) and the prediction yielded by model 4 on the basis of pore-scale information. We focus on a comparison of the longitudinal profiles of concentrations of the three chemical species as rendered by pore-and continuum-scale models and quantify model performance through some global indicators. (Figure 5b,d,f).

As time advances, B displaces A and the two reactants mainly mix in a limited region around the the interface between B and A (i.e., the reaction front), which migrates in the system by advection. Owing to our choice of dimensionless space-time reference, the reaction front is found at locations x ≈ t for each considered time level. The backward tail of A c in our reference pore-scale solution displays an oscillatory behavior characterized by isolated localized peaks (see Figures 5a-b, respectively for x < 70, x < 400). These oscillations are related to the investigated pore structure, i.e. they are chiefly due to the presence of cavities in the porous domain where the solute is trapped for long times (see also Figures 12). As time progresses, the reaction front advances in the porous domain by advection and new immobile (or low velocity) zones with A c > 0 remain isolated and surrounded by the invading solute B (see Figure 2). Thus, the number of localized concentration peaks in the profile of A c increases with time. The effect of these low- velocity zones is also detected on B c and C c profiles (Figures 5c-d ande-f). The section- averaged concentrations of reactants B and A show an anti-correlated behavior, i.e. sharp increases of A concentrations correspond to localized declines of B concentrations, consistent with the pore-scale spatial distributions depicted in Figure 2. The presence of cavities (poorly connected and almost immobile zones) results in localized high concentration peaks that are also visible in the profiles of the reaction product C, which is delayed in low velocity regions. We note that this phenomenon is not observed when considering experimental measurements and numerical simulation of the same reactive process within three-dimensional regular porous media, e.g., glass beads packing [START_REF] Gramling | Reactive transport in porous media: A comparison of model prediction with laboratory visualization[END_REF][START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF]. The considered disordered two-dimensional porous medium includes poorly connected and essentially stagnant regions, giving rise to the local accumulations of species A and C observed in Figure 5.

All continuum-scale models considered in this study fail to reproduce trapping of concentration of species A within immobile regions (see Figure 5b). To quantify the relevance of this inaccuracy, we consider the temporal evolution of A RC , the residual average concentration of A, in a subregion of our computational domain

    0 1 , M X AA M RC t C x t dx X   (33)
XM < L identifying the longitudinal size of a given subdomain we consider. The time evolution of i.e. for t ≈ XM. For t < 300, the average concentration of A progressively decreases to a value approximately equal to 10 -2 . For t > 300, the concentration of A decreases at a lower rate as compared to the above mentioned first stage. In this second stage, the residual mass of A can be found within system cavities and smoothly decreases in time due to reaction. We observe that all four continuum-scale models reproduce almost exactly the (section-averaged) pore-scale results up to a time t ≈ 250, the evolution of RCA in the second stage (t > 300) being misrepresented by all investigated continuum-scale models. Results in Figure 5a-b and 6a suggest that considering (a) effective reaction terms which model incomplete mixing effects (as in models 2-3) and (b) time nonlocalities in both transport and reaction terms (as in model 4) has virtually no effect in improving our ability to reproduce the delayed mass of the reactant A within poorly connected cavities. Figure 6a suggests that all four considered continuum-scale models allow reproducing the residual concentration of the reactant A in the system with comparable accuracy, i.e. up to values of 0.03-0.05 (3-5% of the initial concentration of A).

The reference longitudinal profiles of the concentration of the invading reactant, B c , are well predicted by model 4 for both considered time levels (see Figure 5c-d). The effective models 1-3 yield CB = 1 for x < 200 and t = 448, the solution associated with model 4 tending smoothly to unity at x = 0 and closely following the trend of the pore-scale results (see Figure 5d). Insets in Figure 5c-d representation of mixing at the reactive front, as compared to models 2 and 3, which however provide a reasonably accurate interpretation of the spatial distribution of the concentration of the reaction product C (see Figure 3d).

We consider now the flux weighted concentration of the invading reactant B, corresponding to the solute breakthrough curve (BTCB). The latter is defined as

    , B B M
BTC t C X t  for continuum-scale models and as

  B H B H c udy BTC t udy    (34)
when calculations are performed from pore-scale data. Figure 6b depicts the temporal evolution of the breakthrough curve (BTCB) of B at position x = XM = 300. We observe that model 4 reproduces the pore-scale results more closely than the remaining three effective models. This result is consistent with longitudinal concentration profiles in Figure 5c-d. The observed differences amongst the models can be explained by considering that model 4 accounts for non-Fickian transport through the nonlocal term ( 16), while models 1-3 consider a Fickian dispersion model. Moreover, the volume-averaged (upscaled) reaction term ( 22)-( 23) renders an improved representation of the mixing of the reactants close to the reactive front as compared to the effective reaction terms embedded in models 2 and 3. This is suggested by the close agreement between the results rendered by model 4 and their pore-scale counterparts for low concentration values. As a consequence, Figure 6b shows that, contrary to the results given by the effective models 1-3, consideration of nonlocal transport effects and of the upscaled reactive terms improves our ability to accurately quantify the arrival times of a reactive solute B at a given location in the system. These results are complemented by the temporal evolution of the spatial moments of the reaction product concentration profiles (Figure 7). Figure 7a shows the temporal evolution of the total mass of C in the systems, i.e., the zero order moment of the concentration profile   ( , )

0C C L M t C x t dx    (35)
Figure 7b-c display the evolution of the spreading

     2 0 , CC L C C C x t x x dx t M     ;   0 ( , ) C L C C C x t x dx xt M    (36)
and of the skewness

        3 , C CC C L x x t t C x t dx t           (37)
of the concentration of the reaction product. For the pore-scale simulation, quantities ( 35)-( 37 Figure 7a reveals that effective models 1-3 lead to an overestimation of the total mass of C for early times. These three models yield modest errors for long times (t > 300). This is consistent with the observation that the parameters of models 1-3 are estimated through the concentration of C at t = 448. Model 4 reproduces accurately the total mass of C for early times, while leading to a slight underestimation of the reference result for t > 200. All four continuumscale models underestimate the total reaction product by less than 10% at time t = 448. ). This may be explained by (i) the modeling error embedded in the approximations underlying the volume averaging procedure, and (ii) numerical inaccuracies related to the approximation of the closure problem ( 18)-( 20) and of the integro-differential terms appearing in ( 15)-( 21). With reference to the latter point, we note that while in this work we implement a standard numerical solver, the development of more sophisticated techniques for the treatment of integro-differential terms might improve the global performance of model 4. A detailed analysis of this aspect is beyond the scope of this contribution. The effective models 1-3 underestimate spreading and do not provide a satisfactory interpretation of delayed concentrations of the reaction product C, even as these are explicitly considered in the model calibration procedure (see Figure 5f).

Pore-scale results display some early time oscillations of the skewness coefficient (Figure 7c). The latter attains a maximum value of 0.05 at t = 50. After this time the value of C  associated with pore-scale concentration shows a sustained decrease until it attains negative values for t > 100. This result is consistent with the asymmetric pattern observed in Figure 5e-f.

In particular, Figure 5f shows that for t = 448 the backward tail (x > 450) is characterized by a smaller average (spatial) gradient than the forward one (x < 450). The continuum models 1-3 yield a positive value of C  within the whole considered time window. The non-negligible positive value of the skewness observed at early times for models 1-3 is due to the prescribed set of initial and boundary conditions. We observe that C  tends to zero for long times when we consider the effective models 1-3, i.e. the spatial distribution of C C tends to become symmetric around the peak, consistent with the structure of the implied Fickian dispersion picture (see also . This result is associated with the two features which distinguish model 4 with respect to the effective models 1-3, i.e., (i) the nonlocal nature of the embedded dispersive transport term, and

(ii) the dependence of the reaction term on space and time.

Results of Figure 7 suggest that considering a non-Fickian dispersion model is critical to capture the evolution in time of spreading and of the asymmetry displayed by the longitudinal distribution of the reaction product . This result is consistent with the findings of [START_REF] Edery | Modeling bimolecular reactions and transport in porous media[END_REF][START_REF] Edery | Particle tracking model of bimolecular reactive transport in porous media[END_REF] and [START_REF] Alhashmi | Predictions of dynamic changes in reaction rates as a consequence of incomplete mixing using pore-scale reactive transport modeling on images of porous media[END_REF]. We note that Figure 7a suggests that models based on a Fickian model assumption can interpret the time evolution of the total reaction product mass (corresponding to M0C) with a reasonable accuracy, even under a pre-asymptotic transport regime, as previously shown by [START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF].

Finally, we observe that models 2, 3 and 4 predict very similar peak values of C C at t = 448 (Figure 5e-f), in spite of the very different structure exhibited by the reaction terms in the three models. This is a remarkable result, also considering that previous studies identify the peak concentration of C as a key indicator of incomplete mixing in the considered setting (e.g., [START_REF] Chiogna | Analytical solution for reactive solute transport considering incomplete mixing within a reference elementary volume[END_REF][START_REF] Sanchez-Vila | Interpretation of column experiments of transport of solutes undergoing an irreversible bimolecular reaction using a continuum approximation[END_REF]. We emphasize that the profiles associated with models 2 and 3 are obtained through calibration against data of C c  (see Figure 5f) while model 4 leads to a very similar result directly embedding available information on pore-scale geometry and velocity field, i.e. without the need for any adjustable parameter. The reaction product concentration profile rendered by model 4 shows a change in concavity in the vicinity of the concentration peak for short times (see, e.g., Figure 5e, related to x ≈ 70). This behavior is likely due to the numerical discretization of the integro-differential terms of the model. This observation is also supported by considering that this feature vanishes for long times (see Figure 5f).

CONCLUSIONS

We consider an irreversible homogeneous reaction taking place in a two-dimensional disordered porous medium where pore-scale geometry and velocity are known and study the reactive transport regime characterized by 1 Da Pe

, corresponding to a fast reaction taking place under advection dominated conditions. We simulate the reactive transport process at porescale and we consider the pore-scale numerical results as calibration and validation data for continuum-scale models.

We consider three effective models already presented in the literature which require calibration against observed concentration profiles. Our results show that estimates of effective reaction parameters describing incomplete mixing can be remarkably sensitive to the type of concentrations one employs for model calibration, i.e., parameter estimates can be different depending on whether the concentration of the reaction product or of one of the two reactants is employed. This suggests that none of these models include a robust description of the way the totality of pore-scale processes are transferred to continuum-scale formulations in the porous medium we investigate. The estimated parameter values which model the effect of incomplete mixing in the reaction term are characterized by large uncertainty when we employ the concentration of the injected reactant for model calibration. This result suggests that information on concentration of the reaction product is required for a reliable (continuum-scale) characterization of incomplete mixing in the considered setting.

We numerically solve a continuum-scale nonlocal model resulting from an upscaling of the reactive transport setting considered via volume averaging. This model allows predicting the longitudinal distribution of the reactants and of the reaction product without the need of any adjustable parameter, solely requiring available information on pore-scale geometry and velocity.

This model is computationally intensive, as compared to the remaining three effective models we analyze, but allows interpreting specific features observed from direct averaging of the porescale numerical solution, such as the non-Fickian transport behavior of the invading reactant B and the asymmetric shape of the reaction product concentration profile.

The considered two-dimensional porous system is characterized by a geometry which includes large cavities where the reactant A, initially residing in the system, is trapped for long times. All considered continuum-scale models fail to reproduce residual (normalized) concentrations of A in the system below values of about 3-5%. As a consequence, they are not able to capture the details of the system behavior at large residence times associated with the presence of minute values of trapped mass of solute A. It might be possible that the use of double-or multi-continuum formulations of this reactive transport setting lead to improved interpretation of this specific feature. This aspect will be investigated in future contributions.

ACKNOWLEDGMENT

The financial support of the MIUR Project PRIN 2010/2011 "Hydroelectric energy by osmosis in coastal areas" is gratefully acknowledged. We acknowledge the contribution of Sergey Chaynikov to the pore-scale numerical simulations.

 C  C M 0C

  inward unit vector normal to the boundary. A bimolecular irreversible reaction of the kind A + B

  a characteristic spatial dimension; and Pe and Da are the Péclet and Damköhler numbers, respectively. Note that (1)-(2) are written in terms of the can be readily obtained from algebraic recombination of the mass conservation equations written for the three chemical species A, B and C (e.g.,Porta et al., 2012a and references therein).

  ) are motivated by previous literature approaches to model conservative and reactive transport. Sanchez-Vila et al. (2010) introduce (12) based on the rate-limited mass transfer process observed by[START_REF] Haggerty | What controls the apparent timescales in aquifers and soils? A comparison of experimental results[END_REF].[START_REF] Hochstetler | The behavior of effective rate constants for bimolecular reactions in an asymptotic transport regime[END_REF] introduce the reaction effectiveness factor EF (13) following the definition of the segregation intensity introduced by[START_REF] Kapoor | Bimolecular second-order reactions in spatially varying flows: Segregation induced scale-dependent transformation rates[END_REF].

  D,E) is a spatial fluctuation of the concentration about the mean. The medium we consider is characterized by a periodic structure and the closure variable b can be computed within a unit cell. The following differential problem is satisfied by b

  and the liquid-solid surface ls   . The format of the transport equation for the reactive species B depends on the relative importance of Da as compared to Pe. Here, we consider Da >> Pe so that(Porta et al., 2012a)

  terms (22)-(23) stem from volume averaging of the pore-scale reaction term

  of the calibration datasets are concentrations Bi c  (i = 1 … DB N ), and Ci c  (i = 1 … DC N ), DB N and DC N indicating the number of data included in the datasets Bi c  and Ci c  , respectively. of the unit cell   along the x-direction. In other words, we consider the 20 unit cells constituting the total length L of the porous domain and identify Ci c  as the spatial average of the pore-scale concentrations of the reaction product cC computed within each of the unit cells (i.e, we set 20 DC N  ).

Bc

  are observed only for x < 500, i.e., in the first 8 unit cells. We augment the entries in this calibration dataset by considering additional concentration values at locations close to the reaction front (x ≈ 450), because one of our goals is the characterization of the parameters of the incomplete mixing models which are expected to influence concentration values at such locations. As such, dataset Bi c  includes a total of 16 data which are distributed along the length of the porous domain

J

  where indices k and j respectively indicate the model and the dataset employed for is the Jacobian matrix whose entries are the derivatives of the output state variables Cj with respect to the parameters of model k. The diagonal entries of matrix (32) quantify the uncertainty 2 p  associated with the ML estimate of model parameter p.

  each of the three effective models (9)-(13). The calibrated model results are compared in Figure 4 against the calibration data and the complete profile of section-averaged concentrations   B cx and   C cx rendered by pore-scale simulation at t = 448. Models 2 and 3 lead to smaller values of JCmin than model 1.

Cic

  data are employed. In general, we observe that the estimated values of the effective reaction parameters embedded in models 2 and 3 largely depend on the type of information available for calibration.

Figure

  Figure4dshows that the peak of concentration

  concentration B C rendered by the investigated continuum models attains very small values close to the reaction front. Large variations of the incomplete mixing parameters induce modest variations of B C . Otherwise, the peak concentration of C C is very sensitive to variations of these reaction model parameters, as previously shown by Sanchez-Vila et al. (2010), Chiogna and Bellin (2013), and Ciriello et al. (2015). Therefore, model calibration against Ci c  yields a reduced uncertainty for the reaction parameters which represent incomplete mixing, as compared to results obtained through calibration based on Bi c  observations.

Figures

  Figures 5a-f depict a comparison between the concentration profiles rendered by the

AFigure

  Figure 6a depicts the evolution of

  depict profiles of the concentration of the B reactant in semi logarithmic scale. Model 4 reproduces closely the forward tail of the invading reactant, which is mostly affected by the reactive process. Values of B c at locations where 0underestimated by model 1. This result suggests that the volume averaged formulation (15)-(21) (i.e., model 4) leads to an improved continuum-scale

Figure

  Figure 5e-f depicts the comparison between the longitudinal profiles of concentration of

Figure

  Figure 7b depicts the temporal evolution of the spreading of the reaction product C, as

Figure 5f )

 5f Figure 5f). Otherwise, results of the nonlocal model 4 are characterized by a negative skewness

Figure 1 .Figure 2 .

 12 Figure 1. Distribution of velocity / U  uu within a unit cell   . Areas enclosed by the white

Figure 3 .

 3 Figure 3. Temporal evolution of the effective dispersion coefficient DU (16) (continuous curve,

Figure 4 .Figure 5 .

 45 Figure 4. Model calibration results: concentration profiles resulting from model calibration

Figure 6 .Figure 7 .

 67 Figure 6. Temporal evolution of (a) flux-weighted concentration of reactant B (BTCB) at XM =

  Table 1 lists the values of the estimated parameters and the related standard deviation

p

 , together with the minimum values of (

30

)-(

31

) (respectively denoted as JCmin and JBmin) obtained through calibration against the two types of data

  estimated longitudinal dispersion coefficients are significantly different across models and/or depending on the type pf calibration data employed. Model calibrations based solely on Bi c  yield estimated longitudinal dispersions which are considerably larger than those based on

			c  (see
			Ci
	Table 1). This leads to an increased spreading of the reaction product when these models are
	calibrated against	c  as opposed to what can be obtained by employing	c  as a calibration
		Bi	Ci
	dataset (see Figure 4b and d).
	The estimate of parameter m in model 2 rendered by calibration against

Ci c  is consistent with the value obtained by Sanchez-Vila et al. (2010) by using the same type of information.

Note that Table 1 lists 0  in dimensional units [m 3 mol -1 s m-1 ]. Given that the dimensional units of this parameter depend on the specific value assumed by the exponent m, the values of 0  obtained in Table 1 are hardly comparable to those obtained in Sanchez-Vila et al. (2010). The best estimate obtained for parameter  in model 3 is of the same order of magnitude of the one obtained by Hochstestler and Kitanidis (2013) when

Data


Hermine Chatoux

Noël Richard

François Lecellier

Chrisitne Fernandez-Maloigne

Can fixation-point and key-point coincide on cultural heritage colour paintings?

This article compares the key-point extracted with a colour key-point detector and the location of fixation point thanks to an eye-tracking experiment. We hypothesize the first fixation points should correspond to key-points with the most significant gradients. The colour detector is based on Harris and Stephens corner detector extended to colour. The eye-tracking experiment was realised on medieval art work. We compare the location of both the detected key-points and the fixation points observed. Between 30% to 50% of the key-points coincide with fixation points. A second analysis display the number of matching key-points over the number of fixation points at a given time. The ratio decreases with the observation time which concur our initial hypothesis. Overall, several of the first fixation points correlate with high response key-points detected with our method.

Introduction

Several eye-tracking experiment have been conducted to estimate saliency maps. [Borji and Itti(2012)] offered a wide survey of the different methods used to compare an estimated saliency map. The estimation can focused on top-down analysis (task related) or bottom-up exploration (pre-attentive step). [START_REF] Borji | Computational modeling of top-down visual attention in interactive environments[END_REF] Borji, Sihite, and Itti] proposed classification tasks to extract top-down saliency maps. On the contrary [Le Meur(2005)] used low-level features to perform bottom-up saliency maps. These maps are based on the location of fixation points: a point focused by a human eye for sufficiently long time.

It is important to define what represent a key-point: it is a point or an area allowing to characterize the analysed image. Two major family extract key-point, the corner point detection and the blob detection. The most used detector is the Difference of Gaussian (DoG) presented by [START_REF] Lowe | [END_REF]] belongs the the blob family. It is based on the second derivative and extract uniform areas. [START_REF] Bay | Speeded-Up Robust Features (SURF)[END_REF] Bay, Ess, Tuytelaars, and Van Gool] proposed an adaptation that was computationally faster. The corner detection family can be decomposed in two group: one based on the first derivative and based on pixel comparison. The derivative based group was introduced by [Moravec(1980)] and generalized by [Harris and Stephens(1988)]. They rely on the spatial auto-correlation matrix and study the eigenvalue to determine which pixels are key-points. The second group compares a value set as the centre with values at a given radius from the centrer. It was introduced by [Smith and Brady(1997)]. [START_REF] Rosten | [END_REF]] accelerated the computation with the FAST detector which is often used.

In this article, we wonder about a possible relation between detected key-points and salient fixation point. Key-point detection is based on low-level features. Therefore we expect a better match between key-points and fixation point in the pre-attentive phase. The rest of the observation is linked to the brain analysis reading of the image. So, low-level features are not the main criteria. However, our goal is not to prove both detected key-points and fixation point coincide. But we are interested in estimating what matching we can expect.

The next section will present the two colour key-point detectors used in the comparison. Section 3 presents the eye-tracking experiment: device, images and observers group. The fourth section compares the detected key-points with the fixation points obtained during the experiment. We conclude in the last section.

Colour key-point detector

This vector-key-point detector is based on the same steps as the [Harris and Stephens(1988)] corner detector. The figure 1 presents these steps which are separated in three phases. First, the gradient needs to be measured, then the corner informations are extracted and finally the key-point decision is made.

Gradient extraction

Shaping

Decision criteria

Image

Multi-scale derivative filters Auto-correlation matrix:

M SC ∂x T 1 G∂x 1 ∂x T 1 G∂x 2 ∂x T 1 G∂x 2 ∂x T 2 G∂x 2 Gradient rate extraction det(M SC ) tr(M SC )
Local maxima extraction Key-point relocation

Key-point

Figure 1: Steps to obtain full-vector key-point

Gradient extraction

As an image contains corners at different scale, we used multi-scale derivatives filters. [Mikolajczyk and Schmid(2 proposed multi-scale binomial derivative filter. To be more generic, we construct our filter with a Gaussian and its derivative. Equation (1) illustrates these filters:

∂F ilter ∂x i =    x y    ×    x y    T .
(1)

The Gaussian function is centred, the spatial filter size (S F ) depends only on the standard deviation σ selected: S F = (6σ + 1) × (6σ + 1).

We selected eight following scales form σ = 1 to 8. Usually, the scales are selected with a constant size ratio. We chose a constant six pixels difference between each scale. It allows a redundancy in certain key-point location. The more scale detected at a key-point, the stronger it is. These filters allows to measure, marginally, the gradient on every channel. The next step is to combine these gradients.

Inspired form [Di Zenzo(1986), Koschan and Abidi(2005)], we proposed the full-vector gradient in [Chatoux et al.(2019a)Chatoux, Richard, Lecellier, and Fernandez-Maloigne] to extract gradients in the sensor space. We calculate a spatial correlation matrix (M SC ) depending on all the hyper-spectral channels and their inter-correlations from the Jacobian:

M SC =   ∂I(x) ∂x 1 2 ∂I(x) ∂x 1 , ∂I(x) ∂x 2 ∂I(x) ∂x 1 , ∂I(x) ∂x 2 ∂I(x) ∂x 2 2   , (2)
with

∂I(x) ∂x 1 , ∂I(x) ∂x 2 = ∂I(x) ∂x 1 T • G • ∂I(x) ∂x 2 , (3)
G =      s 0 2 2 s 0 , s 1 2 . . . s 0 , s m 2 s 1 , s 0 2 s 1 2 2 . . . s 1 , s m 2 . . . . . . . . . . . . s m , s 0 2 . . . s m , s m-1 2 s m 2 2      . ( 4
)
The Gram matrix G uses the scalar products defined for the integrable functions. The functions used are the Spectral Sensitivity Functions (SSF) of each channel:

s i , s j 2 = R S i (λ)S j (λ)dλ. (5)
The correlation matrix will allows us to extract key-point as presented by [Harris and Stephens(1988)].

Shaping

An intermediate step is necessary to extract the strongest gradients represented by corners. The gradient rate extraction consist in a simple integrative filter. It is often considered as a denoising step in other algorithms ( [Harris and Stephens(1988), [START_REF] Mikolajczyk | [END_REF]). With the filters we proposed, the smoothing step is already realized. Yet, an integrative filter adds all the gradients from its window. A corner area contains more edges than an edge one, therefore after the integrative filter, the response from the corner area will be stronger than the edge one. This step allows to increase the response on corner areas.

For the appropriate response function, we based our analysis on the principle given by [Harris and Stephens(1988)]:

• 2 small eigenvalues represents uniform area,

• 1 strong eigenvalue represents edge,

• 2 strong eigenvalues imply a corner area.

The response proposed by [Harris and Stephens(1988)] depends on both matrix invariants: the determinant and the trace. The proposed response function is

R H (k) = det(M SC ) -k tr(M SC ) 2 , k being an empirical constant.
To free ourselves from the instability implied by the constant k, we propose a new response:

R F V KP = det(M SC ) tr(M SC ) . ( 6
)

Decision criteria

Corner areas will give high value on the response function. To select only the corner, we extract the local maxima of high responses. These maxima are labelled corners. Once the key-point have been selected, several informations are attached to it apart from its location: the response value and an angle.

The spatial direction θ of the gradient is defined by [START_REF]Improved direction estimation for di zenzo's multichannel image gradient operator[END_REF]Jin, Liu, Xu, and Song] lifting the imprecision of ± π 2 from the initial Di Zenzo expression. Another step can be added to merge key-points. A corner can be detected at several filter sizes. It is not relevant to keep each size for the same location. Therefore, we propose two different extractions of the key-point.

The first one is based on suppressing overlapping key-points. We keep the key-point with the strongest response function when there is a large overlap between two key-points. It will be called Harris S 2 on the rest of the paper.

The second detector consider the set of key-points associated to the same location with different scales. We hypothesize that the more scale are detected the stronger the corner is. This detector, Harris M E , will keep a key-point if it has been detected on 3 scale or more and associate the number of detected scales as the response function. The figure 2 presents the results of both detector on the same image. We observe that the detector Harris M E is more selective than the other. The constraint of several scale of detection allows in this image to remove all key-point associated to a mosaic tile.

Eye-tracking experiment

To compare the detectors with the human eye, we realized a psycho-visual experiment allowing to eye-track the user during the reading of an image.

The eye-tracker An eye-tracker and its associated software analysed the fixation point of each observer. An eye-tracker is a device that allow to follow the eye movement of a user. The device measured the eye movement observing a screen. It is based on infra-red that will be reflected by pupils only.

The eye-tracking device used id the Tobii X-120. It allows to measure two majors eyemovements. Firstly, there is the saccades and then fixation points. The saccades are high speed movements (<50 ms pause between movement). There is no apparent pattern to their trajectory. These saccades allows to create a first image that will allows the brain to extract the area of interest that will be modelled by the fixation points. These points are obtained when the eye is still for a sufficiently long time (>300 ms). It allows the brain to analyse the area. The fixation points are the one we want to compare with key-points detected with the detectors. The eye-tracking device needs to be calibrated to precisely localise where the observer is looking. The Figure 3 presents the experimental diagram and a photo of the settings during the experiment.

As the experiment was to be held in two different environments, we chose to display the images on a 17 inches screen. The eye-tracking devices used was created for larger screens. Therefore, observers were watching small images considering the perception distance forbidding a great sensitivity in the details.

The images We chose to set the test duration short (12 to 15 min) in order to keep the observers attention focused during all the experiment. We selected 20 images of medieval works.

Among these, half of it came from the Romane database created by the CESCM (Centre d'Études Supérieur de la Civilisation Médiévale) where come most of our observers. These images are defined as known by the observers.

Ten other images have been extracted on other medievals work databases. These are considered unknown for the "experts".

Most of the images are wall paintings but there are also mosaics, wood paintings. . . Some of them present a religious theme but other themes have been selected. Images are well preserved other are very degraded. The figure 4 presents some of the used images.

Each images is presented for 30 secondes with 5 secondes of neutral grey in between to rest the eyes. The oberver has to describe the seen image to identify the scene when possible.

Sample group One of the experiment objectives was to verify how different the analysis is done depending on the knowledge of the observer. We have selected two groups to realise the experiment.

The first one is constituted of students from BD to PhD, researchers specialised in the medieval history. 21 observers are in this group.

The second are members of the XLIM laboratory (researchers, administrative or IT personnels). This sample group is unfamiliar with medieval images even if some are image processing users. 16 persons realised the experiment in this group.

Unfortunately some users cannot be kept in the study. Indeed, glasses or moving during the experiment prevented the eye-tracker device to measure sufficiently the eye movement. We kept 31 observers.

Comparing detectors key-points and saliency points

In this study, we are looking for a relation between detectors key-points and saliency points. We have extracted key-points with the colour detectors Harris S 2 and Harris M E on some of the images

First results and works hypothesis

The figures 5 presents results on fixation points for an observer compared to the detected keypoints for two images. On the first image (fig. 5a), points seems dispersed on the central area (coloured apse). On this image other observers focused less on the frieze inducing larger difference with the detector (fig. 5b). On the second observed image (fig. 5c), fixation points are focused in the remaining painted area avoiding the degraded areas while the detectors (fig. 5d) extracts mostly in the edges between painted and plaster areas as they present strong gradients. This shows a first limits in the analysis of a relation between detected key-points ans salient points. A strong gradient does not necessarily induce high saliency. The number of key-point being superior with the other detector, the percentage of key-point on an edge painting/ plaster is smaller hence the correlation with saliency point is better.

It is worth noting, that the saliency central bias is not considered in this study. The detected key-points are extracted on the whole image while human looks preferably at the centre of the image. Considering the observation duration, the observer should have time to watch the whole image.

The table 1 gives a first results of the relation between detected key-points and fixation points. It presents the percentage of detected key-points located near (less than 10 pixels radius) of a fixation point. Overall 35% of the detected key-points are scrutinize by observers. It gives us a first limits to the expected relation, even if four images are insufficient to draw a conclusion.

The figure 6a presents the number of fixation points over time. These curves corresponds the visual course presents in figures 5a and 5c. With the description task given to the observer they observed the scene with an almost constant frequency of 3 to 4 fixation points per second. To compare this visual course to the extracted key-points, we hypothesize the observation order can vary from and observer to another but a similar set of fixation point is observed after a certain duration. We need to consider the different level of perception from pre-attentive vision to semantic driven vision. Therefore, only a part of key-point is observed by humans. We expect a theoretical curve as presented in figure 6b.

Correlation between detected key-points and fixation points

The curves from figures 7a and 7b shows the our initial hypothesis (expected results from figure 6b) seems consistent. Nevertheless, these curves does not allow to define the pre-attentive duration. As the detector Harris S 2 is less selective, it offers higher correspondence rate between key-points and fixation points. Nonetheless, if some key-points are observed by several observers some are never observed.

The figures 7c and 7d present the ratio of observed key-points over all the observers fixation points depending on the time ordered fixation points. The curves decreases with time. The decrease vary in speed and tends toward zero as the number of fixation points increases indefinitely with time while the number of key-points is finite. It reinforces our hypothesis that mainly the Table 1: Numbers of key-points detected by our proposition. The matching percentage is the number of key-point looked by an observer divided by the total key-point number.

Harris pre-attentive phase will match with detected key-points.

The previous figures showed the first fixation points have a higher match rate than the later ones. We can ask ourselves, if the matching key-points corresponds to the stronger ones. The figures 8a and 8b presents the sum of the response functions of the matching key-points depending on the time ordered fixation points. If a relation is difficult to extract form the detector Harris M E , it is clear for Harris S 2 : stronger key-points are observed on the pre-attentive phase.

From another point of view, figures 8c and 8d present the number of observation of a keypoint depending on the key-point ordered by decreasing response functions. For the detector Harris M E , the number of observation is overall constant whatever the response function. On the contrary, for the detector Harris S 2 , the higher response functions key-point have a slightly higher probability of being observed by humans. We cannot extract a definite tendency with these two figures. It corroborates our initial caution when comparing key-points and fixation points, After the pre-attentive phase, the analysis is brain driven to extract the semantic meaning that was pronounced during the experiment.

Conclusion

Overall, several of the first fixation points correlate with high response key-points detected with our method. This comparison supports our conjecture based on strong correlation between the first fixation points corresponding to the pre-brain analysis and the corner key-points extracted with our detector.

The comparison should be more developed to better understand the link between fixation points and key-points detection. In fine, this could guide the future key-point detector to a better harmony with the visual perception and brain analysis.

The fixation points analysed were obtained by looking the images for a long time. Moreover, the task was to identify the scene when possible which is a top-down approach. As said in the introduction key-points are low-level features not immediately related to these approaches.

Comparing key-points and fixation points associated to another task such as looking for details could induce a better matching rate.

Another limit of this analysis is that the key-points are selected to be corners while fixation points can be a corner or the centre of a uniform area. Therefore, this study should be completed with a one comparing a blob detector and the fixation point. This can be a response on the battle corner/blob detection: both are of interest in regard to human vision!

  (a) 10522 points extracted by the detector HarrisS 2 (b) 687 points extracted by the detector HarrisME

Figure 2 :

 2 Figure 2: Examples of key-points extracted by both detectors

Figure 3 :

 3 Figure 3: Experimental protocol for eye-tracking on medieval works.

  (a) Known image of wall painting (b) Known image of deteriorated wall painting (c) Unknown image of mosaic (d) Unknown image of wall painting

Figure 4 :

 4 Figure 4: Examples of images used for the experiment.

  (a) Example of an observer fixation points (b) 657 points extracted by the detector HarrisS 2 (c) Example of an observer fixation points (d) 163 points extracted by the detector HarrisME

Figure 5 :

 5 Figure 5: Examples of extracted key-points for an observer and the key-point detector.

  fig.5afig.5c

Figure 6 :

 6 Figure 6: Preliminary results. The figure 6a presents the relation between the number of fixation point and time. The second one (fig. 6b) presents the expected relation between key-points and fixation point.

  fig. 4a fig. 4b fig. 4c fig. 4d

Figure 7 :

 7 Figure 7: Match between key-point and saliency depending on the detector. The figures 7a and 7b present the number of key-point observed depending on the time (fixation points number × observers). Figures 7c and 7d present the number of key-point observed over the number of fixation point.

  fig. 4a fig. 4b fig. 4c fig. 4d

Figure 8 :

 8 Figure 8: Harmony between key-point and saliency points depending on the intensity of the point. The figures 8a and 8b present the sum of the response function depending on the chronologically ordered fixation point. Figures 8c and 8d present the number of observation ordered by decreasing response function.