Annexes

Annexe I
Éléments de la Mécanique des Fluides

L’analyse de l’écoulement des eaux souterraines nécessite une compréhension des éléments de la mécanique des fluides. Albertson et Simons (1964) offrent une révision abréviée utile, tandis que Streeter (1962) et Yennard (1961) sont des textes de référence. Notre but ici est de présenter les propriétés de base des fluides : la densité de masse, poids, densité, compressibilité et la viscosité, ainsi que d’examiner les concepts de pression de fluide et une charge hydraulique.

Une discussion des principes de la mécanique des fluides doit commencer par une revue de la mécanique des matériaux en général. Le Tableau A1.1 fournit une liste des propriétés mécaniques de base de la matière, avec leurs dimensions et les unités dans le système métrique SI. Le système SI a comme dimensions de base la masse, la longueur, et de l’heure, tandis que les unités SI de base sont le kilogramme (kg), le mètre (m), et la seconde (s). Toutes les autres propriétés sont mesurées en unités qui sont dérivés de cet ensemble de base. Certaines de ces propriétés sont si communes, et leurs dimensions de base si complexe, que des noms spéciaux ont été accordés pour leurs unités dérivées. Comme indiqué sur le Tableau A1.1, la force et le poids sont mesurés en newtons (N), la pression et le stress en N/m2 ou pascals (Pa), et le travail et l’énergie en joules (J).

Tableau A1.1 Définitions, dimension, et unités SI pour des propriétés mécaniques de base

Dimension de l’unité
Propriété Symbole Définition Unité SI Symbole SI Dérivée Fondamentale
Masse M kilogramme kg kg
Longueur l   mètre m m
Temps t   seconde s s
Aire A A = l2 m2
Volume V V = l3 m3
Vélocité v v = l/t m/s
Accélération a a = l/t2 m/s2
Force F F = Ma newton N N kg•m/s2
Poids w W = Mg newton N N kg•m/s2
Pression p P = F/A pascal Pa N/m2 kg/m•s2
Travail W W = Fl joule J N•m kg•m2/s2
Énergie   Travail effectué joule J N•m kg•m2/s2
Densité de masse \rho \rho = M/V kg/m3
Densité de poids \gamma \gamma = w/V N/m3 kg/m2•s2
Stress \sigma, \tau Réponse interne au p externe pascal Pa N/m2 kg/m•s2
Tension \varepsilon \varepsilon = \DeltaV/V Sans dimension
Module de Young E Loi de Hooke N/m2 kg/m•s2

Le Tableau A1.2 fournit une analyse SI des certaines propriétés de fluides ainsi que des termes concernant les eaux souterraines dans ce texte. Une description approfondie se trouve dans le Chapitre 2.

Une grande partie de la technologie associée au développement des ressources en eaux souterraines en Amérique du Nord est toujours basé sur le système d’unités FPS (pied-livre-seconde).

Tableau A1.2 Définitions, dimensions et unités SI pour les propriétés de fluides ainsi que les termes associés avec l’eau souterraine

Dimensions de l’unité
Propriété Symbole Définition Unité SI Symbole SI Dérivée Fondamentale
Volume V V = l3 litre
(= m3 × 10-3)		  \ell  \ell m3
Décharge Q Q = l3/t \ell/s m3/s
Pression de fluide p p = F/A pascal Pa N/m2 kg/m•s2
Charge h   m
Densité de masse \rho \rho = M/V kg/m­3
Viscosité
dynamique		 \mu Loi de Newton centipoise
(= N•s/m2 × 10-6)		 cP cP, N•s/m2 kg/m•s
Viscosité
cinématique		 v v = \mu/\rho centistoke
(=m2/s × 10-6)		 cSt cSt m2/s
Compressibilité \alpha,\beta \alpha = 1/E m2/N
Conductivité
hydraulique		 K Loi de Darcy cm/s
Perméabilité k k = K\mu/pg cm2 m2
Porosité n   Sans dimension
Stockage spécifique SS = pg(\alpha + n\beta) 1/m
Storativité S S = SSb* Sans dimension
Transmissivité T T = Kb* m2/s
*b, épaisseur de l’aquifère captif (voir Section 2.10).

Table A1.3 Facteurs de conversion des unités FPS (pied-livre-seconde) aux unités SI

Multiplier par pour obtenir
Longueur pi 3,048 × 10–1 m
pi 3,048 × 10 cm
pi 3,048 × 10–4 km
mile 1,069 × 103 m
mile 1,069 km
Aire pi2 9,290 × 10–2 m2
mi2 2,590 km2
acre 4,407 × 103 m2
acre 4,407 × 10–3 km2
Volume pi3 2,832 × 10–2 m3
Gal US 3,785 × 10–3 m3
Gal UK 4,546 × 10–3 m3
pi3 2,832 × 10
Gal US 3,785
Gal UK 4,546
Velocité pi/s 3,048 × 10–1 m/s
pi/s 3,048 × 10 cm/s
mi/h 4,470 × 10–1 m/s
mi/h 1,609 km/h
Accélération pi/s2 3,048 × 10–1 m/s2
Masse lbm* 4,536 × 10–1 kg
slug* 1,459 × 10 kg
ton 1,016 × 103 kg
Force et poids lbf* 4,448 N
poundal 1,383 × 10–1 N
Pressure et stress psi 6,895 × 103 Pa ou N/m2
lbf/pi2 4,788 × 10–1 Pa
poundal/pi2 1,488 Pa
atm 1,013 × 105 Pa
en Hg 3,386 × 103 Pa
mb 1,000 × 102 Pa
Travail et énergie pi-lbf 1,356 J
pi-poundal 4,214 × 10–2 J
Btu 1,055 × 10–3 J
calorie 4,187 J
Densité de masse lb/pi3 1,602 × 10 kg/m3
slug/pi3 5,154 × 102 kg/m3
Densité de poids lbf/pi3 1,571 × 102 N/m3
Décharge pi3/s 2,832 × 10–2 m3/s
pi3/s 2,832 × 10 \ell/s
Gal US/min 6,309 × 10–5 m3/s
Gal UK/min 7,576 × 10–5 m3/s
Gal US/min 6,309 × 10–2 \ell/s
Gal UK/min 7,576 × 10–2 \ell/s
Conductivité hydraulique
(voir aussi le Tableau 2.3)		 pi/s 3,048 × 10–1 m/s
Gal US/jour/pi2 4,720 × 10–7 m/s
Transmissivité pi2/s 9,290 × 10–2 m2/s
Gal US/jour/pi 1,438 × 10–7 m2/s
*Un corps dont la masse est de 1 livre de mass (lbm) a un poids de 1 livre de force (lbf) 1 lbf est la force nécessaire pour accélérer un corps de 1 lbm jusqu’à une accélération de g = 32,2 pi/s2. Un « slug » est une unité de masse qui acquiert une accélération de 1 pi/s2 lorsqu’une force de 1 lbf agit dessus.

Le Tableau A1.3 fournit un ensemble de facteurs de conversion afin to transformer les unités FPS en unités SI. La densité de masse (ou densité) \rho d’un fluide est défini comme la masse par unité de volume (Tableau A1.1). La densité de poids (ou poids spécifique, ou poids unitaire) \gamma d’un fluide est défini comme le poids par unité de volume. Les deux paramètres sont reliés par

\gamma = \rho g (A1.1)

Pour l’eau, \rho = 1,0 g/cm3 = 1000 kg/m3; \gamma = 9,8 × 103 N/m3. Dans le système FPS, \gamma = 62,4 lbf/pi3.

Le gravité spécifique G de n’importe quel matériel est le ratio de sa densité (ou poids spécifique) à celui de l’eau. Pour l’eau, G = 1,0; pour la majorité des sols et des roches, G \approx 2,65.

La viscosité d’un fluide est la propriété qui permet aux fluides de résister le déplacement relation et déformation de cisaillement lors de l’écoulement. Plus qu’un fluide est visqueux, plus grand que le stress de cisaillement sera pour n’importe quel gradient de vélocité. Selon la loi de viscosité de Newton,

\tau = \mu \frac{dv}{dy} (A1.2)

Ou \tau est le stress de cisaillement, dv/dy est le gradient de vélocité, \mu la viscosité, ou viscosité dynamique. La vélocité cinétique v est donné par

v = \frac{\mu}{\rho} (A1.3)

Ou \rho est la densité du fluide.

La compressibilité d’un fluide reflète ses propriétés de contrainte-déformation. Le stress est la réponse interne d’un matériau soumise à une pression externe. Pour les fluides, le stress est transmis à travers la pression du fluide. La tension est une mesure de la déformation linéaire ou volumétrique d’un matériau stressé. Pour les fluides, il prend la forme d’un volume réduit (et densité accrue) sous des pressions de fluide croissantes. La compressibilité de l’eau \beta est discuté en détail à la Section 2.9. Il est défini par l’équation (2.44).

La densité, la viscosité et la compressibilité de l’eau sont des fonctions de la température et pression (Dorsey, 1940, Weast, 1972). En général, leur variation n’est pas grande, et pour la gamme de pressions et de températures qui se produisent dans la plupart des applications concernant les eaux souterraines, il est courant de les considérer comme des constantes. A 15,5 °C, \rho = 1,0 g/cm3, \mu = 1,124 cP et \beta = 4,4 × 10–2 m2/N.

La pression du fluide p à tout moment dans un corps d’eau permanent est la force par unité de surface qui agit à ce moment-là. En conditions hydrostatiques, la pression du fluide à un point reflète le poids de la colonne d’eau recouvrant une unité d’aire de la section transversale autour de ce point. Il est possible d’exprimer la pression par rapport au pression zéro absolue, mais plus communément il est exprimé par rapport à la pression atmosphérique. Dans ce dernier cas, il est appelé pression jauge, car c’est la lecture de pression qui est obtenue sur des jauges mises à zéro dans l’atmosphère.

La charge de pression \psi à un point dans le fluide est la hauteur que cette colonne atteindra si nous plaçons un manomètre à ce point. Dans un corps d’eau permanent, \psi est égale à la profondeur du point de mesure en-dessous de la surface. Si p est exprimé comme pression de jauge, \psi est définit par la relation

p = \rho g\psi = \gamma\psi (A1.4)

En effet, la charge de pression \psi est une mesure de la pression de fluide p.

Les pressions de fluides se développent dans l’eau souterraines lors de l’écoulement dans les formations géologiques et sols poreux. Dans la Section 2.2, les éléments de la mécanique des fluides présentées dans cet Annexe sont appliqués dans le développement de la théorie de l’écoulement des eaux souterraines.

Annexe II
Équation du débit pour un écoulement transitoire due à la déformation des milieux saturés

Un développement rigoureux de l’équation de l’écoulement pour l’écoulement transitoire en milieu saturé doit reconnaître que les changements transitoires de la pression du fluide conduisent à déformations dans le squelette granulaire d’un milieu poreux, et ces déformations impliquent que le milieu, ainsi que l’eau, sont en mouvement. Cette prise de conscience crée le besoin de deux raffinements à la dérivation classique présentée par Jacob (1940) et suivi à la Section 2.11 de ce texte. Premièrement, comme l’a reconnu Biot (1955), il est nécessaire de projeter la loi de Darcy en termes de vitesse relative du fluide aux grains. En second lieu, comme l’a reconnu Cooper (1966), il est nécessaire d’envisager la conservation de masse pour le milieu ainsi que pour le fluide dans le volume de contrôle élémentaire. On peut développer les relations de continuité par une des trois façons suivantes : (1) en considérant un volume élémentaire déformant en coordonnées déformantes, (2) en considérant un volume élémentaire déformant en coordonnées fixes, ou (3) en considérant un volume élémentaire fixe en coordonnées fixes. D’après Gambolati et Freeze (1973), nous utiliserons un volume élémentaire fixe en coordonnées fixes. L’approche nécessite l’utilisation de la notation vectorielle et la dérivée matérielle (dérivée totale, dérivé important). Si ces concepts ne sont pas familiers, Aris (1962) et Wills (1958) fournissent des textes introductoires. Le développement sera présenté ici pour un milieu homogène et isotrope à conductivité hydraulique K, porosité n, et compressibilité verticale \alpha. La même approche est facilement adaptée aux milieux hétérogènes et anisotropes.

En notation vectorielle, la forme trois-dimensionnel de la loi de Darcy est

\vec{v} = -K\nabla h (A2.1)

Ou \vec{v} = (v_x, v_y, v_z) est la vélocité relative du fluide par rapport aux grains, et \nabla h = (\partial h/\partial x, \partial h/\partial y, \partial h/\partial z) est le gradient hydraulique. Nous pouvions développer ce vecteur \vec{v} comme

\vec{v} = n(\vec{v}_w - \vec{v}_s) (A2.2)

Ou \vec{v}_w est la vélocité du fluide et \vec{v}_s la vélocité du milieu déformante.

L’équation d’état pour l’eau (Éq. (2.47)) est

\rho = \rho_0 e^\beta^p (A2.3)

Et celui des grains de sol, qui sont incompressibles, est

\rho_s = constante (A2.4)

L’équation pour la continuité de l’eau est

-\nabla \cdot [np\vec{v}_w] = \frac{\partial}{\partial t}[np] (A2.5)

Et celui pour le sol est

-\nabla \cdot [(1 - n)\rho_s\vec{v}_s] = \frac{\partial}{\partial t}[(1 - n)\rho_s] (A2.6)

Dans ces équations \nabla \cdot est l’opérateur de divergence :

\nabla \cdot = \frac{\partial}{\partial x} + \frac{\partial}{\partial y} + \frac{\partial}{\partial z}

Développant l’Éq. (A2.5), nous arrivons à

-\rho\nabla \cdot (n\vec{v}_w) - n\vec{v}_w \cdot \nabla\rho = n\frac{\partial\rho}{\partial t} + \rho\frac{\partial n}{\partial t} (A2.7)

Si nous cancellons \rho_s de l’Éq. (A2.6) et nous réorganisons cette équation, nous obtenons une expression pour \partial n/\partial t. Ceci peut être substitué en Éq. (A2.7) ensemble avec une expression pour n\vec{v}_w obtenue de l’Éq. (A2.2). En divisant par \rho et en réorganisant, l’Éq. (A2.7) devient

-\nabla \cdot \vec{v} - (\frac{v}{\rho}) \cdot \nabla\rho = (\frac{n}{\rho})(\frac{\partial\rho}{\partial t} + \vec{v}_s \cdot \nabla\rho) + \nabla \cdot \vec{v}_s (A2.8)

Nous nous utilisons le dérivé matériel D/Dt = \partial/\partial t + \vec{v}_w \cdot \nabla, l’Éq. (A2.8) peut être exprimée par

-\nabla \cdot \vec{v} - (\frac{v}{\rho}) \cdot \partial\rho = \frac{n}{\rho}\frac{Dp}{Dt} + \nabla \cdot \vec{v}_s (A2.9)

Le premier terme du côté droit de l’Éq. (A2.9) peut être relié à la compressibilité de l’eau \beta par la relation

\frac{D\rho}{Dt} = \rho\beta \frac{Dp}{Dt} (A2.11)

Le dérivé matériel du côté droit de l’Éq. (A2.10) peut être remplacée par un dérivé partiel seulement si l’inégalité suivant est résoute :

\frac{D\rho}{Dt} = \rho\beta \frac{Dp}{Dt} (A2.11)

Nous assumons d’avantage que, du côté gauche de l’Éq. (A2.9),

\vec{v}_s \cdot \nabla p \ll \frac{\partial p}{\partial t} (A2.12)

Ensuite, la substitution des Éqs. (A2.1) et (A2.10) dans l’Éq.(A2.9) nous donne

\nabla \cdot (K\nabla h) = n\beta\frac{\partial p}{\partial t} + \nabla \cdot \vec{v}_s (A2.13)

Dans un champ de stress trois-dimensionnel, le vecteur de vélocité des grains \vec{v}_s = (v_s_x, v_s_y, v_s_z) est relié au vecteur de déformation (ou déplacement de sol) \vec{u} = (u_x, u_y, u_z) par

\vec{v}_s = \frac{D\vec{u}}{Dt} (A2.14)

Dans un champ de stress unidimensionnel,

v_s_x = v_s_y = u_x = u_y = 0 (A2.15)

Si les conditions de l’Éq. (A2.15) sont satisfaits, le dernier terme de l’Éq. (A2.13) peut être développé (Cooper, 1966; Gambolati et Freeze, 1973; Gambolati, 1973a) comme

\nabla \cdot \vec{v}_s = \frac{\partial v_s_z}{\partial z} = \frac{\partial}{\partial z}(\frac{Du_z}{Dt}) = \frac{D}{Dt}(\frac{\partial u_z}{\partial z}) = \alpha\frac{Dp}{Dt} (A2.16)

Ou \alpha est la compressibilité verticale du milieu poreux. La variation du dérivé autour de l’égalité centrale est valide pour un vecteur de position mais pas en général. Le dérivé matériel dans l’expression droit de l’Éq. (A2.16) peut être remplacée par la dérivé partielle si Éq. (A2.11) est satisfait. Dans ce cas, Éq. (A2.13) devient

\nabla \cdot (K\nabla h) = n\beta\frac{\partial p}{\partial t} + \alpha\frac{\partial p}{\partial t} (A2.17)

Étant donné \rho = \rho g(h - z) et que K est une constante, l’Éq. (A2.17) est simplifié à

\nabla^2h = \frac{\rho g(\alpha +n\beta}{K}\frac{\partial h}{\partial t} (A2.18)

Ou bien, en se rappelant que S_s = \rho g(\alpha + n\beta) et en développant la notation vectorielle,

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} = \frac{S_s}{K}\frac{\partial h}{\partial t} (A2.19)

L’équation (A2.19) est identique à l’Éq. (2.75) développé par Jacob (1940). Ce développement plus rigoureux souligne le fait que la validité de l’équation classique de flux repose sur la satisfaction des inégalités des Eqs. (A2.11) et (A2.12) ainsi que la condition de stress de l’équation (A2.15). Il est peu probable que ces conditions soient toujours satisfaites. Gambolati (1973b) a démontré que lorsque le taux de consolidation \vec{v}_s dépasse le taux de fluide percolant \vec{v}_w (qui est possible dans les épaisses couches d’argile de faible perméabilité et haute compressibilité), il se peut que les inégalités ne soient pas satisfaites. En ce qui concerne la condition de stress, le terme \nabla \cdot \vec{v}_s à la fin de l’équation (A2.13) n’est qu’une petite partie qui relie le champ d’écoulement en trois dimensions au champ de contraintes en trois dimensions. Biot (1941, 1955) a d’abord exposé les interrelations et Verruijt (1969) fournit une dérivation claire. Schiffman et al. (1969) fournissent une comparaison de l’approche classique et celui de Biot, et Gambolati (1974) a analysé la gamme de validité de l’équation classique de l’écoulement.

Annexe III
Exemple d’une solution analytique pour un problème de valeur-limite

Considérons l’étude de cas simple présenté dans la Figure 2.25 (a). L’équation de l’écoulement en régime permanent dans le plan xy est

\frac{\partial^2 h}{\partial x^2} + \frac{\partial^2 h}{\partial y^2} = 0 (A3.1)

L’expression mathématique des conditions aux limites est

\frac{\partial h}{\partial y} = 0 sur y = 0 et y = y_L (A3.2)

h = h_0 on x = 0 (A3.3)

h = h_1 on x = x_L (A3.4)

Nous allons résoudre pour h(x, y) en séparant les variables.

Assumons que la solution est une solution de produit ayant la forme

h(x, y) = X(x) \cdot Y(y) (A3.5)

L’équation (A3.1) devient

Y\frac{\partial^2x}{\partial x^2} + X\frac{\partial^2y}{\partial y^2} = 0 (A3.6)

En divisant par XY, nous obtenons

\frac{1}{X}\frac{\partial^2X}{\partial x^2} = \frac{1}{Y}\frac{\partial^2Y}{\partial y^2} (A3.7)

Le côté gauche est indépendant de y. Par conséquent, le côté droit, malgré son apparence, doit aussi être indépendant de y, car il est égal au côté gauche. De même, les deux côtés sont indépendants de x. Puisque les deux côtés sont indépendants, chaque côté doit représenter une constante. Alors,

\frac{1}{X}\frac{\partial^2X}{\partial x^2} = G and \frac{1}{Y}\frac{\partial^2Y}{\partial y^2} = G (A3.8)

La constante G peut être positive, négative, ou nulle. Chacun des trois cas mènent à des solutions de produit, mais seulement le cas G = 0 mène à une solution qui est physiquement significatif pour ce problème. Les Équations (A3.8) deviennent

\frac{1}{X}\frac{\partial^2X}{\partial x^2} = 0 and \frac{1}{Y}\frac{\partial^2Y}{\partial y^2} = 0 (A3.9)

Ce sont des équations différentielles pour lesquelles les solutions sont bien connues

X = Ax + B and Y = Cy + D (A3.10)

La solution de produit devient

h(x, y) = (Ax + B)(Cy +D) (A3.11)

Nous pouvons évaluer les coefficients A, B, C et D et invoquant les conditions aux limites. La différenciation de l’Éq (A3.11) par rapport à y donne

\frac{\partial h}{\partial y} = (Ax + B)C (A3.12)

Et l’invocation de l’Éq. (A3.2) implique que C = 0. De l’Éq. (A3.11), nous retenons

h(x, y) = (Ax + B)D = Ex + F (A3.13)

Invoquant les conditions limites des Éqs. (A3.3) et (A3.4) nous donne F = h0 et E = -(h_0 - h_1)/x_L. Alors, la solution est

h(x, y) = h_0 - (h_0 - h_1)\frac{x}{x_L} (A3.14)

Cette équation est identique à l’Éq. (2.81), présenté sans la preuve dans la Section 2.11.

Il est évident que l’Éq. (A3.14) satisfait les conditions aux limites des Éqs. (A3.3) et (A3.4). La différenciation par rapport à y donne zéro en réalisation de l’Éq. (A3.2). La différenciation deux fois par rapport à x donne aussi zéro, alors la solution Éq. (A3.14) satisfait l’équation de l’écoulement (Éq. (A3.1)).

Annexe IV
L’Équation Debye-Hückel et Tableau Kielland pour les coefficients d’activité ionique

L’expression Debye-Hückel pour les activités ioniques individuels :

log \gamma = \frac{-Az^2\sqrt{I}}{1 + \r{a}B\sqrt{I}}

Valeurs du paramètre pour la taille d’ion å pour des ions communs dans les eaux naturels :

å × 108 Ion
2,5  \ce{NH^+_4}
3,0 \ce{K+}, \ce{Cl-}, \ce{NO^-_3}
3,5 \ce{OH-}, \ce{HS-}, \ce{MnO^-_4}, \ce{F-}
4,0 \ce{SO_4^{2-}}, \ce{PO_4^{3-}}, \ce{HPO_4^{2-}}
4,0–4,5 \ce{Na+}, \ce{HCO^-_3}, \ce{H2PO^-_4}, \ce{HSO^-_3}
4,5 \ce{CO3^{2-}}, \ce{SO3^{2-}}
5 \ce{Sr^{2+}}, \ce{Be^{2+}}, \ce{S^{2-}}
6 \ce{Ca^{2+}}, \ce{Fe^{2+}}, \ce{Mn^{2+}}
8 \ce{Mg^{2+}}
9 \ce{H+}, \ce{Al^{3+}}, \ce{Fe^{3+}}

Paramètres A et B à 1 Bar

Température (°C) A B (× 10–8)
0 0,4883 0,3241
5 0,4921 0,3249
10 0,4960 0,3258
15 0,5000 0,3262
20 0,5042 0,3273
25 0,5085 0,3281
30 0,5130 0,3290
35 0,5175 0,3297
40 0,5221 0,3305
50 0,5319 0,3321
60 0,5425 0,3338

Tableau Kielland de coefficients d’activité ioniques à 25 °C, organisés par la taille des ions (basé sur l’équation Debye-Hückel)

Charge Taille,* å Ions l =
0,0005		 0,001 0,0025 0,005 0,01 0,025 0,05 0,1
1 2,5 \ce{Rb+}, \ce{Cs+}, \ce{Ag+}, \ce{NH^+_4}, \ce{Tl+} 0,975 0,964 0,945 0,924 0,898 0,85 0,80 0,75
3 \ce{K+}, \ce{Cl-}, \ce{Br-}, \ce{I-}, \ce{CN-}, \ce{NO^-_3}, \ce{NO^-_2}, \ce{OH-}, \ce{F-}, \ce{ClO^-_4} 0,975 0,964 0,945 0,925 0,899 0,85 0,805 0,755
4 \ce{Na+}, \ce{IO^-_3}, \ce{HCO^-_3}, \ce{HSO^-_3}, \ce{H2PO^-_4}, \ce{CLO^-_2}, \ce{C2H3O^-_2} 0,975 0,964 0,947 0,928 0,902 0,86 0,82 0,775
6 \ce{Li+}, \ce{C6H3O^-_2} 0,975 0,965 0,948 0,929 0,907 0,87 0,835 0,80
9 \ce{H+} 0,975 0,967 0,950 0,933 0,914 0,88 0,86 0,83
2 4.5 \ce{Pb^{2+}}, \ce{Hg2^{2+}}, \ce{SO4^{2-}}, \ce{CrO4^{2-}}, \ce{CO3^{2-}}, \ce{SO3^{2-}}, \ce{C2O4^{2-}}, \ce{S2O3^{2-}}, H citrate2– 0,903 0,867 0,805 0,742 0,665 0,55 0,455 0,37
5 \ce{Sr^{2+}}, \ce{Ba^{2+}}, \ce{Cd^{2+}}, \ce{Hg^{2+}}, \ce{S^{2-}}, \ce{WO4^{2-}} 0,903 0,868 0,805 0,744 0,67 0,555 0,465 0,38
6 \ce{Ca^{2+}}, \ce{Cu^{2+}}, \ce{Zn^{2+}}, \ce{Sn^{2+}}, \ce{Mn^{2+}}, \ce{Fe^{2+}}, \ce{Ni^{2+}}, \ce{Co^{2+}}, phthalate2– 0,905 0,870 0,809 0,749 0,675 0,57 0,485 0,405
8 \ce{Mg^{2+}}, \ce{Be^{2+}} 0,906 0,872 0,813 0,755 0,69 0,595 0,52 0,45
3 4 \ce{PO4^{3-}}, \ce{PO4^{3-}} 0,796 0,725 0,612 0,505 0,395 0,25 0,16 0,095
9 \ce{Al^{3+}}, \ce{Fe^{3+}}, \ce{Cr^{3+}}, \ce{Sc^{3+}}, \ce{In^{3+}}, et les terres rares 0,802 0,738 0,632 0,54 0,445 0,325 0,245 0,18

*Notez que ces tailles sont des valeurs arrondies pour la taille effective en solution aqueux, et ne sont pas la taille des ions simples non-hydratés. Pour une discussion détaillé, voir l’article original duquel ces valeurs sont tirés [J. Kielland, J. Amer. Chem. Soc., 59, 1675 (1973)].

Les références

BERNER, R. A. 1971. Principles of Chemical Sedimentology. McGraw-Hill, New York, 240 pp.

KLOTZ, I. M. 1950, Chemical Thermodynamics. Prentice-Hall, Eaglewood Cliffs, N.J., 369 pp.

MANOV, G. G., R. G. BATES, W. J HAMER, et S. F. ACREE. 1943. Values of the constants in the Debye-Hückel equation for activity coefficients. J. Amer. Chem. Soc., 65, pp. 1765–1767.

Annexe VI
Développement d’une équation différences-finies pour l’équation de l’écoulement en régime permanent dans un milieu homogène et isotrope

L’équation différentielle partielle décrivant l’écoulement en régime permanent dans une région d’écoulement deux-dimensionnel, homogène, et isotrope est l’équation de Laplace :

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial y^2} (A6.1)

Afin de trouver l’équation différences finies pour un nœud intérieur dans la grille nodulaire utilisée pour discrétiser la région d’écoulement, nous devons remplacer les dérivées partielles de deuxième ordre dans l’Éq. (A6.1) par des différences. Commençons par considérant le premier terme de l’équation. En se rappelant que la définition de la dérivé partielle par rapport à x d’une fonction de deux variables h(x, z) est

\frac{\partial h}{\partial x} = \substack {lim \\ \Delta x \rightarrow 0} \frac{h(x + \Delta x, z) - h(x, z)}{\Delta x} (A6.2)

Sur un ordinateur il est impossible de définir la limite comme \Delta x \rightarrow 0, mais il est possible d’approximer la limite en assignant une petite valeur arbitraire à \Delta x. En fait, nous pouvons effectuer ceci en désignant un réseau nodulaire avec un espacement de maille de \Delta x.

Pour n’importe quelle valeur de z, comme z0, nous pouvons développer h(x ,z0) avec une expansion de Taylor par le point (x0, z0) comme ce qui suit :

h(x, z_0) = h(x_0, z_0) + (x - x_0)\frac{\partial h}{\partial x}(x_0, z_0) + \frac{(x - x_0)^2}{2}\frac{\partial^2h}{\partial x^2}(x_0, z_0) + ... (A6.3)

Si x = x_0 + \Delta x (reconnue comme la différence vers l’avant) et nous abandonnes tous les termes d’ordre supérieure à l’unité, nous pouvons approximer \partial h/\partial x par

\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - h(x_0, z_0)}{\Delta x} (A6.4)

Les termes abandonnés de l’expansion de Taylor représentent l’erreur de troncation dans l’approximation différences-finies.

Nous pouvons obtenir une expression similaire à l’Éq. (A6.4) en substituant la différence vers l’arrière, x = x_0 - \Delta x, dans l’Éq. (A6.3). Ceci donne

\frac{\partial h}{\partial x}(x_0, z_0) = \frac{h(x_0, z_0) - h(x_0 - \Delta x, z_0)}{\Delta x} (A6.5)

Afin d’obtenir l’approximation pour \partial^2h/\partial x^2, nous écrivons l’équation des différences en termes de \partial h/\partial x en se servant d’une expression différence-vers l’avant :

\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{\frac{\partial h}{\partial x}(x_0 + \Delta x, z_0) - \frac{\partial h}{\partial x}(x_0, z_0)}{\Delta x} (A6.6)

Et nous substituons l’expression différence vers l’arrière de l’Éq. (A6.5) dans l’Éq. (A6.6) pour obtenir

\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0 + \Delta x, z_0) - 2h(x_0, z_0) + h(x_0 + \Delta x, z_0)}{(\Delta x)^2} (A6.7)

De même, nous pouvons développer une expression de différence pour \partial^2h/\partial x^2 comme

\frac{\partial^2h}{\partial x^2}(x_0, z_0) = \frac{h(x_0, z_0 + \Delta z) - 2h(x_0, z_0) + h(x_0, z_0 - \Delta z)}{(\Delta z)^2} (A6.8)

Pour une grille carrée, \Delta x = \Delta z; l’ajout des Éqs. (A6.7) et (A6.8) pour former l’équation de Laplace nous donne

\frac{1}{(\Delta z)^2}[h(x_0 + \Delta x, z_0) + h(x_0 - \Delta x, z_0) + h(x_0, z_0 + \Delta z) + h(x_0, z_0 - \Delta z) + 4h(x_0, z_0) = 0] (A6.9)

Si (x0, z0) est le point nodulaire (i, j), l’Éq. (A6.9) peut être réarrangé afin d’obtenir

h_{ij} = \frac{1}{4}(h_{i+1,j} + h_{i-1,j} + h_{i,j+1} + h_{i,j-1}) (A6.10)

Qui est identique à l’Éq. (5.24).

Annexe VII
La solution analytique de Tóth pour l’écoulement régional

Tóth (1962, 1963) a présenté deux solutions analytiques pour le problème valeur-limite représentant l’écoulement en régime permanent dans un champ verticale, deux-dimensionnel, saturé, homogène et isotrope délimité au sommet par une nappe phréatique et des trois autres côtés par des frontières imperméables (la cellule ombrée dans la Figure 6.1, reproduit dans un système de coordonnées xz dans la Figure A7.1).

Figure A7.1 Région d’écoulement pour la solution analytique de Tóth.

En premier lieu, il a considéré le cas où la configuration de la nappe phréatique est une ligne droite ayant une pente constante. Dans ce cas, la région d’écoulement dans la Figure A7.1 est la région ABCEA. Comme il n’est pas possible d’obtenir une solution analytique dans une région ayant la forme d’une trapézoïde, Tóth a approché la région réelle du flux par la région ombrée ABCDA. Il a projeté les valeurs de charge hydraulique qui existent le long de la nappe phréatique réelle AE sur la limite supérieure AD de la région de solution. L’approximation est satisfaisante pour petit \alpha.

L’équation du flux est l’équation de Laplace :

\frac{\partial^2h}{\partial x^2} + \frac{\partial^2h}{\partial z^2} = 0 (7.1)

Pour la région délimitée par x = s et z = z0, les conditions aux limites sont

\frac{\partial h}{\partial x}(0, z) = \frac{\partial h}{\partial x}(s, z) = 0 sur AB et CD

\frac{\partial h}{\partial x}(x, 0) = 0 sur BC (7.2)

h(x, z_0) = z_0 + cx sur AD

Où c = tan \alpha.

La solution analytique, obtenue par la séparation de variables, est de

h(x, z) = z_0 + \frac{cs}{2} - \frac{4cs}{\pi^2}\sum_{m=0}^{\infty} \frac{\cos[(2m + 1)\pi x/s] \cos h[(2m + 1)\pi z/s]}{(2m +1)^2 \cos h[(2m +1)\pi z_0/s]} (7.3)

Cette équation satisfait à la fois l’équation de l’écoulement (A7.1) et les conditions aux limites (A7.2). Lorsque tracé et contourné, il conduit au réseau d’écoulement montré sur la Figure A7.1. Le débit va de la frontière de recharge DF à travers le champ jusqu’à la limite de décharge AF.

Tóth a également considéré le cas où la configuration de la nappe phréatique est spécifiée comme une courbe sinusoïdale superposés sur la ligne AE (pour représenter la topographie irrégulière). La condition limite finale devient alors

h(x, z_0) = z_0 + cx + \alpha \sin bx sur AD (A7.4)

Ou c = \tan \alpha, a = a'/\cos \alpha, et b = b'/\cos \alpha, a’ étant l’amplitude de la courbe sinusoïdale et b’ étant la fréquence.

La solution analytique dans ce cas a la forme de

h(x, z) = z_0 + \frac{cs}{2} + \frac{\alpha}{sb}(1 - \cos bs) + 2 \displaystyle\sum_{m=1}^{\infty} \left[\frac{ab(1 - \cos bs \cos m\pi)}{b^2 - \frac{m^2\pi^2}{s^2}} + \frac{cs^2}{m^2\pi^2}(\cos m\pi - 1)\right] \times \frac{\cos (m\pi x/s) \cos h (m\pi z/s)}{s \cdot \cos h (m\pi z_0/s)} (A7.5)

Dont le réseau d’écoulement décrit par cette fonction ressemble au Figure 6.2 (b).

Annexe VIII
Solution numérique du problème valeur-limite représentant l’infiltration en une dimension au-dessus d’un système d’écoulement d’eau souterraine rechargeant

Considérez un champ d’écoulement unidimensionnelle, verticale et homogène ayant la limite supérieure à la surface du sol et la limite inférieure en-dessous de la nappe phréatique. Comme discuté dans la Section 6.4, l’équation de l’écoulement dans un tel système est de

\frac{\partial}{\partial z}\left[K(\psi)(\frac{\partial\psi}{\partial z}+ 1)\right] = C(\psi)\frac{\partial\psi}{\partial t} (A8.1)

et les conditions aux limites sont

\frac{\partial\psi}{\partial t} = \frac{R}{K(\psi)} - 1 (A8.2)

au-dessus, et

\frac{\partial\psi}{\partial z} = \frac{Q}{K_0} - 1 (A8.3)

en-dessous. La solution est en termes de la charge de pression, \psi(z, t). La position de la nappe phréatique est donnée à tout moment par la valeur de z à laquelle *. La contribution nécessaire comprend les conditions initiales, \psi = 0, le taux de précipitation R, le taux de recharge des eaux souterraines Q, ainsi que les courbes non-saturées caractéristiques K(\psi) et C(\psi). Pour \psi \geq \psi_a, K = K0 et C = 0, \psi_a étant la charge de pression d’entrée d’air.

Le schéma numérique des différences finies dont Rubin et Steinhardt (1963), Liakopoulos (1965b), et Freeze (1969b) se sont servis est le schéma implicite de Richtmyer (1957). Avec cette méthode, le plan (z, t) est représenté par une grille rectangulaire de points j = 1, 2, …, L le long de l’axe z, et n = 1, 2, … le long de l’axe t. La distance entre les nœuds verticaux est \Delta z, et entre les pas de temps, \Delta t. La solution à tout point donné dans la grille (j, n) est \psi_j^n. Sous cette notation, la forme de la différence finie pour l’Eq. (A8.1) pour un nœud intérieur (j = 2 à j = L – 1), après arrangement approprié, est de

C(\psi_{\text{III}})\left(\frac{\psi_j^n - \psi_j^{n-1}}{\Delta t}\right) = \frac{1}{\Delta z}\left\{\left[K(\psi_{\text{I}})(1 + \frac{1}{2\Delta z})[\psi_{j+1}^{n-1} + \psi_{j+1}^{n} - \psi_{j}^{n-1} - \psi_j^n]\right]\right\} (A8.4)

où les valeurs de \psi_{\text{I}}, \psi_{\text{II}}, et \psi_{\text{III}}, qui déterminent les valeurs de K et C pour appliquer à un nœud donné pour un pas de temps donné, sont déterminées par l’extrapolation des étapes précédentes, comme décrit par Rubin et Steinhardt (1963).

Des équations différences-finies similaires, qui incorporent les conditions aux limites (A8.2) et (A8.3), peuvent être écrits pour les nœuds aux limites supérieure et inférieure (j = 1, j = L).

Pour chaque pas de temps, l’approximation différences-finies constitue un système de L équations algébriques linéaires en L inconnues. La forme générale est

-A_j\psi_{j+1}^n + B_j\psi_j^n - C_j\psi_{j-1}^n = D_j (A8.5)

où les coefficients A, B, C, et D varient avec le nœud (j = 1, j = 2 à L – 1, j = L), avec le pas de temps (n = 1, n = 2, n > 2), et la saturation (\psi \geq \psi_a, \psi < \psi_a). Les variables pour laquelle A, B, C, et D dépendent sont les valeurs de \psi du pas de temps précédent, les valeurs de limites R et Q, et les relations fonctionnels K(\psi) et C(\psi).

Les valeurs de * sont calculées de la relation de récurrence

\psi_j^n = E_j\psi_{j+1}^n + F_j (A8.6)

où

E_j = \frac{A_j}{B_j - C_j E_{j-1}}

E_j = \frac{D_j + C_j F_{j-1}}{B_j - C_j E_{j-1}}

Les coefficients E et F sont calculées de j = 1 à j = L, et les \psi’s sont recalculé de j = L à j = 1 en se servant de l’Éq. (A8.6). Pour j = 1, E_1 = A_1/B_1 et F_1 = D_1/B_1.

Si K(\psi) et C(\psi) sont hystérétiques, ils sont généralement intégrés dans un programme informatique sous la forme d’un tableau de valeurs représentant le mouillage principal, le séchage principal, et les courbes de balayage principales. Le programme doit localiser la courbe correcte en sachant si le nœud en question est en voie de mouillage ou séchage, et en scannant l’histoire du nœud pour déterminer la valeur de \psi à laquelle un changement de mouillage au séchage, ou vice versa, s’est produit.

Annexe IX
Développement de l’équation différences-finies pour l’écoulement transitoire dans un aquifère captif horizontale hétérogène et anisotrope

L’équation différentielle partielle qui décrit l’écoulement transitoire dans un milieu saturé anisotropique (Section 2.11) est

\frac{\partial}{\partial x} \left(K_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y} \left(K_y\frac{\partial h}{\partial y}\right) + \frac{\partial}{\partial z} \left(K_z\frac{\partial h}{\partial z}\right) = S_s\frac{\partial h}{\partial t} (A9.1)

Pour un aquifère captif horizontale d’épaisseur b, la forme deux-dimensionnelle de l’Éq. (A9.1) réduit ce dernier à

\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) + \frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial y}\right) = S_s\frac{\partial h}{\partial t} (A9.2)

où Tx et Ty sont les composantes principales du tenseur de transmissivité défini par T = Kb, et S la storativité définie par S = Ssb. Pour trouver l’équation des différences finies pour un nœud intérieur dans la grille nodale utilisée pour discrétiser la région d’écoulement, nous devons remplacer les dérivées partielles dans l’Éq (A9.2) par des différences. En utilisant les définitions développées à l’annexe VI et la notation ijk de la Section 8.8 ainsi qu’au Figure 8.26 (c), nous pouvons écrire l’expression de différence-finie :

\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) \simeq \frac{1}{\Delta x}\left[\left(T_x\frac{\partial h}{\partial x}\right)_{i+1/2,j}^k - \left(T_x\frac{\partial h}{\partial x}\right)_{i-1/2,j}^k\right] (A9.3)

où l’indice (i + \frac{1}{2}, j) signifie que la quantité entre parenthèses est évaluée au point médian entre les nœuds (i, j) et (i + 1, j), et l’exposant k signifie que la quantité entre parenthèses est évaluée au pas de temps k. Nous pouvons avancer d’avantage les termes sur le côté droit de l’Eq. (A9.3) par

\left(T_x\frac{\partial h}{\partial x}\right)_{i+1/2,j}^k = (T_x)_{i+1/2,j}(h_{i+1,j}^k - h_{i,j}^k)\frac{1}{\Delta x} (A9.4a)

\left(T_x\frac{\partial h}{\partial x}\right)_{i-1/2,j}^k = (T_x)_{i-1/2,j}(h_{i,j}^k - h_{i-1,j}^k)\frac{1}{\Delta x} (A9.4b)

Si nous évaluons (T_x)_{i+1/2,j} et (T_x)_{i-1/2,j} par moyennes de la forme

(T_x)_{i+1/2,j} \simeq \frac{(T_x)_{i+1,j} + (T_x)_{i,j}}{2} (A9.5)

Ces expressions peuvent être substituées dans les Éqs. (A9.4), et ces équations Peuvent alors être substitués dans l’Éq. (A9.3) afin de donner

\frac{\partial}{\partial x} \left(T_x\frac{\partial h}{\partial x}\right) \simeq \frac{1}{2\Delta x^2} \left\{\left[(T_x)_{i+1/2,j} + (T_x)_{i,j}\right]h_{i+1,j}^k - \left[(T_x)_{i+1/2,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i,j}^k + \left[(T_x)_{i,j} + (T_x)_{i-1,j}\right]h_{i-1,j}^k \right\} (A9.6)

Similairement :

\frac{\partial}{\partial y} \left(T_y\frac{\partial h}{\partial y}\right) \simeq \frac{1}{2\Delta y^2} \left\{\left[(T_y)_{i+1/2,j} + (T_y)_{i,j}\right]h_{i+1,j}^k - \left[(T_y)_{i+1/2,j} + 2(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i,j}^k + \left[(T_y)_{i,j} + (T_y)_{i-1,j}\right]h_{i-1,j}^k \right\} (A9.7)

Finalement, nous pouvons approximer le côté droit de l’Éq. (A9.2) par

S\frac{\partial h}{\partial t} \simeq S_{i,j}\left(\frac{h_{i,j}^k - h_{i,j}^{k-1}}{\Delta t}\right) (A9.8)

Substituant les Éqs. (A9.6), (A9.7), et (A9.8) pour les trois termes de l’Éq (A9.2) et l’assemblage des termes mène à l’équation différences-finies générale pour un nœud intérieur dans un aquifère hétérogène et anisotrope :

Ah_{i,j}^k = Bh_{i,j-1}^k + Ch_{i+1,j}^k + Dh_{i,j+1}^k + Eh_{i-1,j}^k + F (A9.9)

où

A = \frac{1}{2\Delta x}^2[(T_x)_{i+1,j} + 2(T_x)_{i,j} + (T_x)_{i-1,j}] + \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + 2(T_y)_{i,j} + (T_y)_{i,j-1}] + \frac{S_{i,j}}{\Delta t}

B = \frac{1}{2\Delta y^2}[(T_y)_{i,j} + (T_y)_{i,j-1}]

C = \frac{1}{2\Delta x^2}[(T_x)_{i+1,j} + (T_x)_{i,j}]

D = \frac{1}{2\Delta y^2}[(T_y)_{i,j+1} + (T_y)_{i,j}]

E = \frac{1}{2\Delta x^2}[(T_x)_{i,j} + (T_x)_{i-1,j}]

F = \frac{S_{i,j}}{\Delta t} \cdot h_{i,j}^{k-1}

Si l’aquifère est homogène et isotrope, alors Tx = Ty = T pour tout (i, j) et Si, j = S pour tout (i, j). Sous ces conditions, et pour une grille nodale ayant \Delta x = \Delta y, les coefficients de l’Éq. (A9.9) deviennent

A = \frac{4T}{\Delta x^2} + \frac {S_{i,j}}{\Delta t}

B = C = D = E = \frac{T}{\Delta x^2}

F = \frac{S}{\Delta t} \cdot h_{i,j}^{k-1}

Si nous divisons par T/\Delta x^2, ces coefficients sont considérés comme les mêmes développé de manière moins rigoureuse à la Section 8.8 et présenté en relation avec l’Eq. (8.57).

Trescott et al. (1976) ont suggéré qu’il y a un certain avantage à utiliser la moyenne harmonique plutôt que la moyenne arithmétique dans l’Éq. (A9.5). Cette approche change les coefficients dans l’équation aux différences finies mais ne modifie pas les concepts qui sous-tendent le développement.

L’équation (A9.9) est écrite en fonction des valeurs de la charge hydraulique à cinq nœuds au pas de temps k et un nœud au pas de temps k – 1. Il est connu comme l’approximation différence vers l’arrière. Remson et al. (1971) notent les avantages d’utiliser une approximation de différence centrale, connue sous le nom du schéma Crank-Nicholson, qui utilise des valeurs de charge à cinq nœuds à l’étape k et cinq nœuds à l’instant k – 1. La procédure implicite à direction alternative (ADIP) utilisé par Pinder et Bredehoeft (1968) implique deux équations aux différences finies, un dans le plan xt et un dans le plan yt. Chacun utilise des valeurs de charge à trois nœuds au pas de temps k et trois nœuds au pas de temps k – 1.

Annexe X
Dérivation de l’équation advection-dispersion pour le transport de solutés dans un milieu poreux saturé

L’équation dérivée ici est une déclaration de la loi de conservation de la masse. La dérivation est basée sur celles de Ogata (1970) et Bear (1972). On supposera que le milieu poreux est homogène et isotrope, que le milieu est saturé, que l’écoulement est de régime permanent, et que la loi de Darcy s’applique. Sous l’hypothèse de Darcy, le flux est décrit par la vitesse linéaire moyenne, qui porte le soluté par advection. Si c’était le seul mécanisme de transport opérationnel, les solutés non réactifs transportés par le fluide se déplaceraient comme un bouchon. En réalité, il existe un processus de mélange supplémentaire, soit la dispersion hydrodynamique (Section 2.13), qui est causée par des variations de la vitesse microscopique dans chaque canal de pore, ainsi que d’un canal à l’autre. Si nous souhaitons décrire le processus de transport à l’échelle macroscopique en utilisant des paramètres macroscopiques tout en incorporant l’effet du mélange microscopique, il est nécessaire d’introduire un deuxième mécanisme de transport, en plus de l’advection, pour tenir compte de la dispersion hydrodynamique.

Pour établir l’énoncé mathématique de la conservation de masse, le flux de soluté entrant et sortant d’un petit volume élémentaire dans le milieu poreux sera considéré (Figure A10.1). En coordonnées cartésiennes, la décharge spécifique v a les composants (vx, vy, vz) et la vitesse linéaire moyenne \vec{v} = v/n a des composants (\vec{v}_x, \vec{v}_y, \vec{v}_z). Le taux de transport par advection est égal à \vec{v}. La concentration du soluté C est défini comme la masse de soluté par unité de volume de solution. La masse de soluté par unité de volume de milieu poreux est donc nC. Dans un milieu homogène, la porosité n est une constante, et \partial(nC)/\partial x = n\partial C/\partial x.

Figure A10.1 Mass balance in a cubic element in space.

La masse de soluté transporté dans la direction x par les deux mécanismes de transport peut être représenté par

transport par advection = \vec{v}_xnC dA (A10.1)

transport par dispersion = nD_x\frac{\partial C}{\partial x} dA (A10.2)

Où Dx est le coefficient de dispersion dans la direction des x et dA est la section transversale élémentaire de l’élément cubique. Le coefficient de dispersion Dx est relié à la dispersivité \alpha_x et le coefficient de diffusion D* par l’Éq. (9.4) :

D_x = \alpha_x\vec{v}_x + D (A10.3)

La forme du composé dispersif inclut dans l’Éq. (A10.2) est similaire à la première loi de Fick.

Si Fx représente la masse totale de soluté par unité de l’aire transversale transporté dans la direction x par unité de temps, alors

F_x = \vec{v}_xnC - nD_x\frac{\partial C}{\partial x} (A10.4)

Le signe négatif qui précède le terme dispersif indique que le contaminant se déplace vers le zone de plus basse concentration. De manière similaire, les expressions dans les deux autres directions sont exprimées par

F_y = \vec{v}_ynC - nD_y\frac{\partial C}{\partial y} (A10.5)

F_z = \vec{v}_znC - nD_z\frac{\partial C}{\partial z} (A10.6)

Le montant total de soluté entrant l’élément cubique (Figure A10.1) est

F_x dz \hspace{1mm} dy + F_y dz \hspace{1mm} dx + F_z dx \hspace{1mm} dy

Le montant total de soluté sortant de l’élément cubique est

\left(F_x + \frac{\partial F_x}{\partial x}dx\right)dz \hspace{1mm} dy \left(F_y + \frac{\partial F_y}{\partial y}dy\right)dz \hspace{1mm} dx + \left(F_z + \frac{\partial F_z}{\partial z}dx \hspace{1mm} dy\right)

où les termes partiels indiquent la variation spatiale de la masse de soluté dans la direction spécifiée. La différence du montant entrant et sortant est alors

\left(\frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}\right) dx \hspace{1mm} dy \hspace{1mm} dz

Puisqu’on assume que le soluté est non-réactif, la différence entre le flux entrant et sortant de l’élément représente le montant de soluté accumulé dans l’élément. Le taux de changement de masse dans l’élément est

-n\frac{\partial C}{\partial t} dx \hspace{1mm} dy \hspace{1mm} dz

L’expression complète de la conservation de masse est alors

\frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z} = n\frac{\partial C}{\partial t} (A10.7)

La substitution des expressions (A10.4), (A10.5), et (A10.7) en (A10.7) et l’annulation de n des deux côtés nous donne :

\bigg[\frac{\partial}{\partial x}\left(D_x\frac{\partial C}{\partial x}\right) + \frac{\partial}{\partial y}\left(D_y\frac{\partial C}{\partial y}\right) + \frac{\partial}{\partial z}\left(D_z\frac{\partial C}{\partial z}\right)\bigg] - \bigg[\frac{\partial}{\partial x}(\vec{v}_xC) + \frac{\partial}{\partial y}(\vec{v}_yC) + \frac{\partial}{\partial z}(\vec{v}_zC)\bigg] = \frac{\partial C}{\partial t} (A10.8)

Dans un milieu homogène pour laquelle la vélocité \vec{v} est constante et uniforme (ne varie pas dans le temps et l’espace), les coefficients de dispersion Dx, Dy, et Dz ne varient pas dans le temps et l’espace, et l’Éq. (A10.8) devient

\bigg[D_x\frac{\partial^2C}{\partial x^2} + D_y\frac{\partial^2C}{\partial y^2} + D_z\frac{\partial^2C}{\partial z^2} \bigg] - \bigg[ \vec{v}_x\frac{\partial C}{\partial x} + \vec{y}_x\frac{\partial C}{\partial y} + \vec{v}_z\frac{\partial C}{\partial z} \bigg] = \frac{\partial C}{\partial t} (A10.9)

En une dimension,

D_x\frac{\partial^2C}{\partial x^2} - \vec{v}_x\frac{\partial C}{\partial x} = \frac{\partial C}{\partial t} (A10.10)

Pour quelques applications, la direction unidimensionnelle est considérée comme coordonné curvilinéaire dans la direction de l’écoulement sur une ligne d’écoulement. L’équation de transport devient alors

D_l\frac{\partial^2C}{\partial l^2} - \vec{v}_l\frac{\partial C}{\partial l} = \frac{\partial C}{\partial t} (A10.11)

Pour un problème deux dimensionnel, il est possible de définir deux directions de coordonnés curvilinéaires, Sl et St, ou Sl suit la ligne d’écoulement, et St est orthogonal à la ligne d’écoulement. L’équation de transport devient

D_l\frac{\partial^2C}{\partial S_l^2} + D_t\frac{\partial^2C}{\partial S_t^2} - \vec{v}_l\frac{\partial C}{\partial S_l} = \frac{\partial C}{\partial t} (A10.12)

Où Dl et Dt sont les coefficients de dispersion dans les directions longitudinaux et transverses, respectivement.

Si vl varie le long de la ligne d’écoulement et Dl et Dt varient dans l’espace, l’Éq. (A10.12) devient

\frac{\partial}{\partial S_l}\left(D_l\frac{\partial C}{\partial S_l}\right) + \frac{\partial}{\partial S_t}\left(D_l\frac{\partial C}{\partial S_t}\right) - \frac{\partial}{\partial S_l}(\vec{v}_lC) = \frac{\partial C}{\partial t} (A10.13)

Les équations (A10.8) à (A10.13) représentent six formes de l’équation dispersion-advection pour le transport de solutés dans des milieux poreux saturés. L’équation (A10.11) est identique à l’Éq. (9.3) dans la Section 9.2. La solution à n’importe quelle de ces équations fournira la concentration de soluté C en fonction de l’espace et du temps. Il prendra la forme C(x, y, z, t) pour les Éqs. (A10.8) et (A10.9); C(l, t) pour l’Éq. (A10.11) et C(Sl, St, t) pour les Éqs. (A10.12) et (A10.13). Il y a beaucoup de solutions analytiques bien connues pour les formes plus simples de l’équation de transport. L’équation (9.5) dans la Section 9.2 est une solution analytique à l’équation (A10.11), et l’Éq. (9.6) est une solution analytique à l’équation (A10.9). Dans la plupart des situations de terrain, des analyses en deux ou trois dimensions sont nécessaires. De plus, les vitesses sont rarement uniformes et les dispersivités sont généralement variables dans l’espace. Pour ces conditions, des méthodes numériques adaptées aux ordinateurs doivent être utilisées pour obtenir des solutions.

Le coefficient de dispersion dans un milieu poreux tridimensionnel, homogène, et isotrope tel qu’exprimé dans l’équation (A10.12) est un tenseur symétrique de second ordre avec neuf composants. Dl et Dt sont les termes diagonaux de la forme bidimensionnelle. Les propriétés directionnelles du coefficient de dispersion sont associées à la nature du processus dispersif dans les directions longitudinale et transversale. Si le milieu est lui-même anisotrope, une description mathématique du processus dispersif nécessite beaucoup plus de complexité. Aucune solution analytique ou numérique n’est disponible pour les systèmes anisotropes. Seuls les coefficients de dispersions en milieux isotropes ont été représenté par des méthodes numériques.

Il est possible d’étendre l’équation de transport pour inclure les effets du retard sur le transport de soluté par l’adsorption, les réactions chimiques, les transformations biologiques, ou désintégration radioactive. Dans ce cas, le bilan de masse effectué sur le volume élémentaire doit inclure un terme source-puits. Pour le retard dû à l’adsorption, l’équation de transport dans un milieu homogène dans un système unidimensionnel le long de la direction du flux prend la forme

D_l\frac{\partial^2C}{\partial l^2} - \vec{v}_l\frac{\partial c}{\partial l} + \frac{\rho b}{n}\frac{\partial S}{\partial t} = \frac{\partial C}{\partial t} (A10.14)

Où \rho_b est la densité de masse globale du milieu poreux, n la porosité, et S la masse du soluté adsorbé sur un unité de masse de la partie solide du milieu poreux. Le premier terme de l’Éq. (A10.14) est la terme de dispersion, le deuxième est la terme d’advection, et la troisième est le terme de réaction. Cette forme du terme de réaction est considérée dans la Section 9.2 en connexion avec les Éqs. (9.9) à (9.13). Des solutions analytiques pour l’Éq. (A10.14) sont présenté par Codell et Schreiber (en presse). La Figure 9.12 présente quelques exemples de solutions effectués par Pickens et Lennox (1976) en se servant de méthodes numériques.