Abstract: In this study, we derived the effective mass-transfer coefficient between two-phase flow in a three stratified porous medium, one (water) is flowing and the other (NAPL) is trapped. Traditional volume averaging method was used to calculate the large scale effective properties equation describing NAPL dissolution in one-dimensional stratified systems. One limiting case has been investigated corresponding to local equilibrium dissolution with great mass exchange coefficients. It is interesting to note that this case results in a large-scale non-equilibrium dissolution model. The theoretical developments allow predicting the large-scale effective properties from the Darcy scale characteristics in unit-cell. Darcy scale direct simulations were compared favorably to large-scale predictions obtained from the theoretical models.

INTRODUCTION

The protection of groundwater became a major objective since the years 1970. Indeed, the technological, industrial development and agricultural have for major risks the pollution of this water. Certain pollutants discharged into nature such as the toxic products, heavy metals and hydrocarbons constitute a great risk of contamination of groundwaters. Pollution by hydrocarbons comes as well from accidental discharge as of not controlled infiltration of the rejections. When the hydrocarbons penetrate in underground, they infiltrate initially vertically in the unsaturated zone under the effect of gravity. The hydrocarbon source of long-term pollution is in various forms: nonaqueous, dissolved or gas (Abriola and Pinder, 1985a, b). Then when the hydrocarbon is denser than water, it crosses the water saturated zone if not for hydrocarbons lighter than water, it is spread out horizontally in the capillary fringe. The tablecloths fluctuations due to the rainy infiltrations move the sources of pollution between the saturated zone and the unsaturated zone. In the saturated zone, hydrocarbon ganglia are trapped in NAPL form and because of their very low solubility, they can maintain a pollution of long duration. In practice to model the dissolution of the pollutants in porous media, the assumption of local equilibrium since strong a long time was used in the majority of the industrial codes (Abriola and Pinder, 1985a, b; Corapcioglu and Baehr, 1987; Kaluarachchi and Parker, 1989; Powers et al., 1991; Fried et al., 1979; Miller et al., 1990; Parker et al., 1991). It translates the fact that the average concentrations are close to the concentrations of thermodynamic equilibrium on the interface between the two phases on pore scale. However, the recent studies (Mercer and Cohen, 1990; Powers et al., 1991; Miller and Mayer, 1996; Miller et al., 1998; Hunt et al., 1988a, b; Geller, 1990; Razakarisasoa et al., 1989; Imhoff et al., 1994; Borden and Kao, 1992) tend to show that this assumption, even if it is usable for the small scales, moves away from the experiment on the large scales, especially in heterogeneous medium with weak saturation in NAPL or when the aqueous phase reaches high speeds. These limitations can be induced by heterogeneities at pore scale or by medium heterogeneities or saturations at higher scales or by raised circulation speeds of water or by the dispersion state of the organic phase. In fact for a given configuration, there is a scale to which the transfer is well described by the thermodynamic equilibrium of dissolution (Quintard and Whitaker, 1994). This design assumption is not always valid for REV sufficiently large or at the ground scale (Fried et al., 1979). This is why, of work having for objective the installation of the mathematical models to describe the kinetics of dissolution were developed (Fried et al., 1979; Sleep and Sykes, 1989; Miller et al., 1990; Powers et al., 1991; Brusseau, 1992; Quintard and Whitaker, 1994; Imhoff et al., 1994). The determination of the mass exchange coefficient α (T^-1) differentiates the various models. The evaluation of this coefficient is available in the literature either from experimental correlations (Miller et al., 1990; Powers et al., 1991, 1992; Geller and Hunt, 1993; Imhoff et al., 1994; Kennedy and Lennox, 1997) or from theoretical developments with scale changes (Quintard and Whitaker, 1994; Ahmadi and Quintard, 2000; Ahmadi et al. 2001a, b; Coutelieris et al., 2006). But the using of the heterogeneities effect of the porous medium on dissolution still poses problems, particularly when they are the multi-stratified media.

The objective of this research is to model the pollutants dissolution in three-stratified porous medium at large scale using the numerical code conceived by Mabiala for the dissolution of NAPL in two strata porous medium systems. It is a question of determining from the same assumptions adopted by Mabiala, the effective properties and particularly the large scale mass exchange coefficient which appears in dissolution and saturation equations in one dimensional case for three strata porous medium.

MATERIALS AND METHODS

Among the solutions planned to follow the evolution of pollutants in porous media, appears in good place the grounds pollution modeling. To model the pollution of the grounds by hydrocarbon compounds, such is the objective of this research, amounts modeling a transport phenomenon in porous media. To date, studies (Miller et al., 1998) show that in spite of many research completed, this modeling still poses a certain number of questions of which some have not yet a response from a research point of view. The grounds are naturally heterogeneous porous media at all observation scales, the multi-scale aspect remains still a characteristic of the problems of transfers in porous media (Dagan, 1988; Cushman and Ginn, 1993).

As shown in Fig. 1, one distinguishes several fundamental scales: microscopic scale or pore scale; macroscopic scale or of Darcy scale; generally some characteristic dimensions of the pore; large scale or macroscopic heterogeneities scale of the porous medium. For the first two scales, much work was undertaken. Admittedly, the problems all are not solved but acceptable models in comparison with the precise details requested are available (Abriola and Pinder, 1985a, b; Corapcioglu and Baehr, 1987; Kaluarachchi and Parker, 1989; Quintard and Whitaker, 1994; Ahmadi et al., 2001a, b).

Currently, the using of the heterogeneities effect when the scale of description in the model used is too large, still poses many unsolved problems. For example, the application of these models to very heterogeneous real systems can pose problems when the mesh size is too large. The description of the porous medium at large scale is in general determined by the size of the grids used in the numerical models, this new scaling (Darcy scale to large scale) allows to apply physical laws related to the flows multiphasic multiconstituants in porous media and to highlight the equations and the effective properties associated at large scale which control the dissolution of pollutant in porous medium.

Fig. 1:	Multi-scales aspects in modeling transfers in porous media

The passage of the Darcy scale to large scale is done by means of the volume average method which consists in using theorems of spatial average (Anderson and Jackson, 1967; Slattery, 1967; Marle, 1967; Whitaker, 1967; Gray and Kostin, 1975), in order to integrate the pollutant transports equations obtained on a Darcy scale in Representative Elementary Volumes (REV) of a polyphasic system. The method aims at obtaining the equations of dissolution on the large scale, the effective properties associated and the macroscopic boundary conditions.

Considering the assumptions hereafter; the motionless nonaqueous phase, the porosity of the medium not function of time and constant densities of the various phases, dissolution equations of pollutant in the porous medium commonly met in the literature at Darcy scale (Powers et al., 1992; Imhoff et al., 1994; Quintard and Whitaker, 1994) are:

(1)

(2)

Where, C_βγeq is the equilibrium concentration (in mass fraction) of the pollutant in the liquid phase. One sees appearing on these equations, the mass exchange coefficient which is due to the presence of a pollution source. Several researches were undertaken to estimate this mass exchange coefficient at Darcy scale (Miller et al., 1990; Quintard and Whitaker, 1994; Ahmadi and Quintard, 2000; Coutelieris et al., 2006) on simple and complex unit cells. The average values on a large scale are defined as it follows:

(3)

(4)

(5)

Applying these definitions to Darcy Eqs. 1 and 2, one obtains the large scale dissolution equations of a pollutant in the porous medium hereafter:

(6)

(7)

Where:

Which is the form of the effective properties and particularly that of the large scale mass exchange coefficient in the particular case of this study.

Numerical modeling of large scale dissolution
One dimensional case: Numerical modeling is a scientific alternative in the flow and transport fields of pollutants in the subsoil owing to the fact that the experiments are difficult to realize and their costs are very high. The difficulties of this numerical modeling are due particularly to heterogeneities of the calculation field. These heterogeneities appear on several space scales. At large scale, they are caused by the geological layers changes which involve in their turn the variations of the hydrogeologic properties of the rocks which compose these layers, certainly causing variations of mass exchange coefficient.

One-dimensional model: The porous medium (Fig. 2) which we consider is compared to a stacking of parallel layers with constant thickness and various characteristics (porosity and permeability). This medium is indeformable and contains initially two incompressible liquid phases: oil and water. In this research, the unit cell is a whole of three different strata (ω, η, φ) of length l_ω, l_η and l_φ. The medium contains initially partially miscible pollutant at residual saturation S_γr. At x = 0, at the constant velocity V_β, one injects perpendicular to strata, pure water in the polluted porous medium, pure water takes care while polluting and dissolving the trapped phase and comes out again with a certain level of concentration.

Fig. 2:	ID physical model

In the case of local equilibrium, the dissolution front is located by x_f. The system considered is a binary system (two components, water and hydrocarbon and two phases) and we consider that each stratum has respectively a noted constant porosity ε_ω, ε_η and ε_φ. Under these conditions, the equations of the one dimensional problem are:

(8)

(9)

Initial conditions:

(10)

(11)

Boundary conditions:

(12)

(13)

where, C_βγeq, S*_γr and x_L are respectively the NAPL equilibrium concentration in the liquid phase, NAPL residual saturation and the studied field characteristic length. In addition, the concentrations and saturations obey the following relations:

(14)

(15)

(16)

To solve the system of Eq. 8-12, the diffusion coefficient of water phase is supposed to be constant and pure water is injected at x = 0 at constant velocity. However, the large scale mass exchange coefficient α* is to be determined.

Large scale mass exchange coefficient estimation: The mass exchange coefficient which appears in the large scale dissolution equations of pollutant in the porous medium is influenced by certain complex mechanisms: transfer of the pollutant from interface towards the porous medium, geometrical evolution of the interface due to dissolution. Several researches was undertaken to estimate this coefficient (Miller et al., 1990; Powers et al., 1991; Quintard and Whitaker, 1994; Aigueperse and Quintard, 1995; Ahmadi and Quintard, 2000) in simple unit cells. In this study, we consider initially the local equilibrium conditions at Darcy scale in the unit cell with three strata shown in Fig. 3. We establish theoretically, the large scale masss exchange coefficient for a one-dimensional flow perpendicular to the strata (1D-problem). Initially, the system is at equilibrium under conditions of Fig. 3. Pure water injected at the constant filtration velocity, breaks the equilibrium. Thus, three situations can arise: the dissolution front is in the first stratum; the dissolution front is in the second stratum; the dissolution is in the third stratum.

The dissolution front is in the first stratum: The arrival of pure water at x = 0 modifies the conditions in the unit cell which become those shown in Fig. 4. We obtain the average values of S_γr, ε, C_βγ, S_β and S_γ as it follows:

(17)

(18)

(19)

Fig. 3:	Initial conditions

Fig. 4:	Conditions in the unit cell when the dissolution front is in the first stratum

After some simplifications, we obtain what follows:

(20)

(21)

We evaluate the total mass of NAPL contained in the unit cell of Fig. 4 by supposing that the quantity of water which enters NAPL-phase is negligible i.e., C_βγeq<<1 and C_γγeq≅1. This allows to write:

(22)

Aρ_γε_ωS_γ.rωx_fC_γγeq = 0 because the dissolution front is in the first stratum. Equation 22 becomes:

(23)

The mass flux leaving the nonaqueous phase is equal to the NAPL total mass variation contained in the unit cell of volume Vol = A.l (A being the unit section of the unit cell):

(24)

This enables us to determine the dissolution front velocity:

(25)

By identifying NAPL mass variation given by Eq. 24 per the mass exchange term of the dissolution equation, we can write:

(26)

Thus, the mass exchange coefficient takes the form:

(27)

Where:

The expression of the volume of the liquid phase contained in the unit cell is:

(28)

whose development is written:

(29)

The development of Eq. 20 gives:

(30)

By using Eq. 18 in Eq. 30, we obtain:

(31)

Let us multiply Eq. 31 by the unit section A, this equation becomes:

(32)

According to Eq. 32, Eq. 29 becomes:

(33)

Thus, Eq. 27 is still written:

(34)

Equation 19 can be written:

(35)

In other words:

(36)

According to Eq. 31, one obtains:

(37)

Equation 31 allows to obtain according to Eq. 17 and 18, what follows:

(38)

Let us replace the expression of Eq. 38 in Eq. 37, it is written then:

(39)

Equation 34 is written according to Eq. 39 in the following way:

(40)

This large scale mass exchange is obtained if the dissolution front is in the first stratum i.e., when x_f≤l_ω. This condition can be still written:

(41)

According to Eq. 38, Eq. 41 becomes:

(42)

Equation 42 can be written:

(43)

One defines the volume fraction corresponding to each stratum of the porous medium considered by φ_i = l_i/l with i = ω, η, φ. According to this definition, finally condition x_f≤l_ω is written:

(44)

The dissolution front is in the second stratum: The conditions describing this situation are shown in Fig. 5. In the same way that previously, we obtain the following average properties:

(45)

(46)

(47)

The solubility of NAPL in the aqueous phase being weak, the same assumptions are used. Under these conditions, the total mass of NAPL in the unit cell is:

(48)

The mass flux of NAPL outgoing from the nonaqueous phase is equal to the variation of the total mass of NAPL contained in the unit cell.

(49)

(50)

Identifying Eq. 49 with the dissolution equation, one finds the same relation of Eq. 27:

(51)

Fig. 5:	Conditions in the unit cell when the dissolution front is in the second stratum

Equation 51 can be written:

(52)

Where:

The volume expression of water phase contained in the unit cell is:

(53)

(54)

Using Eq. 18, one obtains:

(55)

Equation 17 gives:

(56)

Multiplying Eq. 56 by the unit section A, one obtains:

(57)

Replacing Eq. 57 in Eq. 55, one obtains:

(58)

However, Aε_ωl_ωS_γ.rω = 0 because the residual saturation of the polluting phase in the second stratum is null, Eq. 58 is written then:

(59)

Let us calculate the expression of Aε_ηS_γ.rη(x_f-l_ω). Equation 47 is written according to Eq. 18:

(60)

The development of the Eq. 60 gives:

(61)

Replacing Eq. 61 in Eq. 59, the volume expression of the aqueous phase becomes:

(62)

Let us calculate the expression of C*_βγ/C_{β γeq}^¯1. Equation 45 can be written:

(63)

According to Eq. 47, one obtains:

(64)

Replacing Eq. 64 and 62 in Eq. 52, one obtains:

(65)

Equation 61 can be still written:

(66)

Replacing Eq. 66 in Eq. 65, the large scale mass exchange coefficient for this case has finally the form:

(67)

This coefficient is obtained if the dissolution front is in the second stratum i.e., when l_ω≤x_f≤l_ω+l_η. This condition can be still written:

(68)

Let us consider initially the first part of the double inequation 68 i.e.,

(69)

In addition, according to Eq. 66, one obtains:

(70)

Knowing x_f, the inequation 69 becomes:

(71)

After simplification, one obtains the following condition:

(72)

Now let us consider the second part of the Eq. 70 i.e.:

(73)

While proceeding in the same way as previously, we can write:

(74)

As x_f is the same one (Eq. 70), the Eq. 74 allows to obtain:

(75)

Finally, the correlation (Eq. 67) of the large scale mass exchange coefficient if the dissolution front is in the second stratum is valid only if:

(76)

Dissolution front is in the third stratum: The conditions after injection of water in this case are shown in Fig. 6. In the unit cell of Fig. 6 from the same manner that previously, we obtain the following average properties:

Fig. 6:	Dissolution front is in the third stratum

(77)

(78)

(79)

Using the same assumptions, NAPL total mass contained in the unit cell as shown in Fig. 6 is:

(80)

The mass flux of NAPL outgoing from the nonaqueous phase is equal to the variation of the total mass of NAPL contained in the unit cell:

(81)

Identifying the variation of NAPL mass given by Eq. 81 per the mass exchange term of the dissolution equation, we can write:

(82)

May be:

(83)

Equation 83 can be still written:

(84)

Where:

The expression of the volume of the liquid phase contained in the unit cell is:

(85)

Developing Eq. 85 and using Eq. 18, one obtains:

(86)

While proceeding in the same way that previously, one determines using Eq. 18 and 78, the expression of ε_φS_γ.rφ(l-x_f). One establishes finally the volume expression of the aqueous phase contained in the unit cell:

(87)

According to Eq. 18 and 77, one obtains:

(88)

Finally, using Eq. 87 and 88, the mass exchange coefficient when the dissolution front is in the third stratum has the form hereafter:

(89)

This coefficient is obtained if the dissolution front is in the third stratum i.e., when x_f≥l_ω+l_η. By multiplying this condition by ε_φS_γ.rφ, one obtains:

(90)

In addition, Eq. 18 and 78 allow to have:

(91)

That is to say:

(92)

Finally, the mass exchange coefficient when the dissolution front is in the third stratum is valid only if:

(93)

In summary, the large scale mass exchange coefficient in the unit cell is:

(94)

There are three different determinations of the large-scale mass exchange coefficient depending on the location of the dissolution front (Fig. 4-6) within the unit cell. This result can be extended to any number of strata. It will be noted that the behaviour of this coefficient is the same that carried out by Mabiala in the case two stratified porous medium.

RESULTS AND DISCCUSION

For the numerical simulations, we use the computation code DISSOL 1D conceived to modeling pollution of the grounds by hydrocarbon compounds in heterogeneous porous media by Mabiala. The large scale mass exchange coefficient used in this model is that estimated by initially supposing the field polluted in Darcy scale local equilibrium (Eq. 94). Numerical tests were carried out to test the theoretical model. One uses a sufficiently large mass exchange coefficient so the dissolution front is sharp and characteristic of local equilibrium dissolution.

The simulations were performed with the characteristics shown in Table 1. The concentration and saturation fields obtained from direct simulations and the averaged equations at different time are shown in Fig. 7 and 8. They illustrate the good quality of the suggested theoretical model.

However, we note that the trends of the saturation and concentration fields obtained at large scale are less stiff than those obtained at Darcy scale, this is characteristic of the large scale local non-equilibrium situations.

Table 1:	Data on a Darcy scale and the large scale for the simulations carried out with α* estimated from Darcy scale local equilibrium assumption
Characteristics of study domain: xL = 20 cm

Fig. 7:	Darcy scale and large scale concentrations at a) t = 2500 sec; b) t = 3500 sec; c) t = 5500 sec

Fig. 8:	Darcy scale and large scale concentrations at a) t = 2500 sec; b) t = 3500 sec and c) t = 5500 sec

Fig. 9:	The concentration at the downstream of zone source at the point xL at large scale

Therefore a large scale local non equilibrium model may be needed even if the flow at the Darcy scale features local equilibrium. The concentration at the downstream of zone source at the point x_L shown in Fig. 9 is characterized by a plate at the equilibrium concentration: the zone source is still significant, this is why the concentration is equal to the equilibrium concentration during a certain time; a non-local equilibrium situation illustrated by the reduction of concentration (smaller than equilibrium concentration) to reach a zero value when the totality of the polluting phase was dissolved. It is noted that the disappearance duration of the zone source by complete dissolution is rather long.

CONCLUSION

From the large scale dissolution equations of a pollutant in porous media available in the literature, we carried out a scaling of these equations in order to evaluate the effective properties useful for the computation code worked out by Mabiala such as the coefficient of large scale mass transfer coefficient for one dimensional flow in a porous medium with three strata. One limit case of dissolution was studied the local equilibrium case. This development led to a large scale non-local equilibrium model which allows to predict the large scale effective properties mentioned earlier from the Darcy scale characteristics in a unit cell.

Numerical simulations were carried out using the computation code DISSOL 1D developed by Mabiala. The results of large scale simulations obtained were compared favorably with the results obtained by direct Darcy scale simulations in terms of saturation and concentration fields while polluting in the phase water. In prospects, it would be interesting to study the phenomena of transfer and transport of the pollutants in the porous media including >3 strata to take account of heterogeneities which appear in the real porous media.

ACKNOWLEDGMENT

We would like to thank Mr Jacques Ngoulou for his remarks and observations.