Les équations de transport de masse

L’équation de conservation s’écrit :

[flux massique entrant] - [flux massique sortant] =
modification du stock massique au cours du temps (en unités de
masse/unité de temps).
δcω (1)
 div(J) 
δt
avec  = porosité
A. LA DIFFUSION (FICK)
L’équation de diffusivité définie par FICK pour la phase liquide seule:
J   D.grad(C) (2)
D° : coefficient de diffusion moléculaire

Quelques coefficients de diffusion moléculaire :

CATIONS D° en 10-6 cm²/s ANIONS D° en 10-6 cm²/s
H+ 93,1 OH- 52,7
Na+ 13,3 F- 14,6
K+ 19,6 Cl- 20,3
Rb+ 20,6 Br- 20,1
Cs+ 20,7 HS- 17,3
Mg+2 7,05 HCO3- 11,8
Ca+2 7,93 CO32- 9,55
Sr+2 7,94 SO4 2- 10,7
Ba+2 8,48
Ra+2 8,89
Mn+2 6,88
Fe+2 7,19
Cr+3 5,94
Fe+3 6,07

Dans un milieu poreux, la phase liquide baigne un milieu poreux

complexe d’où une loi légèrement modifiée.


J  D. grad(C.ω)  τ
v
 (2bis)

avec V : le volume moyen de l’échantillon poreux

 : le vecteur de tortuosité
 : la porosité
Le passage du coefficient de diffusion dans un milieu biphasique où
la tortuosité est mal connue entraînera la définition d’un
coefficient de diffusion effective D’ avec :
2
 DD' D. 
2  (2 -  ) (3)
on donne aussi :
D'  μ .ω D
(4)
(avec 1)
et
D'     (5)

(avec  la tortuosité : rapport du trajet parcouru / longueur directe)

B. LA DISPERSION MECANIQUE
Suite à l’injection d’une masse de polluant sur un créneau de temps très court, on voit (sur
des piézomètres implantés sur l’écoulement) en fonction de la distance une modification des
concentrations.
C/Co
t= t1
C/ Co C/ Co
1
Co =1 =0
t= t2
C/ Co 0 ,5
Co =0
t= t3
C/ Co
Co =0 0
t
Z o n e d e d is p e r s io n début de p o s it io n d u
l'e s s a i fr o n t d e c o n v e c t io n
t0 t1 t2
1 D ir e c t io n
C/Co d e tra n s p o rt

L d is t a n c e L
t1 t2
p o in t d 'in je c t io n

On constate que les concentrations se distribuent selon la loi normale. Si L² est la
variance de distribution des concentrations pour un écoulement mono-dimensionnel
(L²/2t) = constante = DL coefficient de dispersion longitudinale dans la seule direction.
Si la vitesse est constante, alors
2
DL   L .v (6)
2x
x étant la distance au point d’injection.
Dans le cas de l’écoulement bidimensionnel :
2
DT   T (7)
2t
DT= coeff. de dispersion transversale
On voit que le panache de pollution en 2D et un écoulement unidirectionnel aura une forme
elliptique en ce qui concerne la distribution des concentrations est statistiquement normale
selon les directions longitudinale et transversale.
Dans le domaine 3D, le panache de pollution sera un ellipsoïde de révolution si les 2
composantes transversales de la dispersion sont identiques et nettement inférieures à la
dispersion longitudinale.
On peut interchanger la statistique espace/temps par :

D L 2  v 2 .σ t 2 2 2
et D L  v . T (8)
2t
Le coefficient de dispersion longitudinale peut ainsi être estimé à partir d’une courbe de
restitution d’un traceur.
PROBLEMES

1. Soit 2 milieux différents poreux et fracturés :

-
un milieu fracturé avec porosité 3 10-4
-
un milieu poreux avec porosité 0,3
en considérant que chacun des milieux possède un coefficient de perméabilité K = 10-8 m/s.
Déterminer en convection pure le temps nécessaire à un traceur pour parcourir une
distance de 30m dans chacun de ces milieux avec un gradient de 0,01.

2. Soit un écoulement entrant sur la face gauche d’un cube de grès de 1 m de coté avec une
vitesse horizontale de 10-5 m/s. La porosité  = 0,10 et la diffusion effective du traceur qui
entre de 10-10 m2/s. Si le traceur possède des concentrations de 120, 100 et 80 mg/l
respectivement sur la face gauche au centre du cube et à la sortie sur la face droite,
calculer le flux massique passant sur le plan central du cube du à la convection seule
et à la diffusion seule. J en kg/m2.s

1 m

3. Une expérimentation sur une colonne de milieu poreux à l’aide d’un traceur permet de
mesurer la dispersivité à partir d’une courbe de restitution. On obtient une variance estimée
de la concentration en fonction du temps. Soit une colonne de matériaux poreux d’une
longueur de 56 cm. La courbe de restitution montre que le front de convection arrive à 82
minutes à la sortie de la colonne avec t2 = 5,5 min2. Calculer la dispersivité longitudinale
du matériau L (=DL/v).

F r o n t d e c o n v e c t io n
1
C/ Co

56 cm

8 2 m in te m p s

Les approches analytiques

La solution 1D de l’équation de convection dispersion est donnée par
:
 2C C C C(0,t) = C0 C(x,0) (13)
Dx .  vx .  =0
2 xx t
1ère condition ; à x = 0 pour toute valeur de t , C = C0 = constante
2ème condition ; à t = 0 en tout x, C = 0

La concentration au départ de l’expérience sur colonne de matériaux poreux est égale à 0,

d’où :
 C   x  v.t  (14)
C(x, t)   0 erfc 1/2 
ε
 2   2(α x .v.t) 

avec
x : la dispersivité longitudinale
erfc : la fonction complémentaire d’erreur
v : vitesse linéaire de l’eau = vitesse du traceur introduit
 : terme qui a été négligé (pour Pe= (v.x/D) > 500, l’erreur est de 3%.
La quantité (x.v.t)1/2 est remplacée par (D.t)1/2 car D L  α L .v
avec D, le coefficient de la dispersion hydrodynamique.
Si la diffusion est petite par rapport à la dispersion mécanique, le coefficient de dispersion #
x.v, la solution devient :
  x  v.t   x  v.t   (15)
C   v.x 
C(x, t)   0 .erfc    exp .erfc  
 2   2(D.t)

1/2 
  D   
 2 D.t  1/2 

 

Si le nombre de Péclet croit très vite,
on ajoute cette partie

erfc = 1 - erf(B), erfc(-B) = 1+erf(B), eff(-B) = -erf(B)

pour :
B erf(B) erfc(B)
0 0 1
0,5 0,52 0,479
2 0,84 0,157
1,5 0,966 0,033
2 0,995 0,0046
2,5 0,999 0,0004
3 0,9999 0,000002
F o n c t io n
B = 2 (D .t) 1 / 2
2
1 ,5 e rf( B )
1 1
5
0 0
-0 ,5 e rfc ( B )
-1 -3 0 3
B
(avec erfc ; fonction complémentaire d’erreur)
La courbe de restitution du traceur à la sortie de la colonne peut se mettre soit en fonction
du temps (A), soit pour un temps donné en fonction de la distance le long de la colonne.
Lorsque le rapport C/Co atteint la valeur de 0,5, on considère que le front de convection s’est
manifesté.
En ce qui concerne les fonctions erf et erfc, on fera les remarques suivants :
* erf(B) a une valeur comprise entre -1 et +1
* erfc(B) a une valeur comprise entre 0 et +2, mais elle ne dépasse 1 que pour des
arguments B négatifs.
Sur le plan physique, trois positions sont intéressantes à connaître :
- sur le front de convection, C = Co/2
- devant le front de convection, C = 0
- en arrière du front de convection, C = Co
* Sur le front de convection x = v.t, l’argument B = 0 et erfc(0) = 1, aussi C/Co = 0,5
* Devant le front de convection quand x >> v.t, erfc(B) = 0 et il n’y a pas de mélange. On
observe les concentrations originelles du fluide. Pour des raisons pratiques, si B>2 alors
C=0.
* A l’arrière du front de convection, B<0 et erfc(-) = 2, ce qui signifie que C=Co -
Pratiquement dès que B approche la valeur -2, erfc(-2) # 2.
Remarque : la solution analytique d’Ogata-Banks n’est utilisable qu’au laboratoire, car
quelque soit la distance entre le point d’observation et le point d’injection continue, on obtient
une concentration maximale équivalente à celle de l’injection. Cette condition n’est pas
réaliste sur le terrain, car apparaissent des phénomènes de retard liés à l’adsorption.

5. Un traceur est injecté dans une colonne de 0,30m de long avec une vitesse uniforme de
10-2 cm/s. les rapports C/Co de 0,42 et de 0,573 sont notés à 46,6 et 53,3 minutes après le
début du test. Quelle est la valeur la dispersivité longitudinale ?
On utilisera ici :
1  x  vt 
C  erfc 
 2 α x .v.t  
Co 2

A- L’équation du retard Rf
Dans le cas d’une sorption linéaire on pose :
D x . 2 C v x .C C (15bis)
 
R f .x 2 R f .x t

ou, 1 ω  (15ter)
Rf  1   ρs.Kd
 ω 
avec Kd = le coefficient de distribution d’après l’équation linéaire isotherme de Freundlich.

S C
 Kd.
t t
s = densité du solide (=2,65 g/m3 pour le sable)

L’équation d’Ogata-Banks s’écrit alors :

C0  R f .x  v. .t  (16)
C erfc  
C 
 2  x .v. .t.R f
1/ 2
 
Pour Rf = 1, on retrouve l’équation simplifiée; notons ici que V représente la vitesse de
l’eau pure.
La vitesse du contaminant introduit est inférieure à celle de l’eau
souterraine. Mathématiquement, on pose que :
(17)
vω vω
vc  
Rf 1   1  ω  /ω.s.Kd
C’est l’équation du retard qui va freiner l’arrivée du front de convection en sortie de colonne.
Si Kd = 0, alors Vc = V, car Rf = 1
Kd est exprimé souvent en ml/g et s en g/cm3.
Le coefficient de distribution Kd peut être également défini par :
Ks.CEC (18)
Kd 
τ
Avec
Ks : coefficient de sélectivité
CEC : capacité d’échange cationique (meq/g)
 : la concentration des cations en compétition dans la solution (meq/g)
On peut écrire aussi que :
vω ρs 1 - ω .Ks.CEC (19)
 1
vc ω.τ
Cette équation est très utile pour étudier la contamination des nappes par des composés
organiques. On peut sur un site expliquer les différences entre les distributions observées de
certains contaminants avec la prévision du mode de propagation.

6. Soient des dépôts fluvio-glaciaires recouvrant des granites. La partie inférieure est
sableuse et le set argileux. A quelle vitesse circulera un traceur radioactif Sr 90  ? on
donne la vitesse de filtration dans les sables = 1,85 c m/j, s (1-) = 1,9 g/cm3, CEC = 1,4
meq/100g,  = 4,8 meq/l,  = 0,35 et Ks = 1,3.

B- La décroissance radioactive, la biodégradation et l’hydrolyse

La solution 1D est donnée par Bear (1979) :

C0  x    1  4 1/2    x  vt1  4α x /v  1/2 

C exp .1   x
  .erfc   (20)
2  2α x    v     2α x .vt  1/2 
où  représente la constante de désintégration.
 = 0,693/demi vie (R:e0=1 , e-=0) avec Ln(2) = 0,693
Pour  = 0, on retrouve l’équation d’Ogata-Banks.
Lorsque le terme (1+4x )/v s’accroît avec , le terme exponentiel approche de 0 comme
la concentration. Le matériel est réactif et décroît plus vite que le transport. Si la vitesse
devient importante, le terme de l’exponentielle précédent approche de 0 et on s’approche
alors de la solution d’Ogata-Banks. Un transport de masse plus rapide que la dégradation
(la désintégration) rend cette dernière inopérante.
La compétition existe aussi entre le transport et la cinétique en particulier entre la
concentration des espèces chimiques et leur vitesse d’écoulement. Plus la dispersivité
longitudinale x devient forte, plus (1+4x )/v devient grand et la décroissance domine.
Par exemple, pour un nucléide de demi-vie de 1000 ans, une vitesse de 4 m/an, une
distance de 5000 m, si x=30, 100 puis 300 m, le terme exponentiel est égal respectivement
à 0,5, 0,505 et 0,510.

Tableau : Valeurs expérimentales du coefficient de distribution Kd en ml/g à 25°C et 0,1 Mpa (obtenues par
des essais en batch).
Eléments Sels Basaltes Tufs Granites
Tc 2 20 à 0 10 à 0 4
Pu 500 à 50 200 à 100 500 à 40 500 à 100
Np 30 à 7 50 à 3 50 à 3 50 à 1
I 0 0 0 0
U 1 6 4 4
Cs 1 à 800 300 100 300
Ra 5 50 200 50
sr 5 100 100 12
C 0 0 0 0
Am 300 50 50 200
Sn 1 à 50 10 à 100 50 à 500 10 à 500
Ni 6 50 50 10
Se 20 à 200 20 à 5 2 2
Cm 300 50 50 200
Zr 500 500 500 500
Sm 50 50 50 100
Pd 3 50 50 10
Th 50 à 100 500 500 500
Nb 50 100 100 100
Eu 50 50 50 100
Pa 50 100 100 100
Pb 2 25 25 5
Mo 0à5 10 à 4 10 à 4 5à1

L’équation reflète le régime permanent lorsque l’argument de erfc approche la valeur de -2.
L’équation s’écrit alors :
  1  (21)
 x    1  4λλ x  2  
C  C 0 . exp   1   
  
 2α x   v 
  en effet erfc (-2)=2
  
Cela n’est possible que si le traceur est sans cesse renouvelé à l’entrée du système.

Dans l’équation d’Ogata-Banks, la concentration en régime permanent est indépendante de

la dispersivité longitudinale, alors que dans l’équation ci-dessus la concentration dépend de
la convection et de la décroissance.

Le terme v représente la vitesse du contaminant qui peut être égale à la vitesse linéaire de
l’eau souterraine si le contaminant n’est pas retardé. Pour un contaminant qui présente un
certain retard, il faut utiliser V/Rf.

C- La dispersion transversale
Il faut parfois faire intervenir la géométrie de la zone d’injection du polluant: un point, une ligne verticale ou une
surface plane verticale.
y y y
Y
x x x
Z
z
z

C0Q C0Q C max  C 0

C max  C max 

4xπx α yα z 2 1 
2v πα y x 2 1
Q étant le débit en m3/s Q en m²/s C0Q en kg/m.s ou en kg/s

Ces équations démontrent l’importance de la dispersivité transversale :

 pour la ligne verticale d’émission du traceur, plus y est élevé, plus l’extension du
panache de pollution sera étendue latéralement et plus la concentration sera faible.
 pour le point source, les résultats dépendent de la dispersivité latérale y et verticale z.
Les concentrations sont pour ces équations, indépendantes du temps et approchant la
valeur infinie si x tend vers 0. Ces équations sont seulement utiles pour des distances
importantes (x très grand).
La surface plane est le plus pratique avec une dispersion dans les 2 directions y et z, c’est-à-
dire latéralement et verticalement. Développons l’expression (Domenico et Palciaukas,
1982):
-
pour la composante latérale y de la dispersion :
    
C0  y  Y/2   y  Y/2  
C erf     (22)
2 
   
 2α x .vt  1/2   2 α y .x 1/2  
 
où la dimension Y/2 de la source plane figure dans la solution.
-
pour la composante latérale z verticale :
C    z  Z   z  Z  
C  0 erf    (23)
2   2α .x 1/2   2α .x  1/2  
  z   z 
où la taille Z de la source plane intervient.
La dispersion dans les directions y et z est représentée par le produit des 2 relations
précédentes (22) et (23).
Dans le plan de symétrie pour y = z = 0, on obtient :
 Y   Z 
Cmax  C 0 erf  .erf   (24)
 4α c .x  1/2   2α z .x  1/2 

Les dimensions de la source et les dispersivités transversales contrôlent les concentrations

maximales qui peuvent se manifester en régime permanent dans un panache de pollution.

Par exemple, si Y ou Z augmentent en valeur, ou si y et z diminuent, les augmentations

de la fonction erreur s’approchent de +2 et la valeur des fonctions d’erreur de la valeur 1-
Dans ce cas la dispersion est restreinte et la concentration est maintenue autour de Co pour
les valeurs faibles de x.

D- Les modèles de transport multidimensionnels

La solution analytique inclue 3 composantes de dispersion, une vitesse de convection
constante et un terme cinétique :
 2C  2C  2C C r C
Dx  Dy  Dz - Vx. - 
x 2 y 2 z 2 x ω t (25)
r étant défini par une loi mathématique.
Cette équation permet de traiter les entrées instantanées ou continues.

a- Les entrées continues de traceur

Considérons 4 cas :
s o u rce s o u rc e
C 0 C 0

x = v.t x = v.t
C a s a : é c o u le m e n t p is t o n C a s b : d is p e r s io n lo n g it u d in a le
s e u le
x = v.t x = v.t

C a s c : d is p e r s io n t r a n s v e r s a le C a s d : d is p e r s io n lo n g it u d in a le
s e u le e t t r a n s v e r s a le
Dans le cas de la dispersion transversale, les concentrations en régime permanent sont
toujours inférieures à Co, sauf en x = 0.
Pour le cas d’une source plane verticale Y x Z on pose :

C    x  v.t    

 y  Y/2   erf   y  Y/2  .
C x, y, z  t   0 .erf  .erf  
 8 
 
 
 2α x .v.t 1/2    2 α y .x 1/2 
   y 
 2 α .x 1/2  

(26)

   z  Z/2    z  Z/2   
erf    erf  
  2α z .x   2α z .x 
1/2  1/2  
 

Pour x  (H-Z)²/ z et une dispersion verticale vers le haut et le bas (Z/2).

Pour x << v.t, l’argument de erfc approche de -2 et l’on retrouve le produit des relations (22)
et (23)
Pour le plan de symétrie (y = z = 0), on obtient alors :

   
 C    x - v.t     Y  Z 
C 0,0, t    0 erf  .erf .erf   (27)
 2   2α x .v.t  1/2
 

   4 α y .x 1/2 
   2α z .x  1/2  


Remarque : 2* pour 1 direction z, 4 pour 2 direction z.

On peut montrer que si y et z tendent vers zéro, on retrouve l’équation d’Ogata-Banks
(14), car erf () = 1.
Si la dispersion z est négligée, alors les fonctions erf contenant le terme z sont ignorées et
Co/8 devient égal à Co/4.
Si l’aquifère est de faible épaisseur, le panache de pollution va occuper tout l’espace.
la distance x’ sur laquelle la dispersion se développe est approximée par :

x'   H  Z 2 /α z (28)

x'
z
Pour des distances > x’, la distance x au
dénominateur de la fonction erreur en z
H est remplacée par x’.
Toute dispersion selon z est absente
pour x > x’.
 (0 ,0 ,0 ) Y/ 2 a : s o u r c e p o n c t u e lle s a n s
Y/ 2 d is p e r s io n v e r t ic a le d a n s
Z le p la n h o r iz o n t a l

c : d is p e r s io n à p a r t ir
d 'u n e s o u r c e lin é a ir e
v e r t ic a le
(0 ,0 ,0 ) Y/ 2
Y/ 2

b : d e u x d ir e c t io n s d e
d is p e r s io n s e lo n z
s o u r c e d e t y p e s u r fa c e
v e r t ic a le p la n e
Y/ 2
Y/ 2
(0 ,0 ,0 )
Z/2
Z/2

z
y

x
7. Des déchets de diéthyle éther et de tétrachlorure de carbone ont été enterrés dans un
sable aquifère depuis 15 ans. Calculer les concentrations de chaque contaminant le
long d’un plan de symétrie du panache de pollution au point x = 225, y = 0 , z = 0 au
temps de 15 ans (= 4,73 10 8 s) ou 5475 h. la vitesse de l’eau est estimée à 10 -6 m/s. Le
facteur de retard respectif est de 1,5 et 27,4. Les concentrations à la source de
contamination sont successivement de 1.104 g/l et 5.102 g/l. La taille de la source de
contamination est de Y = 25 m et Z = 5 m. Les dispersivités sont estimées à X = 1m , Y =
0,1m et Z = 0,01m. On suppose que le panache de pollution se disperse selon les deux
directions z et les deux directions y. on fait appel ici à solution :
   
  x  vt     
C0  Y Z 
C  x,0,0, t   erfc .erf  .erf  
2  2  
  x .v.t 
4

 
 y .v.t 
 4
  z .v.t   

 
   


On peut introduire la dégradation ou la décroissance dans l’équation de Domenico (1987)

suivante :
 x    1  4λλ x 1/2    x - vt 1  4λλ /v 1/2 
C0  
C x, y, t   exp   1     .erf  x .
8 2α x    v     2α x vt 1/2 

 
(29)
    
  y  Y/2   y  Y/2      zZ   z  Z  
erf  erf .erf    erf  

 
 
2 α yx 1/2

 
2 α yx
1/2 


  2α z x  
1/2
 2α z x   
1/2

Si  = 0, la concentration le long du plan de symétrie est calculée par (26) :
 x   1/2 
C0   1  4λλ x     x - vt1  4λλ /v  1/2 
C x, y, z, t   exp  1     .erfc  x .
2 
 2α x  
  v  
 
 2 α x vt  1/2 
(30)
 x - vt1  4λλαx/  1/2   
   
erfc  .erf 
Y .erf  Z 

 2α x vt  1/2  
 
 
4 α yx 1/2 
  
2 α yx 
1/2 

Si x et y s’approchent de zéro, on retrouve l’équation de Bear (1979), développée
antérieurement (20).
Une relation en régime permanent peut être obtenue si l’argument de la fonction
complémentaire d’erreur est de l’ordre de -2.
Il faut éventuellement corriger la vitesse par le coefficient de retard Rf.

b- Les entrées instantanées

Les deux grands types de modèles sont basés sur:

- une source finie parallélépipède (Hunt, 1978),
- un point source (Beatsle, 1969).
Dans ce dernier cas, le modèle point source est le suivant :
 
 exp   x  vt   y
C 0 .V0  2 2 z2 
C x, y, z, t      λt 
 
 8πt  3/2 . D D D 1/2 
x y z  

4D x t 4D y t 4D z t

avec
C0 : concentration originelle, V0 = volume de départ injecté
DxDyDz : représentent les dispersions (coeff.) hydrodynamiques
v : la vitesse du contaminant  : la constante de désintégration
Si le produit n’est pas radioactif, t est négligé.
 D ir e c t io n d 'é c o u le m e n t
p o in t d 'in je c t io n
( v o lu m e V 0 )

x= 0 x= vt1 x= vt2 x= vt3

Le pic de la concentration maximale apparaît au centre d’un nuage avec y = z = 0 et x = vt.
C 0 .V0 .e λ t
Cmax 
 
8πt  3/2 . D x D y D z 1/2 (32)

Les dimensions du nuage de polluant en les supposant issues d’un point sont données par :

 
3σ x   2D x .t  1/2 , 3σ y  2D y .t 1/2 et 3σ z   2D z .t  1 / 2

où  est l’écart type, donc 3 représentent 99,7% de la masse.

8. Un rejet accidentel de déchets radioactifs d’une activité totale de 1000 ci comprend 400 ci
de Cs137, 400 ci de Sr90, 100 ci Pm147 et 100 ci H3 avec des demi-vies respectives de 33 ans,
28 ans, 2,7 ans et 12,26 ans. Une rivière existe à 100 m à l’aval nappe du point de rejet
dans l’aquifère. Déterminer les temps mis par les contaminants pour atteindre la rivière
et les concentrations maximales atteintes. La vitesse de l’eau est de 10 cm/jour (1,16 10-4
cm/s).
Radioéléments Kd (ml/mg) Rf Dx (cm2/s) Dy = Dz
Cs137 10 47,6 10-4 10-3
Sr90 0,6 2,885 10-3 10-5
Pm147 100 476 10-5 10-5
H3 0 1 10-5 10-5

9. Un contaminant organique est déposé dans une tranchée pénétrant dans une nappe peu
profonde semi-infinie. Calculer la concentration du contaminant (qui est sans cesse
renouvelé) après 3 ans et à 10 m de distance de la tranchée continue.

x = 10m

C/ Co R f = 3 ,0
x = 4 ,3 m
=
100m g/ l
v = 2 ,3 1 1 0 -6 m / s

La programmation
Des programmes en Q BASIC peuvent être aisés à écrire pour les solutions
analytiques 1D, 2D et 3D avec le coefficient de retard, la décroissance en régime
permanent ou en régime transitoire.
2   x
erf(x)   e .dx
π 

10. Dans une zone polluée vous remarquez que les chlorures ont migré de 1000 m alors que
le cadmium ne s’est délacé que de 100 m. En supposant que ces deux éléments chimiques
aient été disséminés en même temps trouvez le coefficient de distribution Kd si a
porosité de l’aquifère est de 0,2 et sa densité s = 2,65.

11. Un épandage accidentel est constitué des radioéléments

suivants :
Éléments Activité(Ci) T1/2 (ans) Kd(ml/g) Dx(cm2/s) Dy Dz
A 400 33 8 10-4 10-5 10-6
B 400 28 52 10-3 10-5 10-7
C 400 2,7 40 10-5 10-5 10-5
D 400 20 0,1 10-6 10-5 10-6

 Quel est le constituant qui se déplace le plus rapidement puis le plus long ?
 Quel constituant montrera la concentration la plus faible puis la plus forte ?
 Quel constituant donnera un panache le plus petit autour du point de
contamination dans le sens d’écoulement de la nappe… puis le panache le plus
grand ?
 Quel constituant donnera un panache de pollution identique dans toutes les
directions ?

12. Une cellule de stockage de déchets dangereux a été construite dans une zone donnée.
Une zone tampon est nécessaire entre les tranchées de stockage et une limite donnée. Plus
de 200 constituants sont stockés avec des coefficients de distribution très variables.
Quelques années plus tard un piézomètre construit sur la limite de la cellule de stockage
montre une pollution de 6 m d’épaisseur et étendue sur une longueur de 60 m. La faible
perméabilité des matériaux ne permet pas de dispersion verticale.
Calculer la concentration maximum qui peut se manifester (en terme de Co) une fois
que le panache arrive sur le bord de la zone tampon qui est à 150 m de la limite du
stockage des déchets. Le coefficient de dispersion transversale est de 10 -3 m2/s et la
vitesse de 10-3 m/s.

13. Soit un écoulement des eaux souterraines à une vitesse linéaire de 15 cm/j. Estimer la
vitesse d’un contaminant organique ayant un coefficient de distribution Kd = 6,6 ml/g
dans un milieu poreux de porosité 0,35 et s = 2,65 g/cm3.
14. Un panache de pollution de bromoforme est issu d’une lagune et montre une longueur de
150 m. La lagune pénètre l’aquifère sur une distance de 5 m alors que l’épaisseur est de 15
m au total. Evaluer la distribution des concentrations dans la nappe au bout de 14 ans
avec :
la vitesse linéaire de l’eau souterraine, v = 1,6 10-3 cm/s ; les dispersivités X = 5m , Y = 0,5m
et Z = 0,1m ; Co = 4000 mg/l ; foc = 0,0045 ; log (Ko) = 2,38 ; porosité = 0,30 ; s = 2,65
g/cm3.
On ferra appel à la relation de Hasset (1983).
 Résolutions numériques de l’équation de transport

Il existe de nombreuses méthodes numériques pour résoudre l’équation de transfert de

masse. Citons par exemples : les différences finies, les éléments finis, les caractéristiques et
la marche au hasard. Nous allons décrire les aspects théoriques, les avantages et les
inconvénients de ces méthodes.

1. méthode des différences finies

la résolution numérique des équations aux dérivées partielles s’est considérablement
développée depuis une cinquantaine d’années avec la mise au point d’ordinateurs de plus
en plus rapides et ayant une taille mémoire de plus en plus importante.
Ces méthodes numériques s’appliquent pour les prévisions de l’évolution de nombreux
phénomènes naturels spécialement dans le domaine des sciences de l’eau. Cette méthode
sera présentée dans le cas d’une configuration bidimensionnelle et radiale (discrétisation,
différents schémas, traitement des conditions aux limites ainsi qu’une présentation de
l’analyse de la stabilité).

1.a) Ecoulement bidimensionnel

Dans le cas d’un écoulement bidimensionnel, on suppose que la vitesse de l’écoulement est
constante. Dans ce cas on injecte uniformément le traceur en un point d’une nappe
homogène et d’épaisseur constante. Le diamètre du puits d’injection est faible par rapport à
la distance horizontale. Il est nécessaire que la distance du point d’injection au point de
prélèvement soit grande en comparaison avec la distance nécessaire pour
l’homogénéisation des concentrations par dispersion transversale.
L’équation aux dérivées partielles de transfert de masse en mode convection-dispersion,
quand l’écoulement est bidimensionnel uniforme et le milieu est homogène, s’écrit :
C  2C  2C C
 DL  DT u
t x 2 y 2 x
avec ( D L  α L .u ) et ( D T  α T .u ) en supposant que les coefficients u, DL et DT sont
constants. La forme discrétisée de l’équation peut s’écrire sous la forme suivante :

C i, j t Δt  C i, j t D  D
 L DD L  T DD T 
u
CC
Δt Δx 2 Δy 2 Δx

avec
  
DD L  δ C i 1 , j 2Ci, j  C i 1, j t Δ t  1  δ  C i 1, j  2Ci, j  C i 1, j t
DD T  δ C i, j1  2C i, j  C i, j1  t Δt  1  δ  C i, j1  2Ci, j  C i, j1  t

    
CC  1  λ  γC i 1, j t Δt  1  γ  C i 1, j t   2λ  1 γC i, j t Δt  1  γ  C i, j t  λ γC i 1, j t Δt  1  γ  Ci 1, j t 
où  = coefficient de pondération temporelle pour le terme dispersif ( 0  δ  1 )
 = coefficient de pondération temporelle pour le terme dispersif ( 0  γ  1 )
 = coefficient de pondération temporelle pour le terme dispersif ( 0  λ  1 )
On obtient un schéma implicite pour ( δ  γ  1 ), un schéma explicite pour ( δ  γ  0 ) et un
schéma centré en temps pour ( δ  γ  0,5 ). De la même manière on obtient des schémas
décentrés en aval ( λ  1 ), des schémas décentrés en amont ( λ  0 ), des schémas centrés
en espace ( λ  0,5 ).
La discrétisation bidimensionnelle par la méthode des différences finies dans le domaine
concerné aux intervalles de temps et d’espace peut être réalisée par des cellules
rectangulaires en supposant des dimensions constantes ( Δx ) et ( Δ y ) respectivement en x
et y.
 X
( i, j+ 1 )

( i-1 , j) n o e u d i, j
Y ( i+ 1 , j)

E p a is s e u r d e
l'a q u ifè r e
( i, j-1 )

z
y

Fig. ???- discrétisation bidimensionnelle d’un aquifère.

Le centre des cellules est appelé nœud. Les indices (i,j) désignent la position d’un nœud. Le
bilan des échanges est déterminé entre le centre de la cellule (i+1,j) et le centre des cellules
voisines aux intervalles de temps (t,t+t), en prenant en compte la convection, la dispersion,
la diffusion et éventuellement les phénomènes physico-chimiques.
La détermination de la masse d’entrée dans une cellule par convection s’obtient en fonction
de la vitesse d’écoulement et du gradient de concentration. La concentration peut être
calculée par une moyenne pondérée entre les cellules voisines. Deux méthodes de choix
peuvent être appliquées ; si le transfert convectif est prépondérant, le coefficient de
pondération doit être choisi sur le coté amont de la cellule et si la dispersion est dominante
une pondération symétrique est préférable.

Analyse de la stabilité
Au regard du choix des différents schémas de type différences finies, il existe plusieurs
critères pour le choix du pas de temps. L’un de ces critères est le critère de courant qui
exprime que la masse d’entrée dans une cellule par convection pendant un intervalle de
temps (t) ne peut pas être supérieure à la masse de sortie (Kinzelbach, 1986) :
Δt.u x Δt.u y
C0 x   1 et C0 y  1
Δx Δy
ou bien C0 x  C0 y  1
ou C0 est le nombre de courant.

1.a) Dispersion numérique

La plupart des solutions de différences finies utilisent le schéma centré en espace et en
temps. Néanmoins, ce schéma crée des oscillations. Pour les éviter, on peut choisir une
discétisation plus fine. Le schéma décentré en amont est sans oscillations, mais engendre
une dispersion artificielle.
Le phénomène de dispersion numérique est observé dans le cas d’un transfert de masse
purement convectif. Si on considère le cas d’un écoulement parallèle à l’axe x, le front de
concentration commence par le nœud (i,j) et progresse dans un intervalle de temps
(tx/ux) pour se trouver entre les nœuds (i,j) et (i+1,j). La masse d’un contaminant est
propagée entre les deux nœuds à cause de la discétisation, c’est à dire que le front n’a pas
encore atteint le nœud (i+1,j) mais ce nœud indique déjà une concentration supplémentaire.
En outre si l’écoulement n’est pas parallèle à l’axe x, il y aura une dispersion numérique
angulaire, c’est-à-dire que le front va s’étendre de la direction 02 du nœud (i,j)
jusqu’aux nœuds (i+1,j) et (i,j+1) à cause du mouvement convectif du contaminant.
Une autre méthode pour analyser la dispersion numérique consiste à considérer ce
phénomène comme une erreur de troncature (Holly, 1975). Si on suppose que la direction
d’écoulement est parallèle à l’axe x, la discrétisation des dérivées premières en espace et en
temps ( C C
x , t ) de l’équation de transfert de masse introduit une erreur appelée
dispersion numérique ou erreur locale de troncature ().

X X

( i, j+ 1 )

u. t
ux. t uy t
Y Y ( i, j)
( i, j) ( i+ 1 , j) ux. t ( i+ 1 , j)

A- déplacement parallèle B- déplacement angulaire

Fig. – Représentation schématique des causes de dispersion numérique par la

méthode des différences finies.
Au premier ordre près, le développement de l’équation de dispersion se traduit par une
expression obtenue en appliquant successivement un développement de Taylor (schéma
décentré en amont) :
 C  1  2 
2  C   ....
C x  Δ x   C x   Δ x    Δ x  
 x  2  x 2 
 
Il s’ensuit que :
C C x  Δx   C x  1   2C 
 - Δx    ....
x Δx 2  x 2 
 
Ar cnséquent e tere

C C t  Δt   C t  1   2C 
 - Δt    ....
t Δt 2  t 2 

C C t  Δt   C t  1 2   2 C 
 - u Δt  ....
t Δt 2  x 2 
 
C C x  Δx   C x - Δx     3C 
  0 Δx  2 . 
t 2ΔΔ   x 3 
  
 1 Δt  1   2 C
ε L  uΔΔ σ    u θ  
 2 Δx  2  x 2
   2C 
ε L  D num  
 x 2 
 
 1 Δt  1 
avec D num  uΔΔ σ    u θ  
 2 Δx  2 
0 σ 1
0 θ 1
2C
 D L  D num   u
C

C
x 2 x 2 t

D*  D L  1    
 1 2 
2uΔu  2u Δt 
 
D*  D L   1 2 
 2u Δ t 
D*  D L  1  
 1  
2uΔ x  2u 2Δ t 

C C S o lu t io n c a lc u lé e a v e c
S o lu t io n e x a c t e d is p e r s io n n u m é r iq u e

X S o lu t io n e x a c t e X
S o lu t io n c a lc u lé e
a v e c o s c illa t io n
u x u y y
Pe x  x et Pe y  ut
D xx D yy
 p u it s d e p ié z o m è t r e
p o m p a g e (Q ) d 'in je c t io n

e1 Q1

e2 Q2

Zone de Nuage de
p o llu t io n p o llu t io n

C A  2C A C u  2C
α L ε α T
t r r 2 r r r 2 θ 2
1 2π
avec C r, t    C r,θ, t  dθ
2π 0
Q
A
2ππω
ε  sign  Q 
ε  1
uA
r
u  2C
αT
r 2 θ 2
C A  2C A C
α L ε
t r r 2 r r

Ci t Δt  Ci t   α L A
DD 
A
CC
Δt rΔΔ2 rΔΔ

avec DD  δ  C i 1  2C i 1  t Δt  1  δ  C i  1  2C i  C i 1 t 
   
CC  1  λ  γ Ci 1t Δ t  1  γ  Ci 1t   2λ  1 γC i t Δ t  1  γ  C i t  γC i 1t Δ t  1  γ  Ci 1t 
0 δ 1
 0 γ 1
0 λ 1
i, j+ 1 i+ 1 , j+ 1

i-1 , j+ 1

i-1 ,j i,j i+ 1 , j

i-1 , j-1 i, j-1 i+ 1 , j-1

: nœ ud
: p a r t ic u le
: m o u v e m e n t c o n v e c t if
: m o u v e m e n t d is p e r s if
X p t  dt  X p t  u xΔt

X p t  dt  X p t  u yΔt
 Cij
Δ  Cij
C p t  dt  C p t  ΔC ij t
 Quelques éléments de calcul d’une EDR santé
DEFINITION DES CONCENTRATIONS D’EXPOSITION

- Concentration accumulées dans les légumes :

C  C sol.BCF sol  plante

- Concentration dans les tissus musculaires de la vache :

C vache  C vache- plante  C vache-sol

C vache - plante  C plante .BTF  plante  vache
C vache- sol  C sol.BTF  sol  vache

- Concentration dans le lait de vache :

C lait - vache  C vache- plante  C vache-sol

C vache - plante  C plante .BTF Plante  lait de vache 
C vache- sol  C sol .BTF Sol  lait de vache

ESTIMATION DES DOSES D’EXPOSITION (DJE)

DJE pour l’ingestion de sol :

C s .IR.0,000001.EF.O.Z
DJE 
BW.365

DJE pour l’inhalation des poussières :

C s .B o .0,000001.EF
DJE 
BW.365
avec B o  C sol .0,07.10  6

DJE pour contact dermique avec le sol :

C s .SA.Fs.SL.J.0,000001.EF.Z
DJE 
BW.365
 4.BW   7
avec SA 
BW  90

DJE pour l’ingestion des légumes cultivés sur les sols pollués :

 Cs . (FR .VR .HG R )   FNR .VNR .HG NR   

DJE 
BW
DJE pour l’ingestion de viande bovine :
 C vache.IR.EF
DJE 
BW.365

DJE pour l’ingestion de lait de vache :

C vache.IR.EF
DJE 
BW.365

ETUDE DES EFFETS TOXIQUES

ESTIMATION DES RISQUES

- Estimation du risque non cancérigène (développement de tumeur, leucémie, …)

DJE
IR 
DJA

Si IR 1  risque toxique considéré comme acceptable

- Estimation du risque cancérigène

ERI  DJE.ERU

avec ERI : excès de risque individuel

ERU : excès de risque unitaire
Si ERI  10  5 (1 sur 1 million)  risque certainement acceptable
Si ERI  10  4 (1 sur 1 million)  risque certainement inacceptable
 MODELE D’ECOULEMENT A REGIME TRANSITOIRE
BIDIMENSIONNEL “ASM” AVEC TRANSPORT

A. Etape d’introduction des données

1ère fenêtre : choix possibles : 8 (MAIN MENU)

2ème fenêtre : elle concerne l’étape n°2

3ème fenêtre : elle

4ème fenêtre :

5ème fenêtre :

6ème fenêtre :

7ème fenêtre :

8ème fenêtre :

B. Etape de calcul
1ère méthode :
2ème méthode: