Modélisation du transport réactif en milieux poreux non saturés avec une méthode ELLAM en maillage variable

Anis Younes

doi:10.1016/j.crte.2004.01.003

Géosciences de surface (Hydrologie–Hydrogéologie)

Modélisation du transport réactif en milieux poreux non saturés avec une méthode ELLAM en maillage variable
[Modelling reactive transfer in unsaturated porous media using a moving-grid FV-ELLAM]

Presented by: Ghislain de Marsily

Anis Younes ^1,²

¹ Laboratoire de génie industriel, université de la Réunion, 15, av. René-Cassin, BP 7151, 97715 Saint-Denis cedex 09, La Réunion, France
² Institut de mécanique des fluides et des solides de Strasbourg, UMR ULP–CNRS 7507, 2, rue Boussingault, 67000 Strasbourg, France

Comptes Rendus. Géoscience, Volume 336 (2004) no. 6, pp. 547-552.

Abstracts

French
English

La modélisation du transport réactif avec des processus biologiques, chimiques ou radioactifs nécessite des méthodes numériques capables de reproduire des fronts de concentration très abrupts. Dans ce travail, nous présentons une nouvelle méthode ELLAM (Eulerian–Lagrangian Localized Adjoint Method) pour résoudre l'équation de transport réactif avec des coefficients non constants. Afin d'éviter les interpolations (source d'erreurs), la méthode utilise des maillages variables pour définir la solution ainsi que les fonctions test. Cette méthode est utilisée pour simuler la propagation d'un front raide en milieu poreux non saturé ainsi que le transport multi-espèces. Les résultats ne présentent, ni oscillations, ni diffusion numérique, et restent précis, même en utilisant de grands pas de temps.

Modelling contaminant transfer with biological/chemical/radioactive processes needs appropriate numerical methods able to reproduce sharp concentration fronts. In this work, we develop a new Eulerian–Lagrangian Localized Adjoint Method (ELLAM) for solving the reactive transport equation with non-constant coefficients. To avoid interpolation (leading to errors), we use a moving grid to define the solution and test functions. The method is used to simulate first the infiltration of solute into a column of unsaturated porous medium and second the multispecies transport. The developed ELLAM gives accurate results without non-physical oscillations or numerical diffusion, even when using large time steps.

Metadata

Received: 2003-12-30
Accepted: 2004-01-12
Published online: 2004-04-01

DOI: 10.1016/j.crte.2004.01.003

Mots-clés : sol, milieux poreux, écoulements non saturés, ELLAM, transport réactif
Keywords: soil, porous media, unsaturated flow, ELLAM, reactive transport

Author's affiliations:

Anis Younes ^{1, 2}

¹ Laboratoire de génie industriel, université de la Réunion, 15, av. René-Cassin, BP 7151, 97715 Saint-Denis cedex 09, La Réunion, France
² Institut de mécanique des fluides et des solides de Strasbourg, UMR ULP–CNRS 7507, 2, rue Boussingault, 67000 Strasbourg, France

@article{CRGEOS_2004__336_6_547_0,
     author = {Anis Younes},
     title = {Mod\'elisation du transport r\'eactif en milieux poreux non satur\'es avec une m\'ethode {ELLAM} en maillage variable},
     journal = {Comptes Rendus. G\'eoscience},
     pages = {547--552},
     publisher = {Elsevier},
     volume = {336},
     number = {6},
     year = {2004},
     doi = {10.1016/j.crte.2004.01.003},
     language = {fr},
}

TY  - JOUR
AU  - Anis Younes
TI  - Modélisation du transport réactif en milieux poreux non saturés avec une méthode ELLAM en maillage variable
JO  - Comptes Rendus. Géoscience
PY  - 2004
SP  - 547
EP  - 552
VL  - 336
IS  - 6
PB  - Elsevier
DO  - 10.1016/j.crte.2004.01.003
LA  - fr
ID  - CRGEOS_2004__336_6_547_0
ER  -

%0 Journal Article
%A Anis Younes
%T Modélisation du transport réactif en milieux poreux non saturés avec une méthode ELLAM en maillage variable
%J Comptes Rendus. Géoscience
%D 2004
%P 547-552
%V 336
%N 6
%I Elsevier
%R 10.1016/j.crte.2004.01.003
%G fr
%F CRGEOS_2004__336_6_547_0

Anis Younes. Modélisation du transport réactif en milieux poreux non saturés avec une méthode ELLAM en maillage variable. Comptes Rendus. Géoscience, Volume 336 (2004) no. 6, pp. 547-552. doi : 10.1016/j.crte.2004.01.003. https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.1016/j.crte.2004.01.003/

Original version of the full text

The full text below may contain a few conversion errors compared to the version of record of the published article.

Abridged English version

Modelling contaminant transfer with biological/chemical/radioactive processes needs appropriate numerical methods able to reproduce sharp concentration fronts. In this work, we are interested in the numerical solution of the non-constant-coefficient transport equation (1).

For the convective-dominated transport, Eulerian methods often generate numerical solutions with artificial diffusion and/or non-physical oscillations [6]. Characteristic methods use a Langragian treatment of the advective part and an Eulerian treatment of the dispersive part of the transport equation [1,7]. The Eulerian–Lagrangian Localized Adjoint Method (ELLAM) is an improved characteristic method [4]. ELLAM uses space–time test functions [2,10]. Conventional ELLAM formulations use a fixed grid and need interpolation at each time step that is known to add numerical diffusion to the solution [8,9]. To avoid this problem, we develop a new Finite-Volume Eulerian–Langragian Localized Adjoint Method using a moving spatial grid to define the solution and the test functions of ELLAM.

The weak formulation of Eq. (1), using space–time test functions in the form ω=ω(x,t) is given in Eq. (2). ELLAM selects ω(x,t) so as to make the adjoint integral of (3) vanish for all x and t. Therefore, ω(x,t) must be exponential along the characteristics. Hence, the transport equation simplifies to Eq. (4).

In this work, we use a moving-grid method. Therefore, subdomains are not fixed in space (see [14] for details). To evaluate the performance of the developed method, we simulate first the infiltration of solute into a column of unsaturated porous medium. We perform three simulations with different time steps to simulate the non-constant-coefficient transport equation (see Table 1). All simulations give results without any oscillation or numerical diffusion (see Fig. 1).

	Δt	Nombre de Δt	p	Nombre final d'éléments	Temps de calcul [s]
Simulation 1	1	1000	1	1020	17
Simulation 2	100	10	1	30	0,16
Simulation 3	1000	1	6	26	0,14

Fig. 1
Profil de concentration après une durée d'infiltration de 1000 s.
Simulated concentration at t=1000 s.

The method is then used to simulate multispecies transport (Eq. (7)). The three species modelled are ammonium, nitrite, and nitrate, respectively. Good agreement between analytical and ELLAM solutions is obtained (see Fig. 2).

Fig. 2
Simulation du transport multi-espèces.
Simulated concentrations for the multispecies transport problem.

The new ELLAM using a moving spatial grid to define the solution and test functions is robust and well adapted for solving reactive transport equations in porous media.

1 Introduction

La contamination des eaux souterraines par des solutés avec processus biologiques, chimiques ou radioactifs est un problème d'environnement important. Dans ces cas, le modèle mathématique qui décrit le transport de contaminants dans le sol est formé par un ensemble d'équations de convection–dispersion–réaction, qui peuvent être couplées à travers les termes de réaction (réaction entre plusieurs espèces).

Ces processus modifient les fronts de concentration, qui peuvent devenir raides. De même, quand la convection est dominante, la solution présente des fronts de concentration abrupts. Dans ces cas, la solution obtenue par les méthodes classiques (différences finies ou éléments finis) présente des oscillations non physiques [6]. Les différences finies amont ou les éléments finis de Petrov–Galerkin permettent d'éviter ces oscillations, mais génèrent un étalement du front [13] qui aboutit à une mauvaise prédiction de la zone de contamination.

Les méthodes Euleriennes–Lagrangiennes (ELM) telles que MMOC (Modified Method of Characteristics) ont rencontré un certain succès ces dernières années. En combinant une méthode lagrangienne (pour la convection) et une méthode eulerienne (pour la dispersion), elles ont notamment permis d'obtenir des solutions sans oscillations et avec une diffusion numérique très limitée [1,7]. Les principaux inconvénients des ELM sont qu'elles ne conservent pas la masse et ne traitent pas de façon rigoureuse les conditions aux limites [13].

Les ELLAM [2,4,10] constituent une classe de méthodes de caractéristiques qui permettent de prendre en compte les conditions aux limites en utilisant des fonctions test qui dépendent de l'espace et du temps.

Comme les autres ELM, les ELLAM nécessitent des interpolations pour trouver la concentration au pied de chaque caractéristique et à chaque pas de temps [2]. Ces interpolations peuvent engendrer de la diffusion numérique [8,9]. Afin d'éviter ces interpolations, nous avons développé [14] une nouvelle approche, où les fonctions test ainsi que la solution sont définies sur une discrétisation spatiale qui se déplace avec le fluide.

Après la présentation du modèle mathématique de transport réactif en milieu poreux, les grands principes des ELLAM sont exposés. Le modèle est ensuite utilisé pour simuler les expériences de Touma et Vauclin [11] pour le transport en milieu poreux non saturé ainsi que pour le transport multi-espèces (ammonium, nitrite et nitrate) dans un sol.

2 Le modèle mathématique à une dimension d'espace

Le transfert de polluant avec des processus biologiques, chimiques ou radioactifs peut être décrit par un ensemble d'équations d'advection–dispersion–réaction, de la forme :

L (C) = \frac{\partial (R (x, t) C (x, t))}{\partial t} + \frac{\partial (q (x, t) C (x, t))}{\partial x} + K (x, t) C (x, t) - \frac{\partial}{\partial x} (D (x, t) \frac{\partial C (x, t)}{\partial x}) = f (x, t)

(1)

avec q(x,t) la vitesse de Darcy (LT⁻¹), C(x,t) la concentration, L(C) l'opérateur différentiel, R(x,t) le facteur de retard et K(x,t) le coefficient de réaction de premier ordre (T⁻¹).

Le coefficient de dispersion D(x,t) est donné par D(x,t)=λ_L|q|+D_mτ, où λ_L est la dispersivité du milieu (L), τ la tortuosité, D_m la diffusion moléculaire (L²T⁻¹) et f(x,t) le terme puits/source.

L'équation (1) est définie sur le domaine espace–temps, $Ω_{xt} = [0, l] \times [0, T]$ avec des conditions initiales et des conditions aux limites de type Dirichlet ou Neumann.

3 Résolution numérique de l'équation de transport

La méthode ELLAM utilise des fonctions test qui dépendent de l'espace et du temps ω=ω(x,t) et qui sont définies uniquement sur l'intervalle [tⁿ,tⁿ⁺¹]. L'écriture variationnelle de (1) donne alors :

\underset{0}{\int^{T}} \int_{Ω} \frac{\partial (R ω C)}{\partial t} d x d t - \underset{0}{\int^{T}} \int_{Ω} C [R \frac{\partial ω}{\partial t} + q \frac{\partial ω}{\partial x} - K ω] d x d t + \underset{0}{\int^{T}} \int_{Ω} \frac{\partial}{\partial x} [(qC - D \frac{\partial C}{\partial x}) ω] d x d t + \underset{0}{\int^{T}} \int_{Ω} D \frac{\partial C}{\partial x} \frac{\partial ω}{\partial x} d x d t = \underset{0}{\int^{T}} \int_{Ω} f ω d x d t

(2)

On choisit ω de telle façon que l'adjoint associé à la partie hyperbolique soit nul [10] :

\frac{D ω}{Dt} = \frac{\partial ω}{\partial t} + V \frac{\partial ω}{\partial x} = λ ω

(3)

avec V=q/R et λ=K/R.

ω est une fonction qui varie d'une façon exponentielle le long de la caractéristique. L'équation (2) devient alors :

\underset{0}{\int^{T}} \int_{Ω} \frac{\partial (R ω C)}{\partial t} d x d t + \underset{0}{\int^{T}} \int_{Ω} \frac{\partial}{\partial x} [(qC - D \frac{\partial C}{\partial x}) ω] d x d t + \underset{0}{\int^{T}} \int_{Ω} D \frac{\partial C}{\partial x} \frac{\partial ω}{\partial x} d x d t = \underset{0}{\int^{T}} \int_{Ω} f ω d x d t

(4)

La résolution de cette équation nécessite la connaissance de la concentration au pied de chaque caractéristique [2]. Afin d'éviter les interpolations (les pieds des caractéristiques ne correspondent pas forcément aux nœuds du maillage), nous utilisons une discrétisation spatiale mobile qui se déplace avec la vitesse du fluide q/R (voir [14] pour les détails).

4 Simulation de cas d'école

4.1 Transport en milieu poreux non saturé

Nous simulons l'expérience de Touma et Vauclin [11]. Il s'agit d'une colonne de sol verticale, initialement non saturée, de longueur 20 cm, dans laquelle on considère une infiltration continue d'eau et de soluté à partir de la surface. Le milieu poreux est initialement en équilibre avec une pression imposée de −100 cm d'eau en bas de la colonne. Le flux à la surface est de $8, 3 cm h^{- 1}$ .

Les coefficients de l'équation de transport dans ce cas ne sont pas constants, puisque la vitesse et la teneur en eau varient dans l'espace et dans le temps. Pour obtenir la répartition de la vitesse et de la teneur en eau au cours du temps, nous devons préalablement résoudre l'équation de l'écoulement. Le domaine est discrétisé en 20 éléments identiques et l'équation de Richards est résolue avec la méthode des volumes finis en utilisant un pas de temps fixe de Δt=1 s. La teneur en eau est reliée à la pression par la relation [12] :

θ = (θ_{s} - θ_{r}) / {|1 + {(α h)}^{n}|}^{1 - 1 / n} + θ_{r}

(5)

avec [11] : θ_s=0,312, θ_r=0,0265, α=0,044 cm⁻¹ et n=2,2 et la conductivité (en cm h⁻¹) est approchée [11] par :

K = 18 130 θ^{6, 07}

(6)

La dispersivité longitudinale est λ_r=0,1 cm. La durée simulée est de 1000 s. Trois simulations (voir Tableau 1), avec des pas de temps différents (Δt=1 s, 100 s et 1000 s), sont effectuées.

À chaque pas de temps, p nouveaux éléments arrivent par la frontière amont et r éléments quittent le domaine par la frontière aval. Le nombre p peut être fixé par l'utilisateur. Le nombre r est nul dans le problème étudié, car la vitesse de l'eau en bas de la colonne est nulle pendant toute la période de simulation. Le Tableau 1 présente le nombre p d'éléments ajoutés pendant chaque pas de temps ainsi que le nombre final d'éléments et le temps de calcul pour chaque simulation.

La Fig. 1 présente les profils de concentration obtenus pour les trois simulations. La solution de référence, obtenue avec la méthode des éléments finis standard, avec une discrétisation fine en espace et en temps (Δx=0,05 cm, Δt=0,1 s), nécessite un temps de calcul de 43,5 s. Ces résultats montrent la robustesse de la méthode développée. Même avec un seul pas de temps de Δt=1000 s et un maillage de seulement 26 éléments, les résultats obtenus restent précis et ne présentent, ni oscillations, ni diffusion numérique. En effet, la méthode ELLAM permet l'utilisation de grands pas de temps, vu que la variation de la concentration est beaucoup moins importante le long des caractéristiques que la variation dans le temps à une position fixée. L'utilisation d'un seul pas de temps peut entraı̂ner, cependant, une sous-estimation de la dispersion [3].

4.2 Transport multi-espèces

Le modèle proposé est ensuite utilisé pour simuler le transport multi-espèces (ammonium, nitrite et nitrate) avec sorption et dégradation d'ordre 1, dans un écoulement monodimensionnel. On considère le cas où l'espèce 1 est adsorbée à la surface des grains et se dégrade pour donner l'espèce 2, qui elle-même se dégrade pour donner l'espèce 3. La solution analytique à ce problème a été développée dans [5].

Le système d'équations à résoudre est le suivant : (7)

Les conditions initiales et aux limites sont :

C_{1} (0, t) = 1, C_{2} (0, t) = 0, C_{3} (0, t) = 0 C_{1} (z, 0) = C_{2} (z, 0) = C_{3} (z, 0) = 0

Les valeurs des paramètres sont données dans le Tableau 2.

Vitesse [ $cm h^{- 1}$ ]	u	1,0
Coefficient de dispersion [ ${cm}^{2} h^{- 1}$ ]	D	0,2
Facteur de retard	R ₁	1,0
Coefficient de réaction de l'espèce 1 [h⁻¹]	k ₁	0,01
Coefficient de réaction de l'espèce 2 [h⁻¹]	k ₂	0,1
Coefficient de réaction de l'espèce 3 [h⁻¹]	k ₃	0,0

Le domaine représente une colonne de 60 cm de long. La discrétisation spatiale est uniforme avec un pas d'espace Δx=1 cm. Le pas de temps fixe est Δt=5 h. La Fig. 2 présente les profils de concentration pour chaque espèce pour une durée de simulation de 50 h. Les résultats numériques obtenus avec les ELLAM sont en très bonne adéquation avec la solution analytique donnée dans [5].

5 Conclusion

Le transport de contaminants avec des processus biologiques, chimiques ou radioactifs en milieu poreux peut présenter des fronts de concentrations raides, qui nécessitent des méthodes numériques adaptées. La nouvelle méthode ELLAM, où les fonctions test ainsi que la solution sont définies sur des maillages variables, reproduit ces phénomènes sans oscillations ni diffusion numérique. Les expériences numériques effectuées ont permis de mettre en valeur la précision et la robustesse de cette méthode pour la résolution des équations de transport réactif en milieux poreux non saturés.

References

[1] J.-P. Benque; J. Ronat Quelques Difficultés des Modèles Numériques en Hydraulique (R. Glowinski; J.-L. Lions, eds.), Computing Methods in Applied Sciences and Engineering, North-Holland, Amsterdam, 1982, pp. 471-494

[2] P. Binning, Modeling unsaturated zone flow and contaminant in the air and water phases, Ph.D. Thesis, Department of Civil Engineering and Operational Research, Princeton University, Princeton, NJ, 1994

[3] P. Binning; M.A. Celia A forward particle tracking Eulerian–Lagrangian localized adjoint method for solution of the contaminant transport equation in three dimensions, Adv. Water Resour., Volume 25 (2002), pp. 147-157

[4] M.A. Celia; T.F. Russell; I. Herrera; R.E. Ewing An Eulerian–Lagrangian localized adjoint method for the advection–diffusion equation, Adv. Water Resour., Volume 13 (1990), pp. 187-206

[5] C.M. Cho Convective transport of ammonium with nitrification in soil, Canad. J. Soil Sci., Volume 51 (1970), pp. 339-350

[6] E.B. Diaw; F. Lehmann; P. Ackerer Modélisation du transport d'un soluté réactif en milieux poreux non saturés, C. R. Acad. Sci. Ser. IIa, Volume 333 (2001), pp. 129-132

[7] J. Douglas; T.F. Russell Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures, SIAM J. Numer. Anal., Volume 19 (1982) no. 5, pp. 871-885

[8] R.W. Healy; T.F. Russell A finite-volume Eulerian–Lagrangian localized adjoint method for solution of the advection–dispersion equation, Water Resour. Res., Volume 29 (1993), pp. 2399-2413

[9] A. Oliveira; A.M. Baptista On the role of tracking on Eulerian–Lagrangian solutions of the transport equation, Adv. Water Resour., Volume 21 (1998), pp. 539-554

[10] T.F. Russell; R.V. Trujillo Eulerian–Lagrangian localized adjoint methods with variable coefficients in multiple dimensions, Computational Methods in Surface Hydrology, Proc. 8th Int. Conf. on Computational Methods in Water Resources, Venice, Italy, 1990, pp. 357-363

[11] J. Touma; M. Vauclin Experimental and numerical analysis of two-phase infiltration in a partially saturated soil, Transp. Porous Media, Volume 1 (1986), pp. 27-55

[12] M.T. Van Genuchten A closed-form equation for predicting the hydraulic conductivity in soils, Soil Sci. Soc. Amer. J., Volume 44 (1980), pp. 892-898

[13] H. Wang; R.E. Ewing; M.A. Celia Eulerian–Lagrangian localized adjoint method for reactive transport with biodegradation, Numer. Methods Partial Differential Equations, Volume 11 (1995), pp. 229-254

[14] A. Younes, An accurate moving grid Eulerian–Lagrangian Localized Adjoint Method for solving the one-dimensional variable-coefficient ADE, Int. J. Numer. Meth. Fluids, submitted for publication

Comments - Policy