Fabrice Golfier

Cesar Zarcone

Brigitte Bazin

Ronald Lenormand

Didier Lasseux

Michel Quintard

On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Acid injection is a process widely used for stimulation of petroleum wells to increase rock permeability [START_REF] Rowan | Theory of acid treatment of limestone formations[END_REF][START_REF] Williams | Acidizing Fundamental[END_REF][START_REF] Chang | A novel self-diverting acid. Developed for matrix stimulation of carbonate reservoirs[END_REF]. The dissolution of the porous matrix is an unstable process similar to viscous fingering that leads to the formation of dissolved channels called wormholes. When a part of the medium is etched, the local permeability is increased and consequently the mass transfer is also increased, leading to more dissolution. The result is the formation of ramified channels with typical sizes ranging from the pore size to a few mm. Many experimental and theoretical studies have attempted to describe the roles of the various factors on the formation of the wormholes: injection rate, acid volume, permeability, reaction kinetics, etc. Today, most of the mechanisms can be considered as understood but the models available in the literature are still qualitative. The purpose of this paper is to propose a numerical model for predicting the development of wormholes at the scale of a core sample.

The dissolution kinetics at the level of the solid walls of the grains is well understood. The dissolution process can be separated into three successive steps: acid transport by diffusion and advection to the solid surface, chemical reaction at the solid surface, and transport of the products of the reaction away from the surface. If the chemical reaction characteristic time is very short compared to the mass-transfer kinetics, the process is called mass-transfer limited. This is the case for limestone dissolution with HCl at temperatures greater than 0 • C [START_REF] Lund | The dissolution of calcite in hydrochloric acid[END_REF]. On the other hand, if mass-transfer kinetics is slow, then the process is reaction-rate limited. It is the case for dolomite dissolution with HCl at temperatures less than 50 • C [START_REF] Lund | The dissolution of dolomite in hydrochloric acid[END_REF]. Both mass-transfer and reaction-rate-limited mechanisms may lead to Darcy-scale non-equilibrium dissolution, which means that acid is not immediately consumed as soon as it enters into the porous medium. The dissolution zone within the porous domain is characterized by a variable porosity. On the other hand, local-equilibrium dissolution corresponds to immediate consumption of the acid at the fluid-porous medium interface, and this leads to a sharp interface between the fluid region and the porous medium. In this paper we will study at the Darcy-scale local-equilibrium and non-equilibrium dissolution by considering only mass-transferlimited dissolution processes. This mechanism has been studied both experimentally and theoretically.

Several experiments on dissolution have been performed in a variety of fluidmineral systems, and for mass-transfer-limited or reaction-rate-limited processes. Fredd, Hoefner and Fogler [START_REF] Hoefner | Pore evolution and channel formation during flow and reaction in porous media[END_REF]Fredd & Fogler 1998a, 1999) used limestone core samples and HCl injection. The wormhole structure is visualized by injecting a low-melting-point alloy and then dissolving the porous medium or by neutron radiography. In addition to HCl-limestone systems, Daccord, Lenormand & Lietard (1993a) also studied the water-plaster-of-Paris system in various geometries (disks and cylinders). [START_REF] Wang | The optimum injection rate for matrix acidizing of carbonate formations[END_REF] investigated the effect of temperature, acid concentration, rock mineralogy and injection rate for some Indiana or Glenn Rose limestone and dolomite cores. Bazin and co-workers [START_REF] Bazin | Experimental investigation of some properties of emulsified acid systems for stimulation of carbonate formations[END_REF][START_REF] Bazin | From matrix acidizing to acid fracturing: A laboratory evaluation of acid/rock interactions[END_REF] have conducted acidizing experiments using both limestone and dolomite and various acid fluids including acids in emulsion. All the experiments led to the same conclusion of different dissolution regimes with compact dissolution, ramified wormholes and uniform dissolution. A more detailed description of the regimes will be given with the presentation of our experimental data. Another important result is the determination of an optimum injection rate [START_REF] Hoefner | Pore evolution and channel formation during flow and reaction in porous media[END_REF][START_REF] Wang | The optimum injection rate for matrix acidizing of carbonate formations[END_REF][START_REF] Bazin | A laboratory evaluation of acid propagation in relation to acid fracturing: Results and interpretation[END_REF]. It corresponds to the maximum of penetration of the wormhole for a given volume of acid injected. It is reported that the optimum conditions are related to the formation of a dominant wormhole with little branching through the core. The optimum flow rate depends on several parameters, including rock mineralogy (calcite or dolomite for carbonate formations), temperature and acid concentration. The influence of the various parameters on the acid injection process has been studied in detail by [START_REF] Bazin | Experimental investigation of some properties of emulsified acid systems for stimulation of carbonate formations[END_REF].

Several models have been developed to represent the wormhole formation and growth:

(i) The simplest approach assumes that a cylindrical wormhole already exists [START_REF] Hung | A mechanistic model of wormhole growth in carbonate matrix acidizing and acid fracturing[END_REF][START_REF] Wang | The optimum injection rate for matrix acidizing of carbonate formations[END_REF][START_REF] Buijse | Mechanisms of wormholing in carbonate acidizing[END_REF][START_REF] Huang | Reaction rate and fluid loss: The keys to wormhole initiation and propagation in carbonate acidizing[END_REF][START_REF] Huang | Prediction of wormhole population density in carbonate matrix acidizing[END_REF]. In the first models, the solid is not porous and dissolution of the wall of the tube is studied by solving the transport equation inside the tube. The velocity field (radial and longitudinal component C(r, x)) is calculated by using Stokes equations inside the tube (Poiseuille flow). The concentration field is determined with the real dissolution kinetics condition at the wall (zero concentration in the case of a mass-transfer-limited process) and introducing the molecular diffusion in the liquid phase. This simple model gives a good understanding of the dimensionless numbers that govern the dissolution. The other interest of the model is to illustrate the notion of upscaling. If the concentration is averaged along the longitudinal axis, the average concentration is no longer in equilibrium with the 'wall' concentration, and a kinetic term appear in the one-dimensional transfer equation. Improvements have been introduced by considering porous walls (fluid loss) and a wormhole population density. However this model is based on a pre-existing shape of the wormhole and cannot describe the ramified structure observed in the experiments.

(ii) The network approach considers the dissolution of the grains of the solid in a two-or three-dimensional structure. Generally, the space between two pores is represented by a cylindrical capillary. Contrary to the previous model, velocity and concentration are not solved by the Stokes equation but replaced by the averaged properties. The average flow rate is assumed to be proportional to the pressure difference (Poiseuille's law) and the dissolution proportional to the average velocity and average concentration [START_REF] Hoefner | Pore evolution and channel formation during flow and reaction in porous media[END_REF][START_REF] Daccord | Carbonate acidizing: Toward a quantitative study of the wormholing phenomenon[END_REF]. Network simulations are able to reproduce, at least qualitatively, the trends observed experimentally for the dissolution regimes and obtain an optimum injection rate. The major difficulty with the network approach lies in the application to largescale systems, because of the limited number of pores considered in the simulators. In addition to the problem of size, the network approach does not account for the pore merging, and this may lead to a lack of representativity of the model in the fluid zone. This may explain some discrepancy observed in comparing experimental data and simulation results, especially for the optimum flow rate as reported by Fredd & Fogler (1998b). Attempts to obtain a more physically representative network have been undertaken (Fredd & Fogler 1998b). While a better representation of the geometry and the physics was achieved, this kind of model still requires an approximate solution for the evolution of the geometry as well as for the solution of the Navier-Stokes equations. Moreover, this has been developed so far for packed beds only.

(iii) A different approach, used in fluid mechanics or thermal processes, is based on correlations at the scale of the core sample itself (Daccord, Lenormand & Lietard 1993a, b;[START_REF] Frick | Modeling of fractal patterns in matrix acidizing and their impact on well performance[END_REF]. All the physical parameters involved in the problem are introduced in the form of dimensionless numbers. The relationships between observed parameters and dimensionless numbers are determined either from experiments or using theoretical considerations. For instance, Daccord et al. (1993a) used the properties of self-similarity of the structure of the wormholes in the ramified regime to write the scale dependence as a power law characterized by the fractal dimension. These approaches have been partially validated for experiments on the water-gypsum system. Although these works give some qualitative results (dissolution behaviour diagram), they cannot be used as fully predictive tools. A similar approach was used by [START_REF] Hoefner | Pore evolution and channel formation during flow and reaction in porous media[END_REF], [START_REF] Fredd | Optimum conditions for wormhole formation in carbonate porous media: Influence of transport and reaction[END_REF] and [START_REF] Fredd | Dynamic model of wormhole formation demonstrates conditions for effective skin reduction during carbonate matrix acidizing[END_REF] to estimate a Damköhler number corresponding to the optimum flow rate. It reproduces the trends observed experimentally for the dissolution regimes, but the optimum injection rate cannot be predicted without a model of wormhole growth. Moreover, the concentration effect on the wormhole development is not taken into account.

(iv) A model based on continuum equations written at Darcy's scale has been developed by [START_REF] Lhuillier | La viscosité effective dans le modèle de Darcy-Brinkman: Approche self-consistante. 4èmes Journées d'Etude sur les Milieux Poreux de Nancy[END_REF] and [START_REF] Chen | Reactiontransport simulation of matrix acidizing and optimal acidizing strategies[END_REF]. The resulting simulator solves the Darcy-scale equations of fluid flow, acid transport and solid-fluid reaction. Although fingering and wormholing phenomena were captured, a complete description of the different dissolution regimes was not provided and prediction for optimum injection conditions was lacking.

All these efforts have led to the understanding of the physics of porous media dissolution and wormhole development. However, most of the existing models fail in describing quantitatively the coupled nature of flow and reaction without assuming a wormhole geometry, or, otherwise, are limited to porous domain of small size or of particular shape (packed bed structure for instance). In addition, the use of Darcy-scale models has not led to a thorough investigation of all dissolution patterns. The purpose of this paper is to build a Darcy-scale numerical simulator based on improved physics and test it with existing or new experiments that will be presented.

The first part of this paper describes the physical model used to account for the coupling between flow and dissolution. In the second part the numerical implementation of our model is presented. The numerical schemes used to solve the set of equations, boundary and initial conditions are detailed. Finally, results are presented and compared to experimental results from porous media made of salt, in which solutions of under-saturated salt water were injected at a constant flow rate. The main result is a good agreement in terms of dissolution regimes and prediction of wormhole evolution.

Physical model

The different scales of our problem are represented in figure 1 from the pore scale to the core scale. In this paper our interest is focused on the intermediate scale, or Darcy scale (also identified as the local scale) where the properties are averaged over a large number of pores. Since the development of a Darcy-scale dissolution model is not straightforward, we need to start with the pore-scale description in order to understand the limits and the problems associated with Darcy-scale continuum description.

Pore-scale equations

At the pore-scale, the solid phase and the fluid phase are identified as the σ-phase and the β-phase respectively. In the following equations, A βσ represents the area of the interface between the two phases and n βσ the unit normal vector directed from the β-phase towards the σ-phase. We denote respectively as c Aβ and c σβ the mass concentration of the acid species and the solid species in the β-phase, and c Bβ is the water mass concentration.

We consider that the concentrations of the diffusing species in the β-phase are small, and that the acid (or dissolved rock) transport in water does not change the liquid density ρ β and viscosity µ. Therefore, the flow velocity can be solved independently from the reaction-transport problem (for a given porosity field, ε β ). Furthermore, inertial terms are neglected in the momentum equations and we use Stokes' equations. We denote as v Aβ and v β the velocity of the acid species and the average velocity in the β-phase, respectively.

In the general case, the overall reaction between acid and rock can be expressed by the following boundary condition at the interface A βσ :

-D∇c Aβ • n βσ = k c (c Aβ ) n at A βσ , (2.1)
where D is the molecular diffusion coefficient, k c is the reaction rate constant and n is the order of the chemical reaction [START_REF] Fogler | Elements of Chemical Reaction Engineering[END_REF]. For mass-transfer-limited reactions, i.e. HCl-calcite, the dissolution process is very fast and the boundary condition at the fluid-solid interface, (2.1), becomes

c Aβ = 0 at A βσ . (2.2)
This condition will be used in the rest of the paper. However, there would be no mathematical difficulty in developing a similar calculation for the general case.

Given all the above assumptions, the final pore-scale boundary value problem can be written as follows:

momentum balance

-∇p β + µ∇ 2 v β = -ρ β g, (2.3) BC 1 v β -n βσ n βσ • v β = 0 at A βσ ;
(2.4) species and total mass balance

∇ • v β = 0, (2.5) ∂c Aβ ∂t + ∇ • (v Aβ c Aβ ) = 0 in the β-phase, (2.6) ∂c σβ ∂t + ∇ • (v σβ c σβ ) = 0 in the β-phase, (2.7) ∂ρ σ ∂t = 0 in the σ-phase, (2.8) BC 1 ρ β (v Bβ -w) • n βσ = 0 at A βσ ,
(2.9) where w • n βσ is the velocity of the A βσ interface. Here, we have written the boundary condition for the passive fluid constituent (usually water), instead of a boundary condition for the mixture, since it will be useful in developing the dissolution model: BC 2 c Aβ = 0 at A βσ ;

(2.10) BC 3

-ρ σ w • n βσ = c σβ (v σβ -w) • n βσ at A βσ = -βc Aβ (v Aβ -w) • n βσ .
(2.11) Equation (2.11) represents the relation between the dissolved solid flux and the consumed acid flux, with β being the stoichiometric coefficient of the chemical reaction.

Solving the set of equations at this scale is a tremendously complicated task and would require knowledge of the pore structure; we adopt in this paper a more macroscopic point of view, using Darcy-scale equations.

Darcy-scale equations

It is beyond the scope of this paper to present a complete derivation of the Darcyscale equations from the above pore-scale problem. Equations have been introduced heuristically in chemical engineering (see for instance a general introduction about chemical reaction engineering in [START_REF] Fogler | Elements of Chemical Reaction Engineering[END_REF]. Attempts have been made to obtain the Darcy-scale equations using averaging techniques, and a discussion of the associated approximations can be found in [START_REF] Quintard | Dissolution of an immobile phase during flow in porous media[END_REF]. However, it is important to understand the limitations and problems associated with any Darcy-scale dissolution model. On this basis and in the light of a quite complete analysis of the different approaches used to describe the fluid-solid reactions in porous media provided by [START_REF] Sahimi | Statistical and continuum models of fluid-solid reactions in porous media[END_REF], we shall recall below the main aspects of the difficulties in deriving such a macroscopic model.

We will follow the notation associated with the volume-averaging theory presented in [START_REF] Quintard | Dissolution of an immobile phase during flow in porous media[END_REF]. Darcy-scale averaged quantities are defined according to [START_REF] Anderson | A fluid mechanical description of fluidized beds[END_REF], [START_REF] Marle | Ecoulements monophasiques en milieu poreux[END_REF], [START_REF] Slattery | Flow of viscoelastic fluids through porous media[END_REF] or [START_REF] Whitaker | Diffusion and dispersion in porous media[END_REF]. For instance the superficial velocity is given by

V β = v β = 1 V V β-phase v β dV , (2.12)
where V β-phase represents the volume of the β-phase contained within the averaging volume, V , represented in figure 1, and the intrinsic average pressure and concentration are defined by

P β = p β β = 1 V β-phase V β-phase p β dV , (2.13) C Aβ = c Aβ β = 1 V β-phase V β-phase c Aβ dV .
(2.14)

Liquid flow and momentum equations

At the Darcy scale shown in figure 1, three regions can be identified: a fluid zone created by complete dissolution, the initial porous medium and a transient zone characterized by a porosity gradient which may eventually develop under local non-equilibrium conditions. Fluid and porous zones will be modelled using either the Stokes equations or Darcy's law. This approach is only valid if we assume that the interface velocity is small, which is currently acceptable for most dissolution problems. At the dissolution front, a macroscopic matching must be satisfied between these two models through a macroscopic interface between the fluid zone and the porous medium in the case of local equilibrium, or through a zone with a finite porosity gradient in the case of local non-equilibrium.

The treatment of a fluid-porous medium interface has received much attention [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieu poreux. Application à la convection naturelle[END_REF][START_REF] Vafai | Analysis of flow and heat transfer at the interface region of a porous medium[END_REF][START_REF] Vignes-Adler | Transport processes along fractals. the Cantor-Taylor brush[END_REF][START_REF] Larrea | Etude par simulation numérique de l'hydrodynamique au voisinage d'une interface poreuse[END_REF][START_REF] Larrea | Sur les conditions aux limites à la paroi pour le calcul d'un écoulement laminaire dans un canal à parois poreuses[END_REF][START_REF] Naimi | Sur la formulation et simulation de l'effet d'inertie dans l'hydrodynamique en milieu poreux[END_REF]. The first approach consists in writing appropriate boundary conditions at the interface and the most popular form corresponds to the work of [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF]. For a one-dimensional flow in the x-direction of a fluid above a porous medium, this condition is written in the following form:

∂V Stokes x ∂y = α Beaver K 1/2 (V Stokes x -V Darcy x ), (2.15)
where α Beaver is a dimensionless parameter which depends on the porous medium, V Stokes is the fluid velocity, V Darcy is the superficial velocity in the porous medium, K is the permeability and x and y are coordinates along the fluid-porous medium interface and its normal respectively. The introduction of such boundary conditions has also been discussed in Ochoa-Tapia & Whitaker (1995a, b). Upscaling considerations (see a recent discussion in [START_REF] Goyeau | Averaged momentum equation for flow through a non homogeneous porous structure[END_REF] show that the Stokes equations, Darcy's law and boundary condition (2.15) can be approximated, for porous media with moderate or low permeability, by a continuous formulation corresponding to the so-called Darcy-Brinkman equation [START_REF] Brinkman | A calculation of the viscous force exerted by a flowing fluid on dense swarm of particles[END_REF]) usually written as

µ ε β ∆V β -∇P β + ρ β g -µK -1 • V β = 0, (2.16)
which can be expressed, in dimensionless form, as

1 ε β ∆V -∇P + Re Fr g |g| - 1 N D V = 0, (2.17)
where the dimensionless variables are defined as follows:

V = V β |v 0 | , P = P β l µ|v 0 | , Re = ρ β |v 0 |l µ , Fr = |v 0 | 2 |g|l , N D = K l 2 .
(2.18)

Here, N D represents the Darcy number which can be used to compare the relative magnitude of the Darcy term to the Stokes viscous term, whereas the ratio between the Reynolds number Re and the Froude number Fr defines the relative influence of gravity and viscosity effects. In (2.18) l is a pore-scale characteristic length and v 0 is the inlet velocity. The error made in using this continuous formulation is of the order of the grain size [START_REF] Vignes-Adler | Transport processes along fractals. the Cantor-Taylor brush[END_REF].

The model was originally developed as a correction to Darcy's law in the case of large porosities [START_REF] Brinkman | A calculation of the viscous force exerted by a flowing fluid on dense swarm of particles[END_REF]) and was further justified on a theoretical basis by [START_REF] Tam | The drag on a cloud of spherical particles in low Reynolds number flows[END_REF], [START_REF] Levy | Loi de Darcy ou loi de Brinkman[END_REF] and [START_REF] Rubinstein | Effective equations for flow in random porous media with a large number of scales[END_REF]. The coefficient µ/ε β in (2.16), often referred to as the effective viscosity, μ, has been the subject of debate in the literature, essentially because direct experimental measurements are very difficult. For instance, Brinkman, suggested using μ = µ whereas more complex forms involving ε β were proposed by [START_REF] Lundgren | Slow flow through stationary random beds and suspension of spheres[END_REF], [START_REF] Durlofsky | Analysis of the Brinkman equation as a model for flow in porous media[END_REF] and [START_REF] Lhuillier | La viscosité effective dans le modèle de Darcy-Brinkman: Approche self-consistante. 4èmes Journées d'Etude sur les Milieux Poreux de Nancy[END_REF]. Although this is still a matter of controversy [START_REF] Nield | The boundary correction for the Rayleigh-Darcy problem: Limitations of the Brinkman equation[END_REF], the use of such a continuous approach is also appealing for describing the flow in a fluidporous system like the one under consideration in the present work [START_REF] Beckermann | Natural convection flow and heat transfer between a fluid layer and a porous layer inside a rectangular enclosure[END_REF][START_REF] Beckermann | Natural convection in vertical enclosures containing simultaneously fluid and porous layers[END_REF][START_REF] Gartling | A finite element formulation for incompressible, conjugate, viscous/porous flow problems[END_REF]). In the case of local equilibrium, it avoids the use of special conditions like the one proposed by Beavers & Joseph required to match the Stokes model and Darcy's law. Moreover, the diffusive viscous term can properly describe the boundary layer type of flow that develops in the porous zone close to the interface and this has been validated by Vafai and co-workers [START_REF] Vafai | Boundary and intertia effects on convective mass transfer in porous media[END_REF][START_REF] Vafai | Analysis of flow and heat transfer at the interface region of a porous medium[END_REF] and by [START_REF] Lhuillier | La viscosité effective dans le modèle de Darcy-Brinkman: Approche self-consistante. 4èmes Journées d'Etude sur les Milieux Poreux de Nancy[END_REF]. In the fluid zone where N D is very large (theoretically infinite), the Stokes model is recovered while in the porous zone for which N D is small, the viscous diffusive term becomes negligible and the model reduces to Darcy's law.

In the case of local non-equilibrium where a zone of finite length with porosity gradients appears, more complicated questions arise if the continuous Brinkman model is used. In this case, it can be shown theoretically that additional terms associated with porosity gradients are required when upscaling the momentum porescale equation (Ochoa-Tapia & Whitaker 1995a, b;[START_REF] Goyeau | Averaged momentum equation for flow through a non homogeneous porous structure[END_REF][START_REF] Whitaker | The Method of Volume Averaging[END_REF]. Although further work is required for more precision about the contribution of such terms, studies on real structures with varying porosities indicate that this model is very precise when used with a varying permeability provided these porosity gradients are not too large [START_REF] Goyeau | Numerical calculation of the permeability in a dendritic mushy zone[END_REF]. As a first approach, this encourages the use of the Brinkman equation (2.16) for the momentum balance in the different situations under consideration in this work. Moreover, the transition zone, when present, between the Stokes and Darcy types of flow, is probably of very small extend. Results will be used a posteriori to check the validity of this hypothesis. Now we will study the relationship between permeability and porosity. Permeability is only well defined when the geometry of the fluid-solid interface is known. However, this interface may evolve during the dissolution process in a complex manner leading to an evolution of the phase indicator and hence of the porosity. A direct relationship between these two macroscopic quantities may depend on the history of the dissolution, due to the coupling between flow and dissolution processes. This has been investigated on the basis of pore-scale simulations by [START_REF] Bekri | Dissolution of porous media[END_REF]. Similar problems arise with dendritic mushy zones [START_REF] Poirier | Permeability for flow of interdentritic liquid in columnar-dentritic alloys[END_REF][START_REF] Goyeau | Numerical calculation of the permeability in a dendritic mushy zone[END_REF]. For a given class of dissolution problem and a given sample, it is reasonable to assume a direct relationship between the permeability and the porosity, i.e. K(ε) instead of K(t) and ε(t) [START_REF] Fogler | Predicting the flow and reaction of HCl/HF acid mixtures in porous sandstone cores[END_REF][START_REF] Lhuillier | La viscosité effective dans le modèle de Darcy-Brinkman: Approche self-consistante. 4èmes Journées d'Etude sur les Milieux Poreux de Nancy[END_REF]. Here, we will adopt this point of view, remembering that it may be the cause of discrepancies between observed results and predictions. It should be emphasized however that the dissolution model proposed here does not depend on a particular choice for the permeability-porosity relationship. Any choice could be made, depending on the amount of information available for the system under consideration. In the section concerning the comparison between experiments and theoretical predictions, we will use the classical Kozeny-Carman relation [START_REF] Ergun | Fluid flow through packed columns[END_REF]) to represent such a relationship, mainly because no information on the real pore-scale structure is available. If we consider that the tortuosity of the medium evolves linearly according to the particle specific area during the dissolution, we can write

K l 2 β = k cte ε β 1 -ε β 3 , (2.19)
where l β is the characteristic length at the pore-scale (mean size of the grains) and k cte is a constant depending on the nature of the medium. Moreover, although the use of the Kozeny-Carman relationship is a strong approximation, it plays a minor role in local equilibrium dissolution where a sharp front delimits fluid and porous zones (no transition zone for porosity or permeability). Its influence is sensitive only for the ramified or uniform regime (local non-equilibrium dissolution) where a small change in porosity may result in a change in permeability of several order of magnitude during wormhole formation.

It will be shown in Appendix A that the overall mass balance equation for the fluid phase may be approximated by ∇ • V β = 0, (2.20) and this completes the formulation for the flow problem.

Dissolution and mass balance equations

From the literature review either for the tube model or the network, we have seen that the mean concentration in the porous medium is not necessarily zero and for this reason a model of local mass non-equilibrium is adopted. Doing so yields a rather general model able to reproduce any kind of dissolution regime from compact (sharp fluid-porous medium interface), to uniform (small porosity gradient). It must be recalled here that the assumption of a mass-transfer-limited reaction, which leads to a boundary condition of zero acid concentration at the interface between the liquid and the solid at the pore scale, should not be confused with a local mass equilibrium dissolution, which corresponds to a condition of zero acid concentration at the fluidporous medium interface at the Darcy scale. Both are distinct mechanisms at different scales, and a mass-transfer-limited reaction may lead to either a local equilibrium or a local non equilibrium process. This is a key distinction that must be emphasized since it has often been a source of confusion in the literature, especially in capillary tube models. Non-equilibrium models are frequently used for the problem of heat transfer (non-local thermal equilibrium). Results by [START_REF] Carbonell | Fundamentals of Transport Phenomena in Porous Media[END_REF], Zanotti & Carbonell (1984a, b), [START_REF] Kaviany | Principles of Heat Transfer in Porous Media[END_REF] and [START_REF] Quintard | Two-medium treatment of heat transfer in porous media: Numerical results for effective properties[END_REF] have shown that this assumption requires the introduction of two macroscopic equations. With regard to the dissolution, this type of assumption has been used for the modelling of pollutant transport [START_REF] Dezabala | A non-equilibrium description of alkaline water flooding[END_REF][START_REF] Powers | Theoretical study of the significance of nonequilibrium dissolution of nonaqueous phase liquids in subsurface systems[END_REF][START_REF] Miller | Multiphase flow and transport modeling in heterogeneous porous media: Challenges and approaches[END_REF], or for sandstone acidizing, for example (Fogler et al. 1967;[START_REF] Schechter | Oil Well Stimulation[END_REF].

In this paper, we do not discuss thoroughly the comprehensive development of the dissolution Darcy-scale model. Several points have already been discussed in the literature. However, since we believe it is important to understand the limitations of the Darcy-scale approach, the reader will find in Appendix A a summary of the different steps and assumptions made in order to obtain the macroscopic equations outlined below. When the change in the interface shape is slow enough, a good approximation of the problem leads to the following form of the local equations:

ε β ∂C Aβ ∂t + V β • ∇C Aβ = ∇ • (D * β • ∇C Aβ ) -αC Aβ , (2.21) ∂ε β ∂t = βαC Aβ ρ σ , (2.22)
in which some additional terms have been discarded because of their negligible effect.

The key parameter is the mass exchange coefficient α, and it should be noted that if its value is large enough, the resulting averaged concentration will be close to the equilibrium concentration, i.e. C Aβ = 0. Therefore, this model can be used to simulate both local equilibrium and local non-equilibrium conditions.

In the theoretical development presented in [START_REF] Quintard | Dissolution of an immobile phase during flow in porous media[END_REF], two porescale 'closure problems' are described relating the pore-scale physical characteristics to the 'effective' properties in (2.21). These closure problems can be used to estimate the β-phase dispersion tensor, D * β , and the mass transfer coefficient, α, for a given shape of the fluid-solid interface. As outlined by those authors, these closure problems may be used in two different ways:

(i) To compute the effective transport parameters for a series of given interface geometry. For instance, one may assume homogeneous dissolution of the interface.

(ii) To compute the effective transport parameters, starting with a given interface, and then, for specified Darcy-scale concentration fields, these closure problems may be used to calculate the dissolution rate at each point of the interface. This will allow us to compute the new interface shape, and this process could be repeated iteratively. In that way, pore-scale closure problems, and Darcy-scale equations would be solved in a coupled manner, and this would account for very complex dissolution processes. This second approach is of course much more sophisticated than the first, and it is of some importance in terms of computational time to see under which conditions it is necessary to maintain this complexity.

For permeability, effective properties may depend on the dissolution process in an historical manner (see for instance the pore-scale simulations developed by [START_REF] Bekri | Dissolution of porous media[END_REF]. In this work, we adopt the traditional approach, and we approximate the historical behaviour of the system by a direct relationship with the porosity of the system, and a pore-scale Péclet number: D * β (ε β , P e cell ) and α(ε β , P e cell ). However, similarly to the discussion on the permeability-porosity relationship, the choice of a particular correlation is of secondary importance in the model. Any choice can be made without causing difficulties for the model implementation. Values of the mass transfer coefficient reported in the literature were determined experimentally either on annular flow reactors [START_REF] Mumallah | Factors influencing the reaction rate of hydrochloric acid and carbonate rock[END_REF][START_REF] Mumallah | Effective HCl diffusion coefficients from correlations of HCl-limestone reactions[END_REF][START_REF] Mumallah | Reaction rates of hydrochloric acid with chalks[END_REF] or rotating disk devices [START_REF] Lund | The dissolution of calcite in hydrochloric acid[END_REF]. Nevertheless, since the value of the α-coefficient strongly depends on the geometry of the domain, it is difficult to extrapolate from these works an effective value for the dissolution of a real porous medium. In the work of Fredd & Fogler (1998b), the mass transfer coefficient was extracted from the analytical solution of the convection-diffusion equation in the laminar regime in a capillary tube as derived by [START_REF] Levich | Physiocochemicals Hydrodynamics[END_REF]:

K mt = 1.86D 2/3 u w d w L w 1/3 , (2.23)
where K mt is the Levich's mass transfer coefficient, d w and L w are the diameter and length of the capillary tube, respectively, and u w is the superficial velocity in the capillary tube. This correlation has the advantage of being dependent on the diffusion coefficient and velocity. However, it is important to recall that it is obtained by neglecting the second-order derivative of the concentration along the tube axis, implying that the diffusive part is not taken into account in the expression for the mass transfer coefficient. As a consequence, this form remains valid for large enough Péclet number and its use in the conical regime where the Péclet number is small may underestimate the diffusion effects. Moreover, there is confusion between the different scales when using the wormhole properties instead of pore-scale properties in order to define the mass transfer coefficient. As expressed by (2.23), K mt is a core-scale mass transfer coefficient and not a Darcy-scale mass transfer coefficient. In the Darcy-scale approach, such Darcy-scale non-equilibrium processes are modelled directly by the interplay of the different diffusion and convection mechanisms.

In the absence of experimental data for salt dissolution, the mass transfer and dispersion coefficient used in our simulations were obtained first from a numerical solution of the closure problems proposed in [START_REF] Quintard | Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF], 1999). The unit cell geometry used in our computations (see Appendix B) can be as complex as necessary in order to capture the features of a real system and can account for the complexity of the velocity field due to the non-local effects in the porous medium. Furthermore, results presented in the next section show that the qualitative behaviour of the model, in terms of dissolution regimes, does not significantly depend on a particular choice of the unit cell. Because of the lack of information on the porous medium used in the experiments, the unit cell is relatively simple. However, for the comparison with the experimental data and the prediction of the optimum flow rate, it was necessary to use an estimation of the effective coefficients closer to the real structure, and this is where the above-mentioned upscaling techniques may provide some quantitative information. The reader can find a presentation of these 'closure problems' and the outline of the calculations in Appendix B.

At this point, we have a complete system of macroscopic equations able to describe acid transport and dissolution in the domain at the Darcy scale. These equations can be made dimensionless by introducing the classical dimensionless variables as follows:

D = D * β D , V = V β |v 0 | , x i = x i l , t = |v 0 |t l , C = C Aβ c 0 , (2.24)
where |v 0 | and c 0 correspond to inlet flow and concentration conditions, and l is a characteristic length. This leads to

ε β ∂C ∂t + V • ∇C = 1 P e ∇ • (D • ∇C ) -Da C , (2.25) ∂ε β ∂t = (1 -ε β ) ε β Da N ac C . (2.26)
It is clear from these two equations that, in addition to Re, Fr and N D already mentioned, the wormhole propagation is a function of the following dimensionless numbers classically reported in the literature [START_REF] Buijse | Mechanisms of wormholing in carbonate acidizing[END_REF]Daccord et al. 1993a, b):

the Péclet number which represents the magnitude of the convective transport versus the diffusional transport,

P e = |v 0 |l D ; (2.27)
the Damköhler number defined as the ratio of acid consumed and the acid transported by convection,

Da = αl |v 0 | ; (2.28)
the kinetic number defined as the ratio between the acid flux consumed by the chemical reaction and the diffusive flux,

K i = Da P e = αl 2 D ; (2.29)
the acid capacity number which represents the dissolving power of the acid solution, that is to say the mass of solid dissolved per unit mass of the rock,

N ac = ε β c 0 β (1 -ε β )ρ σ .
(2.30)

These dimensionless numbers include a characteristic length, l, which has been chosen arbitrarily equal to √ K, a representative length at the pore scale. In the next section, the numerical model developed to solve the complete system of equations is presented.

Numerical model

The equations are discretized on a three-dimensional uniform Cartesian grid. The geometry of the medium under consideration is approximated by taking constant values for permeability and porosity on every mesh of the domain at each time step. Various boundary conditions have been implemented: symmetry, impermeable wall, periodic conditions and, for the inlet, injection at constant pressure or constant flow rate.

Momentum equations

The numerical model uses a multiple-step approach to account for the various mechanisms. Due to the tracer assumption, the stationary Darcy-Brinkman problem is solved first for a given porosity field using a Predictor-Corrector method coupled with a Uzawa algorithm [START_REF] Robichaud | An iterative implementation of the Uzawa algorithm for 3-D fluid flow problems[END_REF]. The difficulty of this linear system is the coupling between a first-order and a second-order operator, which prevents the use of a classical resolution method for the Stokes equations, such as the Augmented Lagrangian method. We decompose the pressure and velocity fields as follows:

P β = P * β + Pβ , (3.1) V β = V * β + Ṽ β , (3.2)
where P * β and Pβ represent respectively the predicted and corrected pressure field, and V * β and Ṽ β the predicted and corrected velocities respectively.

Predictor step

We first calculate the predicted pressure and velocity fields on the basis that they will lead to a good approximation of the Darcy part in the equations. Therefore, P * β and V * β are sought as solutions of the following classical Darcy problem:

∇P * β = ρ β g -µK -1 • V * β , (3.3) ∇ • V * β = 0.
(3.4) Equations (3.3) and (3.4) are combined to yield

∇ • K µ • (∇P * β -ρ β g) = 0. (3.5)
Classical schemes to solve one-phase Darcy flow in heterogeneous media can be applied. Here, a finite volume formulation is used over the Cartesian grid and cellto-cell fluxes are approximated by using a harmonic averaged transmissivity ensuring an order-2 scheme [START_REF] Ahmadi | Large-scale properties for two-phase flow in random porous media[END_REF]. This results in a linear system that is solved using an iterative algorithm like bi-cgstab which is found to be 15% to 20% more effective than the classical conjugated bigradient (Van Der Vorst 1992). The V * β -field is computed from the pressure field P * β by making use of (3.3).

Corrector step

After the predictor step, we solve the equations for the corrected pressure and velocity, which are written as

µ ε β ∆ Ṽ β -∇ Pβ -µK -1 • Ṽ β = - µ ε β ∆V * β , (3.6) ∇ • Ṽ β = 0.
(3.7) We solve this system with a Uzawa iterative algorithm [START_REF] Fortin | Méthodes de Lagrangien Augmenté[END_REF] which can be summarized as follows Given P n β , find P n+1

β and Ṽ n β solutions of µ ε β ∆ Ṽ n β -µK -1 • Ṽ n β = ∇ P n β - µ ε β ∆V * β , (3.8) P n+1 β = P n β -ρ∇ • Ṽ n β ,
(3.9) where ρ is the descent step parameter (with 0 < ρ < 2).

The convergence of the scheme is imposed by the following condition:

| P n+1 β -P n β | 6 ε 0 , (3.10)
corresponding to a velocity field with divergence lower than ε 0 /ρ. A staggered mesh for pressure and velocity is used in order to avoid pressure oscillations with respect to the space coordinates.

Mass transport equations

Effective coefficients in the transport equations are computed next for a given porosity and velocity field. A time-splitting method is used which splits the equation into a hyperbolic part (convective terms) and an elliptic part (diffusive term). A scheme proposed by [START_REF] Takacs | A two-step scheme for the advection equation with minimised dissipation and dispersion errors[END_REF] is used for the hyperbolic part in order to reduce numerical diffusion [START_REF] Bruneau | An accurate finite difference scheme for solving convection-dominated diffusion equations[END_REF]. It is a TVD (total variation diminishing), oscillation-free, and L ∞ -stable scheme. The diffusive part is solved using a classical implicit discretization. The scheme is written in three dimensions by direction. Numerical tests were performed to check the sensitivity of the scheme to the grid orientation, which is found to be very weak. The resulting acid concentration profile is then used to solve the dissolution equation, (2.22), providing a new porosity field. From this point on, the overall procedure is repeated iteratively to update the dissolution pattern.

Instability formation

In order to obtain instability formation, it is necessary to introduce some disturbances in the system. Several possibilities are available: (i) perturbation of the initial concentration field, (ii) perturbed boundary condition, (iii) perturbed porous medium characteristics. We chosed the third by adding to the homogeneous permeability field a random noise using a uniform distribution law of zero correlation length. In a real medium, these heterogeneities always appear at all scales. Small fluctuations, ε, added to the permeability field can be seen as a consequence of the pore-scale heterogeneities. The influence of the initial perturbation on the dissolution pattern and dynamics has been tested by varying the amplitude. For small amplitudes (around 10%), a very small impact on the wormhole development has been observed. This is the range of pertubations used in all our quantitatives studies. However, we observe few lateral ramifications in this case, the wormhole shape being more schematic. In order to obtain more detailed dissolution pattterns, it is necessary to increase the amplitude of the permeability perturbations (about a factor 5). This type of perturbation was only used for the presentation of the dissolution patterns.

Computations have been performed for many different conditions on twodimensional domains, and results are now discussed in the next section and compared to the experimental data.

Results and discussion

In this section, we first present the experimental data that will be used in the subsequent discussion. We explain why the studied system is of interest and the results obtained.

In most experiments on dissolution available in the literature, the dissolution pattern is three-dimensional and detailed evolution of the dissolution pattern with time is not available. For direct comparison with two-dimensional numerical simulations presented in this work, detailed experimental results are needed and for this reason, simple dissolution experiments were performed on two-dimensional beds of packed grains of salt.

4.1. Experimental study 4.1.1. Choice of the system

The experimental results are obtained by the injection of an under-saturated salt solution in a porous medium made of salt grains. This system is comparable to limestone dissolved by acid since the dissolution is limited by the mass-transfer mechanism. The main advantage is that it is possible to take pictures as the injection proceeds. In addition, experiments are fast and easy to perform. Comparison between the acid-limestone and the water-salt pack systems is easy to perform on the basis of the dissolution power quantified with N ac . The fraction of rock dissolved by the acid, βC Aβ , corresponds to the difference between the solubility limit of NaCl in water (C sat = 360 g l -1 at 20 • C) and the salt concentration C NaCl . We have

C Aβ = C sat -C NaCl β . (4.1)
Since the model is written for acid injection, this relation is used for the modelling of the experimental results where an equivalent acid concentration is associated with every salt concentration. Note that in this system an increase in salt concentration is equivalent to a decrease of the acid concentration in the acid-limestone system.

Experimental procedure

The experimental setup is represented in figure 2. Salt is packed in a transparent Hele-Shaw cell of 25 cm in length, 5 cm wide, and 1 mm in depth which can therefore be considered as two-dimensional for comparison with numerical simulations. Salt grains are about 300 µm in mean size, having the form of small plates. Once packed, the porosity is about 0.36. At the beginning of the experiment, the medium is completely saturated with salt water at 360 g l -1 and permeability is measured. The average value of the permeability is 1.5 × 10 -11 m 2 . A solution of given concentration is injected at constant flow rate at ambient temperature while the pressure drop between the inlet and outlet is measured with a pressure transducer. The evolution of the dissolution pattern is recorded using a video camera. Both camera and pressure transducer are computer-monitored. Experiments have been performed for salt concentrations between 50 and 340 g l -1 , and for injection rates between 5 and 350 cm h -1 . Reproducibility of the experiments has been checked and is very good as illustrated by the results in figure 3 where we have reported the wormhole length versus time for two experiments performed at the same conditions.

Dissolution patterns

Qualitative results (figure 4) are identical to those of the classical acid-calcite system, as reported in the literature [START_REF] Hoefner | Pore evolution and channel formation during flow and reaction in porous media[END_REF]). According to [START_REF] Fredd | Validation of carbonate matrix stimulation models[END_REF] (e) Uniform dissolution: Finally, at very high injection rates, the reactant is forced into all pores in the porous matrix. The dissolution front is spread over the whole core length and this corresponds to a non-equilibrium dissolution. It is remarkable that a local non-equilibrium leads to another stable displacement. This latter regime could not be reached in our experimental setup due to the use of an unconsolidated porous material. It was impossible to perform injection at very high velocities to investigate the uniform dissolution. Finally, an example of the wormhole growth is presented in figure 5, where L w is the wormhole length.

Optimum flow rate

The experiments confirmed the existence of an optimum injection rate, i.e. minimum quantity of injected solution to break through the core. The experimental results obtained for salt concentrations equal to 150 and 230 g l -1 are presented, in terms of the number of pore volumes required to extend the dissolution channels to a specified length, here the core length L, in figures 6(a) and 6(b). This volume is usually referred to as the pore volume to breakthrough, V BT . Clearly, these figures show the existence of an optimum injection rate and are characterized by two different behaviours. Breakthrough times decrease as the injection rate increases until the optimum flow rate is reached. Then the breakthrough times increase again for injection flow rates higher than the optimum.

In order to emphasize the existence of an optimum flow rate, results from figure 6 were combined, along with results obtained with C NaCl = 50 g l -1 , in figure 7 where we have represented the pore volume to breakthrough as a function of the injection rate. The optimum flow rate is reported in log-log coordinates. We obtained experimentally optimum flow rates of 100, 80 and 70 cm 3 h -1 for the salt solutions at 50, 150 and 230 g l -1 respectively. A slope of 1/3 for injection rates higher than the optimum is observed, which is in perfect agreement with results obtained for HCl-calcite systems [START_REF] Daccord | Carbonate acidizing: Toward a quantitative study of the wormholing phenomenon[END_REF][START_REF] Bazin | From matrix acidizing to acid fracturing: A laboratory evaluation of acid/rock interactions[END_REF].

Considering the effect of the salt concentration on the optimum flow rate, we observe an increase of the optimum flow rate when the salt concentration decreases. This result is consistent with results obtained with the limestone-acid system since decrease of the salt concentration in the salt-undersaturated salt solution system corresponds to an increase of the acid concentration in the acid-limestone system [START_REF] Wang | The optimum injection rate for matrix acidizing of carbonate formations[END_REF][START_REF] Bazin | Acid filtration in dynamic conditions to mimic fluid loss in acid fracturing: SPE 38168[END_REF].

Finally, breakthrough times decrease with NaCl concentration, following the trend observed with the acid-limestone system.

As mentioned above, we have observed trends in the experimental results in agreement with results obtained by previous investigators. They can be interpreted a priori with a two-dimensional model. In the next section, results obtained from numerical simulations are compared to these experimental data. Given the equivalence between salt and acid dissolution models, we will keep the terminology of acid injection for both phenomena.

Numerical simulations and comparison with experimental data

A large number of simulations have been performed, at various flow rates and other parameter values. They are presented below.

Dissolution regimes

It is first of interest to check whether the model is able to reproduce the dissolution regimes observed experimentally. The two-dimensional configuration under study is illustrated figure 8 and conditions for the simulation are given in table 1. We set the mass transfer coefficient and the diffusion coefficient at constant values to eliminate Da and P e variations during the simulation. Our numerical model works as well with non-constant α, which will be used in the rest of the paper. However, keeping α constant provides a good illustration of this Darcy-scale model, in terms of dissolution instabilities, and this opens a large field of investigations for further studies.

Figure 9 shows the different dissolution patterns obtained numerically for different acid injection rates. The concentration fields are represented in figure 10. Figure 11 shows the porosity field during the wormhole growth in the wormholing regime. The obtained dissolution patterns are remarkably similar to the experimental dissolution patterns shown figure 4. They are used to clarify the dissolution mechanisms. The concentration fields show that regimes (a-c) correspond to local equilibrium conditions, while regimes (d -e) feature local non-equilibrium dissolution.

At this point, we see that the Darcy-Brinkman non-equilibrium model proposed in this paper has all the required features to capture the complex observed physical behaviour. It is shown that all the dissolution structures can be described with a mass transfer coefficient, i.e. without considering a reaction-limited kinetics, and a flow description which allows the Stokes and Darcy equations. We observed that the Stokes flow is prominent in the wormhole zone when the porosity is close to 1, while Darcy's law becomes dominant in the surrounding porous medium. We can conclude that the Darcy-Brinkman model is well suited to describe fluid loss through the wormhole. Considering mass transfer, local-equilibrium and nonequilibrium situations are described as a result of the interplay between fluid flow and the mass transfer coefficient.

Dissolution diagrams

We take advantage of the physics introduced in the numerical model to map the dissolution structures in behaviour diagrams (figures 12 and 13). The P e-Da and P e-Ki diagrams have been drawn from simulations where some different and independent constant values for the effective coefficients α and D have been taken for each run. This offers an interesting general perspective on the theory of such unstable displacements, which we will use here. Of course, it might be difficult in practice to keep all the parameters independent, and this will be the case when the nonlinear relations for α and D * β are used, as presented in Appendix B. This will be illustrated, for instance, in the discussion about figure 17. Boundaries have been drawn between the different regimes based on a qualitative observation of the dissolution structures and some uncertainties are associated with the transition lines. Moreover, we have fixed N ac and the value of the core length, L. It must be emphasized that the modification of these values leads to a displacement of the transition lines, as it will be shown in the last part of the paper. Although the choice of the α-correlation will have an impact on the boundary position, the use of a constant value for the mass transfer coefficient gives an accurate enough idea of transitions between the regimes.

The P e-Da diagram (figure 12) shows that the transition between the compact, conical wormhole and dominant wormhole regimes does not depend on the Damköhler value. It is only determined by the value of the Péclet number. Decreasing the Damköhler number will lead to more local non-equilibrium effects and consequently to transitions towards uniform dissolution. This is also illustrated in the P e-Ki diagram, figure 13, which shows that it is possible to go from the compact to the uniform regime by increasing only the injection flow rate.

The P e-Da diagram (figure 12) provides other information concerning the wormhole formation: the instabilities, i.e. wormholes, start developing above threshold values of the Damköhler and Péclet numbers. For P e 6 2 × 10 -3 , the diffusion mechanism is dominant and the dissolution patterns are in compact or uniform regimes. In the same way, for Da 6 8.10 -4 , the dissolution front is completely spread over the whole width of the core, the instabilities are eliminated and the uniform regime is observed.

Finally, we compare our P e-Ki diagram to the diagram previously proposed by Daccord et al. (1993b), figure 14, for the simple geometry of a straight capillary. They consider a capillary of radius R, in a soluble mineral, in which a reactive fluid of concentration C is injected at a constant rate, q. The dimensionless numbers are expressed as

P e Daccord = q πRD , ( 4.2)
and

Ki Daccord = k r C aD , (4.3)
where k r is the reaction rate constant and a the specific molar area of the mineral constitutive of the capillary. The comparison of figure 14 and figure 13 shows many similarities. The same slopes between the regimes of compact/conical wormholes-dominant wormholes and compact/conical wormholes-ramified wormholes/uniform are obtained. The difference in slope for the limit between the ramified wormholes/uniform regime and dominant wormholes (1/3 instead of 1 in the present work) could be explained by the fact

C acid = 545 kg m -3 ρ σ = 2160 kg m -3 µ = 1 × 10 -3 Pa s ε β = 0.2 β = 1.62 D β = 2 × 10 -9 m 2 s -1 K = 10 -11 ± ε m 2
Table 2. Numerical data for constant flow rate injection.

that Daccord et al. (1993b) do not distinguish the uniform regime from the ramified wormhole regime.

Optimum flow rate

Numerical experiments were performed to calculate the wormhole growth curves for various injection flow rates. In these simulations the mass transfer coefficient is varied according to the local porosity and velocity. The expression for the mass transfer coefficient is α = Aα o (P e cell , ε β ), (4.4) where α o is given by figure 25 (Appendix B). The A coefficient is introduced arbitrarly as a shape factor to account for various pore structures and to test the effect of α on the breakthrough times and optimum flow rates. Values of A are 0.1, 0.2 and 0.5. The numerical data used for the simulations are shown in table 2. The core is 25 cm in length and 5 cm in width. Note that the numerical estimation of the wormhole length depends on the choice of the porosity threshold used to determine the tip of the channel. Arbitrarily, we have considered that the wormhole ends when porosity reaches 0.75, but a different choice would slightly modify the value of the pore volume to breakthrough. For instance, time to breakthrough may be about 5% longer if we had chosen 1 as the final porosity value. Resulting data are plotted in figures 15(a), 15(b) and 15(c). Results showed qualitatively a wormhole growth comparable to the experimental results and the existence of an optimum flow rate. The optimum flow rates are 10, 20 and 50 cm 3 h -1 for A values of 0.1, 0.2 and 0.5 respectively. The optimum flow rate increases linearly with α. As expected, the mass transfer coefficient is the key parameter of the model, and a small variation in the estimation of α leads to an error of the same order of magnitude in the prediction of the optimum flow rate. Therefore, this coefficient must be determined carefully.

Figure 16 shows that the α value has a substantial impact on wormhole growth at flow rates lower or higher than the optimum. For flow rates higher than the optimum, the wormhole growth increases with the value of the mass transfer coefficient whereas it is the opposite for the flow rates lower than the optimum. This result may be explained physically by considering the dissolution mechanisms. As demonstrated with the dissolution patterns, the dissolution occurs at the tip of the wormhole for high flow rates whereas for small injection velocities, a larger quantity of reactant is consumed at the walls of the wormholes, thus decreasing the wormhole tip propagation.

Optimum Damköhler number

The results of the numerical experiments performed for different values of the mass transfer coefficient are reported in a P e-Da diagram (figure 17). It is found that optimum flow rates actually correspond to a unique constant Damköhler number, Da opt , in agreement with previous work (Fredd & Fogler 1998b). This is evident as we claim in the preceding paragraph that the optimum flow rate varies linearly with α. This optimum Damköhler number is located near the transition zone between a dominant wormhole and ramified wormholes, an idea also put forth by [START_REF] Wang | The optimum injection rate for matrix acidizing of carbonate formations[END_REF]. Once again, it should be recalled that all the results are obtained for a given core length and acid concentration, as expressed through N ac . Fredd & Fogler (1998b) obtain Da Fredd-opt = 0.29, whatever the system, with a definition of the Damköhler number involving the characteristic length of the domain. However, Fredd & Fogler's characterization of the optimum Damköhler number has been verified for several experiments performed only at a constant acid capacity number N ac and a given core length (concerning the influence of the acid capacity number, they claim that it slightly modifies the 0.29 value obtained for the Da Fredd-opt ). They use the following definition:

Da Fredd = 4L w k tot d w u w , (4.5)
where u w is the interstitial velocity in the wormhole, d w and L w are the diameter and length of the wormhole respectively (their evaluation is based on final wormhole dimensions), and k tot is the dissolution rate constant (k tot = K mt for a mass-transferlimited reaction). With the formalism used in this work, K mt corresponds to αl, l being a microscopic characteristic length. In the case of a mass-transfer-limited reaction, where q w is the flow rate in the wormhole. According to this definition, the optimum flow rate must increase linearly with the core length and this is not consistent with other experimental results [START_REF] Bazin | From matrix acidizing to acid fracturing: A laboratory evaluation of acid/rock interactions[END_REF]. This apparent excellent agreement between the two values of the optimum Damköhler number must however be put into context since this results from the choice of the characteristic length. For instance, the choice of l β (the unit cell length of figure 24), as the pore characteristic length, would give

Da opt = 4L w αl β d w u w ≈ 14.8. (4.8)
Yet the choice of √ K as the pore characteristic length seems the best one in the context of Fredd & Fogler's model. Their estimate of the mass transfer coefficient is based on a tube model and we have already discussed the limitations of Levich's solution, especially at very low Péclet number. If we leave this discussion open and apply this estimate to the conditions of our numerical experiments (in the case A = 0.5, for example), we obtain

K mt = 1.86D 2/3 u w d w L w 1/3 = 8.06 × 10 -6 m s -1 (4.9)
for the value of the Levich's mass transfer coefficient K mt . This can be compared to the value of our mass transfer coefficient, i.e.

α √ K = 6.32 × 10 -6 m s -1 .

(4.10)

The two values are very close. Among various consequences of these results, it is interesting to point out that the variation of the acid diffusion coefficient in the fluid phase does not directly modify the value of the optimum injection rate. Actually, the Damköhler number depends only on the injection velocity and the mass transfer coefficient and diffusion effects are incorporated into this latter coefficient. However, injection of weak acid or microemulsion [START_REF] Hoefner | Role of acid diffusion in matrix acidizing of carbonates[END_REF][START_REF] Bazin | Experimental investigation of some properties of emulsified acid systems for stimulation of carbonate formations[END_REF][START_REF] Conway | A comparative study of straight/gelled/emulsified hydrochloric acid diffusivity coefficient using diaphragm cell and rotating disk[END_REF][START_REF] Lynn | A core based comparison of the reaction characteristics of emulsified and in-situ gelled acids in low permeability, high temperature, gas bearing carbonates[END_REF], for which diffusion coefficients are two orders of magnitude lower than in HCl aqueous solutions, were revealed to be efficient acidification processes with optimum conditions for lower velocity injections, and this may be explained by the relation between the mass transfer coefficient, α, and D (see Appendix B). Since α decreases when D decreases, the optimum Damköhler number is reached for small values of the flow rate in the case of weak diffusion of the acid. As a consequence, cores stimulated with weak acids exhibit a more uniform dissolution. This is consistent with the fact that, for a given flow rate, decreasing the Damköhler number, leads to dissolution patterns closer to ramified wormholes. The same argument can be put forth to explain the effect of temperature variations. When temperature increases, the diffusion coefficient increases, and, therefore, the optimum flow rate increases while dissolution patterns are closer to conical wormholes. Experimental results [START_REF] Wang | The optimum injection rate for matrix acidizing of carbonate formations[END_REF][START_REF] Bazin | Acid filtration in dynamic conditions to mimic fluid loss in acid fracturing: SPE 38168[END_REF]) are in agreement with these trends. Injection flow rate (cm 3 h -1 ) Pore volume to breakthrough, V BT /V P (×10)

Effect of concentration

We shall now test the model by looking at the effect of the acid concentration on the optimum flow rate. Some simulations were performed at C acid = 326 g l -1 and A = 0.5. Note that there is a priori no reason to change the α value when the concentration changes since α is a parameter which depends only on the porosity and the velocity. Results obtained for the different concentrations are illustrated in figure 18 in terms of pore volume to breakthrough as a function of the injection rate. Contrary to the experimental results, we do not observe a significant variation of the optimum flow rate. This could be explained by the small increase of the concentration. However we obtain an increase of the pore volume to breakthrough, as it has been observed in the literature.

Finally, the effect of the concentration on the P e-Da diagram is illustrated figure 19. The results are reported in the behaviour diagram, where the symbols refer to the highest concentration. There is an appreciable shift of the boundaries between the various dissolution structures. As an example, conditions leading to conical dissolution structures lead to compact structures when the concentration in the injected fluid increases. Consequently, a shift towards a higher optimum flow rate is expected when the acid concentration is increased. This is in agreement with the results observed experimentally.

Representation of the experimental data and calibration of the model

We shall now investigate the ability of the Darcy-scale model to correctly reproduce the experimental data and, in particular, to predict the pore volume to breakthrough, i.e. the dissolution kinetic. Two parameters have a strong influence on the breakthrough time: (i) α-correlation and (ii) the permeability-porosity correlation [START_REF] Civan | Scale effect on porosity and permeability: Kinetics, model and correlation[END_REF]. Therefore the calibration of the model reduces to adjusting these two parameters to the experimental results. While the procedure to determine the mass transfer coefficient may be applied to any pore-scale structure, the current unit cell used to solve the closure problem is a crude representation of the real medium with strong Possible routes would be to use the full coupling between the Darcy-scale model and the microscale closure problems, as outlined in [START_REF] Quintard | Dissolution of an immobile phase during flow in porous media[END_REF]. However, a detailed characterization of the pore scale would be necessary. It is beyond the scope of our work to do such a task for the salt porous domains. In the absence of pore-scale simulations on unit cells more representative of the real medium (three-dimensional calculation, no periodic boundary condition, geometry variation), and given the lack of detailed information on the pore-scale geometry for our porous system, the model is calibrated through the choice of the shape factor A. Only one series of simulations is necessary to adjust the parameter α. Once α is determined, it is in principle possible to predict the optimum flow rate and dissolution pattern of the system under consideration for other injection conditions: variation of the injected acid concentration or the core length. However, we have shown in the theoretical development that there is a priori no intrinsic value of α because of historical effects, and this is also true for other effective properties. It is therefore expected that the α-correlation may require adjustments under experimental conditions far from the range of parameters studied. We do not have a comprehensive analysis of this problem for our experimental setup; however, some elements of this discussion are presented below.

The impact of the choice of the permeability-porosity correlation on the time to breakthrough has been tested with two different correlations. The numerical simulations were performed either with the Kozeny-Carman relationship used previously, or with a correlation obtained from the volume averaging method based on the work of [START_REF] Whitaker | Flow in porous media I: A theoretical derivation of Darcy's law[END_REF]. The unit cell used to solve the closure problem associated with the calculation of the permeability is the same as for the mass transfer coefficient or the dispersion tensor. Comparison between these correlations is represented figure 20. Permeability value K f for the fluid region, i.e. the permeability for a near unity Numerical data used for the simulations are shown in table 3. The core is 25 cm in length and 5 cm in width. The injected acid concentration, C acid , corresponds to a salt concentration of 150 kg m -3 .

A first series of simulations was performed in order to calibrate the mass transfer coefficient with respect to the optimum flow rate. Shape factor A was fixed to 1.1 for the Kozeny-Carman K-ε β relationship, and to 3 for the second correlation. The impact of the K-ε β correlation on the numerical simulations is represented in figure 22 in terms of extension of the dissolution pattern versus pore volume for two injection flow rates. One corresponds to a local equilibrium dissolution (Q = 1 cm 3 h -1 , conical regime), the other corresponds to a local non-equilibrium dissolution (Q = 80 cm 3 h -1 , optimum flow rate). As expected, the K-ε β correlation plays a minor role in local equilibrium dissolution where a sharp front delimits the fluid and porous zones (no transition zone for porosity or permeability). On the other hand, there is a great sensitivity of the model in the local nonequilibrium situation and a precise calibration of the model is therefore necessary.

Numerical simulations were performed with K estimated from the closure problem and A = 3. Comparison between experimental and numerical results in terms of pore volume to breakthrough versus injection flow rate is presented in figure 23. Considering the assumptions made about the pore-scale geometry of the medium, the results are in fair agreement with an accuracy for the time to breakthrough of about 25% in the region near the optimum flow rate. Improvements would require a better characterization of the pore-scale structure, which is beyond the scope of this paper.

Conclusion

A three-dimensional coupled reaction-transport Darcy-scale model has been presented in this paper. Compared to previous dissolution models, its originality lies in the fact that it is based on a combination of a Darcy-Brinkman model and a local non-equilibrium dissolution model. A specific numerical procedure has been developed to solve for the coupled equations. The Darcy-scale model, and the associated 'effective' properties are linked to pore-scale properties. Historical behaviour is disregarded by using simple representative unit cells and a simple pore-scale dissolution framework.

Several numerical simulations have been performed. They show that this mathematical model can reproduce all the dissolution regimes observed in the literature. Not only does the model capture the different dissolution regimes, but it can be used to give a physical interpretation of the dissolution mechanisms. It also can be used to propose quantitative diagrams representing transitions between these different regimes. In addition, important features such as optimum flow rates are recovered, provided that correct local-scale physical characteristics are determined.

New experimental results have been obtained in the saltpack-undersaturated salt solution system. They show the same dissolution patterns as previously published experiments in the acid-limestone system. For a quantitative analysis, it appears that the most sensitive parameter in the model is the mass transfer coefficient, α, since the Damköhler number plays a key role as a controlling parameter for the wormhole propagation. The permeability-porosity relation is also very important, especially for local non-equilibrium regimes.

The link between the mass transfer coefficient and the Damköhler number allowed us to define an optimum Damköhler number, denoted Da opt . By a comparison between numerical and experimental data, we were able to determine the most probable α coefficient, for the system under consideration.

Our Darcy-scale model offers a very interesting tool to understand the physics involved, and to predict the wormhole propagation in a quantitative manner. It may be used to study the effect of many external parameters, such as the variation of temperature or concentration, or process modifications like the use of acids in emulsion, which give lower diffusion effects. Such a model can also be used to obtain information on tricky questions, such as dissolution patterns for porous media of very large extent, since most of the pore-scale physics is incorporated in the effective coefficients appearing in the Darcy-scale description, which are difficult to handle with pure pore-scale models. This could be useful in understanding the transient evolution of unstable dissolution patterns around a wellbore for instance, or in semi-infinite media. Finally, a Darcy-scale model presents a simpler mathematical form, and it can be be used very efficiently for further stability studies for instance, which in turn may provide valuable information. these circumstances, this equation can be replaced by the classical form reported in (2.21).

The hydrodynamic dispersion tensor D * β is assumed to be a function of the interstitial velocity through the dispersivities α L , α T in the longitudinal and transverse directions respectively, and a coefficient of effective diffusion D eff taking into account molecular diffusion and tortuosity of the medium [START_REF] Bear | On the tensor form of dispersion in porous media[END_REF]. It is a second-rank tensor, which for an isotropic medium is often expressed as

D * ij = (α T |V β | + D eff )δ ij + (α L -α T ) V i V j |V β | , ( B 17)
where δ ij is the Kronecker delta function.

The numerical procedure designed for the resolution of the closure problems is described in the literature [START_REF] Quintard | Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF][START_REF] Quintard | Fracture-matrix mass exchange coefficients[END_REF]). As a starting point for a comparison with our experiments, i.e. NaCl grains, a crude representation of the porous medium structure was adopted and is represented by where l β is the size of the unit cell as shown in figure 24. All the results show a diffusive regime followed by a dispersive regime, more pronounced for the longitudinal dispersion coefficient. If we assume that the pore-scale dissolution process corresponds to a uniform dissolution of the solid interface, we can replace the historical dependence of the effective parameters by these correlations with porosity and velocity. Figure 28 shows the influence of the geometry of the unit cell. As shown in this paper, the mass exchange coefficient plays an important role in the quantitative analysis of the wormhole formation, and this is supported by previous reports in the literature. Since we had no means to better characterize our porous medium structure, we scaled the correlation obtained to fit experimental data.

Figure 1 .

 1 Figure 1. The different scales of the problem.

Figure 3 .

 3 Figure 2. Experimental setup.

Figure 4 .

 4 Figure 4. Example of experimental patterns of dissolution: (a) face dissolution, Q = 5 cm 3 h -1 , C NaCl = 150 g l -1 ; (b) conical wormhole, Q = 10 cm 3 h -1 , C NaCl = 230 g l -1 ; (c) dominant wormhole, Q = 50 cm 3 h -1 , C NaCl = 150 g l -1 ; (d ) ramified wormhole, Q = 250 cm 3 h -1 , C NaCl = 150 g l -1 .

Figure 5 .

 5 Figure 5. Dissolution history in the wormholing regime; Ca NaCl = 150 g l -1 , Q = 50 cm 3 h -1 .

Figure 6 .Figure 7 .

 67 Figure 6. Number of pore volumes to breakthrough for a NaCl concentration of (a) 150 g, (b) 230 g l -1 .

Figure 8 .

 8 Figure 8. Representation of the case studied.

Figure 9 .

 9 Figure 9. Porosity fields representative of the dissolution structures obtained numerically. (a) Face dissolution, P e = 8.32×10 -4 , Da = 120; (b) conical wormhole, P e = 4.14×10 -3 , Da = 24; (c) dominant wormhole, P e = 1.66, Da = 6.01 × 10 -2 ; (d ) ramified wormhole, P e = 83.2, Da = 1.2 × 10 -3 ; (e) uniform dissolution, P e = 832, Da = 1.2 × 10 -14 .

Figure 10 .

 10 Figure 10. Concentration fields representative of the dissolution structures obtained numerically; (a-e) same as figure 9.

Figure 11 .

 11 Figure 11. Porosity fields representative of the dissolution history in the wormholing regime.

Figure 12 Figure 13 .

 1213 Figure 12. P e-Da diagram: N ac = 4.66 × 10 -2 , L = 10 cm.

Figure 14 .

 14 Figure 14. Behaviour diagram according to Daccord et al. (1993).

Figure 15 .

 15 Figure 15. Optimum flow rate at C acid = 545 g l -1 and (a) A = 0.1, (b) A = 0.2, (c) A = 0.5.

QFigure 17 .

 17 Figure 16. Effect of α on wormhole propagation.

  Using our numerical data and √ K as the microscopic characteristic length scale l, we can estimate the optimum Damköhler number according to the definition of Fredd

Figure 18 .

 18 Figure 18. Pore volume to breakthrough for different concentrations.

  Figure 19. Concentration effect on the regime diagram; L = 10 cm.

Figure 20 .Figure 21 .

 2021 Figure 20. Comparison of the permeability-porosity correlations; K fluid = 8.33 × 10 -8 .

Figure 22 .Figure 23 .

 2223 Figure 22. Impact of the permeability-porosity correlation on the time to breakthrough for C NaCl = 150 g l -1 in (a) local equilibrium dissolution (Q = 1 cm 3 h -1 ) and (b) local non-equilibrium dissolution (Q = 80 cm 3 h -1 ).

  201 × 101 nodesC acid = 210 kg m -3 ρ σ = 2160 kg m -3 µ = 1 × 10 -3 Pa s ε β = 0.36 β = 1 D β = 2.10 -9 m 2 s -1 K = 1.5 × 10 -11 ± ε m 2

Figure 25 .Figure 26 .

 2526 Figure 25. Mass transfer coefficient α as a function of the Péclet number and for different porosities.

Figure 27 .Figure 28 .

 2728 Figure 27. Transversal dispersion coefficient D * yy .

  the two-dimensional periodic unit cell shown in figure 24. Results for the different effective parameters are represented versus porosity in figures 25, 26 and 27 in terms of the cell Péclet number defined by P e cell = |V β |l β ε

Table 1 .

 1 Numerical data for constant inlet pressure injection.

Table 3 .

 3 Numerical data; constant flow rate injection.

Appendix A. Dissolution and mass balance equations

In this Appendix, we make explicit the relationship between the pore-scale physics and the Darcy-scale model used in our previous work [START_REF] Quintard | Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF], 1999).

Taking the volume average of the fluid mass balance equation, (2.5), and after using the averaging theorem [START_REF] Marle | Ecoulements monophasiques en milieu poreux[END_REF][START_REF] Howes | The spatial averaging theorem revisited[END_REF][START_REF] Gray | Mathematical Tools for Changing Spatial Scales in the Analysis of Physical Systems[END_REF], we obtain

From (2.9) we have n βσ • (v Bβ -w) = 0 at A βσ .

(A 2) Moreover, from the tracer case assumption used for species A and σ, we can write

and the use of (A 3) along with the boundary condition in (A 2) provides the following relation:

(A 4) Finally, (A 1) can be approximated by

At this point, we can make the assumption that the dissolution mechanism is slow enough to be considered as a quasi-steady-state process, and this is consistent with the assumption associated with (A 4). As a consequence, we can couple the Darcy-Brinkman equation with the following relation:

(A 6)

We can now average the acid mass balance equation (2.6)

Applying the general transport theorem and the spatial averaging theorem [START_REF] Whitaker | A simple geometrical derivation of the spatial averaging theorem[END_REF] to this equation, yields

The assumption of a dilute solution for species A leads to

By introducing this relation in (A 8), we can write

Here this equation is coupled with the pore-scale boundary value problem, and in particular with (2.10). This latter equation associated with rapid pore-scale transport may lead to a local equilibrium model characterized by

A local non-equilibrium model has been discussed in the literature, and the reader will find a discussion in [START_REF] Quintard | Dissolution of an immobile phase during flow in porous media[END_REF]. When the change in the interface shape is slow enough, a good approximation of the problem leads to the following form of the local equation:

We are ready now to move on to the dissolution rate of the σ-phase. The upscaling method applied to (2.8) leads to

Since ρ σ is constant, we have ρ σ σ = ρ σ and using the straightforward relationship

We can now make use of this result and the boundary condition on A βσ given by (2.11) to express (A 14) as

in which we recognise the mass exchange flux of species A. Given the approximations made in the theoretical development, we obtain

Appendix B. Correlations for the effective coefficients

The correct estimation of the effective coefficients (α and D * β , equation (2.21)) is a fundamental step in order to properly model the acid propagation into the porous medium. In the absence of experimental data for salt dissolution, we have used in our simulations a numerical estimation obtained from pore-scale calculations using closure problems proposed in [START_REF] Quintard | Convection, dispersion, and interfacial transport of contaminants: Homogeneous porous media[END_REF], 1999). A brief summary of the introduction of these problems is given below.

We introduce Gray's decomposition of velocity and concentration [START_REF] Gray | A derivation of the equations for multiphase transport[END_REF] as

)

) where ṽβ and cAβ are velocity and concentration spatial deviation, respectively.

Use of the spatial decomposition in (A 10) provides [START_REF] Carbonell | Dispersion in pulsed systems II: Theoretical developments for passive dispersion in porous media[END_REF])

At this point we need to make use of (A 1). We multiply this relation by C Aβ and subtract the result from (B 3) to obtain [START_REF] Chella | Multiphase, multicomponent fluid flow in homogeneous and heterogeneous porous media[END_REF])

Finally, use of (A 15) in (B 4) leads to

In order to develop an equation for the spatial deviation concentration, we divide (B 5) by ε β and subtract the result to (2.6) to construct the following result:

The spatial decomposition applied to the boundary condition (2.10) completes the closure problem:

In these equations, non-homogeneous terms can be viewed as source terms for cAβ and this suggests that cAβ can be represented by

where b β and s β are the mapping variables of the following associated problems:

where u β in (B 10a) is given by

)

where the mass transfer coefficient α (B 12a) is defined by

In these equations, we make use of periodicity conditions to complete the closure problems, and this means that the effect of strong gradients close to the macroscale dissolution interface have been neglected.

Using these forms in (B 5) leads to the following macroscopic equation for acid transport: The two last terms on the left-hand side of (B 14) are usually negligible, and under