Two-Dimensional Modeling of the NAPL Dissolution in Porous Media : Heterogeneities Effects on the Large Scale Permeabilities and Mass Exchange Coefficient

In this paper, we are interested by the dissolution of NAPL (Non-Aqueous Phase Liquid) contaminants in heterogeneous soils or aquifers. The volume averaging technique is applied to 2D systems with Darcy-scale heterogeneities. A large-scale model is derived from a Darcy-scale dissolution model in the case of small and large Damkhöler numbers, i.e., for smooth or sharp dissolution fronts. The resulting models in both cases have the mathematical structure of a non-equilibrium dissolution model. It is shown how to calculate the resulting mass exchange and relative permeability terms from the Darcy-scale heterogeneities and other fluid properties. One of the important finding is that the obtained values have a very different behavior compared to the Darcy-scale usual correlations. The large scale correlations are also very different between the two limit cases. The resulting large-scale models are compared favorably to Darcy-scale direct simulations.


Introduction
The mechanisms of migration in aquifers of partially miscible pollutants to water (oil, hydrocarbons…), later denoted NAPL for Non-Aqueous Phase Liquid, are well-known.Two-and three-phase flow mechanisms lead to a more or less dispersed distribution of the organic phase in the aquifer.The prediction of the contaminant plume requires in general solving the transport equations completed with proper initial and boundary conditions.In practice, two classes of model are available.The first class makes use of the assumption of local equilibrium, which translates the fact that the average, Darcyscale concentrations are close to the pore-scale thermodynamic equilibrium concentrations at the interface between the two phases.Indeed, concerning the modeling of dissolution (water/NAPL transfer), many laboratory works ( [1] [2] [3] [4]) show situations where thermodynamic equilibrium makes it possible to describe satisfactory the measured data.However, recent studies ( [5]- [15]), show that this assumption, even it is generally usable for small scales, does not fit very well experimental data in the case of large-scale heterogeneous media, or in the case of low NAPL saturations, or when the aqueous phase reaches high velocities.These limitations of the transfer can be induced by heterogeneities at the pore scale, by heterogeneities of the medium or saturations at higher scales, by high velocities of water circulation, by the dispersion of the organic phase, etc.In fact, for a given configuration, there is a scale at which the transfer is well described by the thermodynamic equilibrium assumption [16].Other works, ( [17] [18]) affirm that heterogeneities of the porous medium can modify the trajectory of the pollutants significantly and thus lead to a flow which can generate a local non-equilibrium situation.This assertion is defended by Mayer et al., [19] who estimate that the modeling of dissolution depends on the knowledge of the heterogeneous distribution of the pollutant in the medium considered.The works of ( [20] [21] [22]) maintain that the disorder and heterogeneities of the porous medium are responsible at the pore scale of the spatial variations of the mass transfer and that the geometry of the medium significantly affects residual saturation, the volume of distribution of the ganglion and the interfacial surface between the aqueous phase and the non-aqueous phase.
Recently, experimental works of NAPL dissolution in heterogeneous porous media [23] were carried out to measure the spatial and temporal distribution of globules of pollutants.The results obtained by these authors show that the concentrations of NAPL decrease with the reduction of the globules of pollutant.Numerical works of modeling of dissolution ( [24] [25] [26] [27]) were carried out.These numerical studies show that: 1) the effectiveness of the dissolution of NAPL depends on the heterogeneity of the aquifer which controls the trapping of NAPL; 2) the deviation of the aqueous phase flow resulting from the heterogeneity and the reduction of relative permeability due to the trapping of NAPL prolongs times of dissolution ( [26]); 3) the morphology of the NAPL source zone, determined by heterogeneity of the porous medium controls dissolution, in other words the precise characterization of the location of the trapped clusters is essential to predict the dissolution flux [27]; 4) when the initial residual saturation is discontinuous, pure water entering in the NAPL source zone can come out again with a certain concentration different or non-equilibrium concentration.Even if in the zones with presence of pollutant, dissolution is of the local equilibrium type, the concentration at the exit can be significantly weaker than equilibrium concentration.Local nonequilibrium situations can thus appear following various mechanisms: effect of heterogeneities at the pore scale and effect of heterogeneities of the medium or initial satura-tions ( [24] [28] [29]).
In recent laboratory experiments, Yra [30] studied the influence of microscopic and macroscopic heterogeneities of two types of porous media (natural macroscopic homogeneous sandy media or macroscopic heterogeneous media formed by alternation of sand and sandstone strata) on the dissolution of a pollutant, with a special interest on the mass transfer between the aqueous phase and the NAPL trapped at residual saturation.Her results show the importance of the porous medium microscopic and/or macroscopic heterogeneities and that of the distribution of the pollutant.The evolutions of the saturation and concentration fields she obtained reveal a slow rate of dissolution and she concluded that these experimental results show the need to use a macroscopic scale local non-equilibrium model.From these different evidences, one can advance the following conclusion: the differential dissolution in the various media can lead to zones with low pollutant saturation beside zones very rich in the pollutant phase.Even if, at the small scale, the mechanisms of dissolution can be described by a model of local equilibrium type, it is not necessarily the same at the large scale (block of a numerical model including heterogeneities).At this scale, water polluted at equilibrium mixes with non-contaminated water or water polluted a non-equilibrium, and the result gives a concentration smaller than the equilibrium concentration at the exit of the source zone, which is typical of a local non-equilibrium situation.It is therefore important to be able to derive a large-scale NAPL dissolution model from the Darcy-scale description, including the heterogeneity effects.
The one-dimensional problem was the subject of a first study [31].In this paper, the two-dimensional problem is considered: a large-scale NAPL, local non equilibrium dissolution model is built using a volume averaging technique and effective properties such as the large scale relative permeability to water and the large scale mass exchange coefficient are estimated for situations of the study.

Physical Large Scale Dissolution Model
In this paper, we consider a stratified model of a heterogeneous porous medium typical of geological systems.The medium is composed of stacks of the elementary pattern described in Figure 1, and, because of periodicity, we will consider only the evolution of the different physical variables within this two-stratum elementary pattern.The media Two different strata ω and η are filled initially with water and NAPL.The typical system is made of two different strata (ω, η) of the same length l.Each stratum contains initially NAPL at residual saturation orw S and or S Pure water is injected at the domain inlet, on the axis (0, y), with a constant Darcy velocity.It is contaminated in pollutant while dissolving the trapped NAPL, and goes out of the domain at a certain concentration, which is one of the important value to be determined by the model for risk analysis.The dissolution front is located by xf.The studied system is a binary system for which we adopt the following notations: the two phases are water, w, and oil, o, and contain each one the two components water (w) and NAPL (o).The strata have a constant porosity ε ω , ε η and constant permeability k ω , k η .
We do not present here the development of the Darcy-scale dissolution model since very detailed and complete developments are available in the literature ( [8] [16] [32]- [37]).Two particular Darcy-scale cases are considered: local equilibrium situation and local non-equilibrium situations.Large scale model equations were established in a previous study [38] as follows: Dissolution equations: ) Large scale Darcy law: Boundary conditions are: at the north and south boundaries of the domain, ( ) ( ) at the east boundary, ( ) ( ) at the west boundary one imposes, the Dirichlet condition as well for the pressure as for the concentration: lated in the case of a stratified heterogeneous system [36] by: with φ i , the volume fraction and K i the permeability of the i-stratum respectively, I the tensor identity and j the unit vector along the y-axis.The resolution of Equations ( 1) to (6'') requires the knowledge of the evolution of the large scale average mass exchange coefficient, a * , and of large scale average dispersion tensor like that of the elements of rw * k .Because of weak variations of saturation, one neglect the dispersion phenomena, one deals with in this work only the closing mass exchange coefficient problems and that of the large scale average relative permeability.

Estimate of the Large Scale Average Effective Properties
The mass exchange coefficient in the transport equation for the pollutant is a simplified representation of complex pore-scale mechanisms: diffusion of the pollutant at the interface towards the porous media, coupling with convection, and geometrical evolution of the interface due to dissolution.Several works were undertaken to estimate this coef- with two strata (Figure 2 and Figure 3), we theoretically build the large scale average mass exchange coefficient.We consider two cases: 1) the cases of a unit cell with strata laid out perpendicular to the flow and 2) the case where the flow is parallel to the strata.

Local Equilibrium Hypothesis: Estimate of Large Scale Average Mass Exchange Coefficient
In previous works ( [29] [31] [38]) on the NAPL dissolution in stratified porous media, the estimates of the mass exchange coefficient was made for flow perpendicular to the strata presented above in two particular cases: 1) very large Damkhöler numbers and 2) Figure 2. Darcy scale conditions when the dissolution front is in the first stratum.
very small Damkhöler numbers.In this work, we extend these investigations to the cases of flow parallel to the strata.One presents below the calculation for the Darcy-scale local equilibrium cases.It is supposed here that the porous medium is initially at Darcy-scale local equilibrium in the two strata illustrated in Figure 4.One injects pure water parallel to the strata.The typical situation is as follows: the velocity of dissolution V f is higher in one of the strata.Here one will suppose that the dissolution front is faster in the ω-stratum.If a certain volume of stratum is taken, one has the two following situations described below.
• Here we assume that the dissolution front is in the first stratum of the unit cell and there is not yet dissolution in the 2nd stratum.This situation is illustrated by Figure 2. We calculate the large scale effective properties and the mass of pollutant contained in the cell in order to estimate the evolution law of the large scale average mass exchange coefficient which controls dissolution in this cases.The large scale effective properties are:  ( ) ) If we suppose, wo C small enough and 1 ooeq C ≈ , which is compatible with the usual data for NAPL contaminants, one can then calculate the mass of NAPL contained in the unit cell ( ) The variation of this mass of NAPL in the cell is: where is the velocity of the dissolution front in the first stratum.This variation of mass is equal to the mass flux of NAPL going out of the unit cell through section A (Cf. Figure 2), i.e.: ( ) where w V * is the filtration velocity of the liquid phase and 16) and ( 17), one deduces the velocity of the dissolution front V f : Following the equation of dissolution, one writes: where v w is the volume of the water phase.Combining Equations ( 17), ( 18) and ( 19), one still writes: ( ) However, the volume of water in the cell is evaluated as follows ) ( ) In addition, according to (11) and ( 13) one has ( ) Moreover, according to ( 8) and ( 13), one can express or or Combination of the relations ( 13), ( 23) and ( 24) in (22) gives finally This correlation is very different from those available in the literature ( and is valid for o S * satisfying the condition: • A situation when the dissolution front has left the first stratum and dissolution starts in the second stratum.This second situation is illustrated in Figure 3.We proceed in the same way that previously and calculation of the effective properties gives The mass of NAPL contained in the unit cell, assuming small wo C and This mass varies in the cell as follows: ( ) ( ) the velocity of the dissolution front in the 2nd stratum.This variation is equal to the mass flux of NAPL going out of the unit cell through section A (Cf. Figure 3).From ( 17) and ( 31), one deduces velocity f V ( ) By the same argument as in the preceding case, one obtains the correlation of mass exchange α * for this situation ( ) In this case, the volume of the water phase can be calculated as follows Moreover, ( ) And ( ) Taking account Equations ( 35) and (36), Equation (33) becomes

{ }(
) This correlation is valid for o S * satisfying the condition below: Ultimately, one finds the same relations as in 1D [29], and when the medium is ho-

Local Non Equilibrium Hypothesis for 2D Darcy-Scale Problems: Calculation of the Large Scale Permeabilities
The problem considered is that to estimate the relative permeabilities during the dissolution of NAPL in a stratified porous medium.The flow, parallel with the strata (Cf. Figure 5), is initially at residual oil saturation.
The coefficient of mass exchange in this case is available in the literature [38].We present here an estimate valid for all dimension of the studied system.For very small Damkhöler numbers, we suppose a certain evolution of saturation at large scale, and we build a large scale average correlation of α * in following differential form: 2) we impose an increment of average dissolution, i.e.,   3).The oil phase being trapped, the flow in each stratum is comparable to a monophasic flow, but for the fact that the permeability will change with the evolving oil saturation.The problem is thus identical to that schematized on Figure 5.Each stratum is characterized by the permeability tensor K wω and K wη .The scaling of the problem makes it possible to homogenize the system with the large scale averaging volume V ∞ and to estimate the average large scale { } rw Kk .Neglecting gravity effects, we write as follows the local scale total flow problem: The scaling of Darcy law for heterogeneous systems is largely developed by Quintard et al., ([40] [41]) and this is known to lead to a large-scale Darcy's law, which we write as: where w * V and w P * are respectively the large-scale filtration velocity and pressure of the water phase and the large-scale permeability tensor of the water phase may be written under the form [38]: In which * K is the average permeability tensor at large scale and rw * k the matrix of the relative permeabilities at large scale whose elements depend on saturation.
Initially the stratified system satisfies the conditions of Figure 4.By adopting the assumption of small Damkhöler numbers, we calculate using the algorithm deduced from approximation (40) rewritten below, the effective properties such as the relative permeabilities and the large scale mass exchange coefficient.
We present in Table 1, the iterative step of calculation which we use in the numerical model.
In Table 1, o S * is the large-scale oil saturation, ow S and o S η are respectively oil saturations in the stratum ω and in the stratum η for this large-scale oil saturation α ω and α ωη the corresponding values of the mass exchange coefficients in the stratum ω and in the stratum η respectively, α * the large scale mass exchange coefficient calcu- lated with approximation of the system of equations 50, k rwω and k rwη the relative permeabilities of water in each stratum, w K * the large scale average permeability of water for the domains ω and η, rw k * is the large-scale water relative permeability and φ ω and φ η are the domain volume fractions such as φ ω + φ η = 1 with φ ω = V ω /V ∞ and φ η = V η /V ∞ .
For a flow which is parallel to the strata, all these characteristics are calculated for the second iteration by the formulas below.
( )   with * K the average permeability of the stratified medium (here the flow is parallel to the strata).Iteration 1 corresponds to the initial conditions in the two strata.We continue the iterations by going each time at the step 1 as indicated by the algorithmic approach above used for the calculation of Moreover, the relative permeabilities at large scale depend on saturation according to Equation (51).In conclusion in this paragraph, it is interesting to stress that the correlations on the large scale mass exchange coefficient obtained in this study differ from the traditional correlations of mass transfer available in literature ( [4] [8] [19]), even for the large scale systems, and who generally take the form where n < 1 and 0.6 ≤ m ≤ 0.9.

Numerical Experiments
With these estimates of the effective properties, the large-scale equations of dissolution are now under a closed form, i.e., purely expressed in terms of large-scale variables.
They are solved numerically in this section for the particular situations considered in

Calculations for Darcy-Scale Local Equilibrium
The numerical experiments were carried out using the mass exchange coefficient of the relations ( 25) and (37).This coefficient corresponds to the dissolution of a NAPL in a stratified system characterized by a unit cell including two strata and having each one, a rather large mass exchange coefficient so that the dissolution front is sharp and characteristic of a local equilibrium dissolution.Simulations at the two scales were carried out with the data of Table 2.
• Concentration fields at large scale Water injection at x = 0, under an entry pressure of Pe = 10 5 Pa in the direction of the flow (the pressure at the exit of the domain is taken to be zero), starts the dissolution phenomenon.Figure 6 illustrate the large-scale averaged concentration fields (in mass fraction) obtained at various times.One observes that: the zone of dissolution increases with time.In the interior of this zone which separates pure water from polluted water, the concentrations are lower than the equilibrium concentration; the concentration at the downstream of this zone remains equal to the equilibrium concentration.Moreover, water with weak concentration in pollutant advances in front of the dissolution front in the strata removed from pollutant in order to contribute to decrease the concentration of pollutant in water polluted in the downstream; the process is very long.
• Large scale saturation fields Figure 7 illustrate the averaged saturation fields obtained at various times.The pressure of entry is one bar and that of exit considered null.One notices in the case of the saturation fields the same behaviors as for the concentration fields.Figure 8 and Figure 9 illustrate the evolutions 2D of the average fields and those at Darcy-scale respectively of the concentration (in mass fraction) and of saturation at

Calculations for Darcy-Scale Local Non-Equilibrium
In this paragraph, we presents the results of the 2D calculations carried out by supposing that the Damkhöler numbers are very small; i.e., the mass exchange coefficient used is that developed in paragraph 3.2.The modeled situation corresponds to the dissolution of NAPLs in a stratified porous medium (Cf. Figure 5).The flow is parallel to the strata.In each stratum, dissolution is controlled by a rather small mass exchange coefficient so that the dissolution fronts are of local non-equilibrium type at Darcy scale.Simulations were carried out with the data of Table 3.
The experiments carried out made it possible to highlight the evolution following the     is higher in the stratum η least permeable.This situation explains certainly the difficulties of starting of dissolution due to the importance of heterogeneities of the medium.
• Evolution of large scale relative permeabilities We illustrate in Figure 11 the evolution of the relative permeabilities to water at the two scales.It is represented on this figure, the relative permeabilities to water for each stratum k rwω and k rw and the large scale relative permeability rw k * .In addition, the large scale permeability is near the permeability of the most permeable medium.
• Effects of heterogeneities on the large-scale effective properties A fine analysis of the effects of heterogeneities of the medium on the large-scale effective properties is essential.Indeed, in the case being studied, NAPL mass flux are controlled at the same time by the differential distribution of the residual pollutant in the medium and by the heterogeneity.
The essential effective properties that we established in this study to model NAPL dissolution in heterogeneous systems are the large-scale mass transfer coefficient and the large-scale effective permeability (Cf.Equations and system of Equations ( 25), (37),  52)).We examine in this paragraph, how the heterogeneities of the medium influence these parameters.
-Effects of heterogeneities on the large-scale relative permeabilities Figure 12 represents the evolutions of the large scale relative permeabilities as a function of saturation for various values of the ω-region volume fraction 0.3, 0.5, 0.7 ω φ = .
( ) φ φ + = and the same porosity (Cf.Table 2).This figure shows that the relative permeability is more significant at the end of the dissolution in the medium containing less pollutant, whereas at the beginning of dissolution, it is rather  the medium containing more NAPL which is more permeable with water.
-Effects of heterogeneities on the large-scale mass transfer coefficients Figure 13 illustrates the evolution of the mass exchange coefficients as a function of the oil saturation, for various values of the ω-region volume fraction, simulated using the algorithm described in Section 3.2.It is observed that the mass exchange coefficient increases with average saturation.
The transfer is more significant when the volume fraction of the stratum ω increases, this means increase of mass of pollutant according to the data of Table 4.This result confirms those obtained previously in the case of the "local equilibrium" assumption for the flows perpendicular to the strata ( [31] [38]), although the correlations of the mass exchange coefficients used for each case are very different.In addition, these mass exchange coefficients grow in both cases with average saturation for the three volume fractions of region considered.
As in the former studies [38], one notices that the mass transfer during the dissolution of NAPL in heterogeneous porous medium is amongst other things controlled by: quantity of residual pollutant contained in poral space; installation of this pollutant in the porous medium, certainly marked by the influence of heterogeneities or hydrodynamic instabilities.

Conclusions
In this paper we have proposed a large-scale model describing NAPL dissolution for stratified heterogeneous systems.
Two limit dissolution cases were studied: a local equilibrium case and a local nonequilibrium case.Both situations lead to a large-scale model of local non-equilibrium type.The construction of the large-scale effective parameters in the local equilibrium case is specific of flows parallel to a stratified system.While the results for the local non-equilibrium case is not specific to 2D or 3D systems, nor a given type of heterogeneity, our application was also presented for flows parallel to a stratified system.
The proposed models have been validated by a comparison of the results of direct Darcy-scale simulations to homogenized simulations.Our results confirm the interest of using local non-equilibrium models for large-scale simulations of contaminated aquifers.In addition, the correlations for the obtained large-scale effective parameters, at least for the 2D stratified systems investigated here, are rather different from most of the correlations suggested by laboratory experiments at, therefore, a small scale.This should be taken into account when designing large-scale models for practical situations.
In perspective, it would be interesting to carry out 2D/3D calculations for more complex heterogeneous media such as nodular media.We believe the results for small Da numbers can be extended in a straightforward manner.However, it must be understood that our results for the case of large Da numbers, i.e., sharp fronts propagating through the system, cannot be in principle extended to more general heterogeneities.V unit cell volume.
x position vector locating the centroid of the averaging volume.
the above equations, a * is the large scale mass exchange coefficient and rw * k is the relative permeability tensor of the w-phase, function of the NAPL saturation.
ficient ([2] [4] [10] [16] [32] [39]) on simple and complex unit cells representative of the pore-scale geometry, and the results obtained lead often to very different estimations.These uncertainties show that the estimate of the mass exchange coefficient remains an open research question.From the Darcy-scale characteristics in a unit cell

Figure 3 .
Figure 3. Darcy scale conditions when the dissolution front is in the second stratum.

Figure 4 .
Figure 4. Initial conditions in the polluted porous medium with parallel flow with the strata.
* + using Equation (40); 6) we reiterate by going to the step 1. However a scaling of the Darcy law is essential to close Equation ( a more general way, if one gives oneself a certain evolution at local Darcy scale of saturation at large scale, of saturation in each stratum, of the coefficient of exchange in each stratum, of the average exchange coefficient at large scale and of the relative permeability of Corey type, all the relations (50) can be evaluated at the time step n + 1.The relative permeability and the large scale mass exchange coefficient are respectively estimated by the following expressions: this work.A comparison is made between Darcy-scale and large-scale results in order to validate the proposed models.The calculations were carried out using a proprietary program DISSOL 2D conceived for this study and based on a classical finite volume formulation [38].The flow is parallel to the strata.The numerical results relate the two theoretical situations analyzed in this work, namely: 1) assumption of local equilibrium at Darcy scale; 2) assumption of local non equilibrium at Darcy scale.

Figure 6 2 .
Figure 6.Concentration fiels wo C * (in mass fraction) at fields at Darcy scale have a very complex spatial evolution certainly due to local heterogeneities strongly present at this scale.In addition, these figures show that the averaged fields obtained from a local equilibrium situation at Darcy scale, are characteristic of a local non equilibrium situation.

Figure 8 .
Figure 8. Evolution of the averaged concentration fields (in mass fraction) and those at Darcy scale at t = 5•10 7 s.

Figure 9 .
Figure 9. Evolution of the averaged saturation fields (in mass fraction) and those at Darcy scale at t = 5•10 7 s.

Figure 10 .
Figure 10.Evolution of the mass exchange coefficients as a function of large-scale saturation.

Figure 11 .
Figure 11.Evolution of the relative permeabilities as a function of large-scale saturation.

Figure 12 .
Figure 12.Evolution of the large scale relative permeabilities as a function of saturation for different volume fractions of the stratified porous medium.

Figure 13 .
Figure 13.Evolution of the large scale mass exchange coefficient as a function of saturation for different volume fractions of the stratified porous medium.

woCCDDU
averaged mass fraction of o-species in the w-phase ( ) equilibrium mass fraction of the pollutant species at the w o − interface.C  deviation of averaged concentration.D dispersion tensor.βdispersion tensor for the β -phase ( ) large dispersion tensor.S β saturation of the β -phase.S  deviation of averaged saturation.orS NAPL residual saturation.NAPLMtotal mass of NAPL in the unit cell.βV average filtration velocity of the β -phase, a vector ( ) average interstitial velocity of the β -phase, a vector ( )1 m s − ⋅ .U  deviation of averaged velocity.fiv dissolution front velocity in i-region ( ) , i ω η = .Pe cell Péclet number.Da Damkhöler number.c l unit cell length.

v
volume of liquid phase in the unit cell.

K
permeability tensor [L2].ri k relative permeability of i-stratum ( ) , i ω η = .ri k large scale relative permeability of i-stratum ( ) , i ω η = .i K permeability of i-stratum ( ) , i ω η = .y position vector locating points in one phase relative to the centroid of the averaging.

Table 1 .
Iterative step for calculation of large scale average characteristics.

Table 3 .
Data at Darcy scale and large scale for the simulations carried out by applying the hypothesis of local non equilibrium at Darcy scale.