Interstices


  Approfondir

L’eau sous nos pieds

Document paru initialement dans la revue Accromath.

Comment peut-on s’assurer que les déchets résidentiels ou industriels n’altéreront pas la qualité de l’eau de la nappe phréatique ? La modélisation mathématique de l’écoulement des eaux souterraines s’avère cruciale dans l’étude de l’impact environnemental de sites d’enfouissement ou de projets de développement industriel.

Comme nous en sommes de plus en plus conscients, l’activité humaine peut modifier profondément l’environnement. Et c’est souvent là où les interactions sont les moins visibles que les impacts peuvent être les plus grands. C’est le cas des eaux souterraines où des matières polluantes venant de sites d'enfouissement ou d’activités industrielles peuvent s’infiltrer et être transportées. Ce phénomène peut entraîner une contamination de l’eau destinée à la consommation ainsi que des cours d’eau, lacs et terres humides en contact avec cette eau souterraine.

Il revient normalement aux ingénieurs de s’assurer qu’un projet de développement ou d’utilisation d’un site tienne compte des exigences, pratiques et règlements, en matière d’environnement et de développement durable. Cela implique notamment la considération des impacts environnementaux potentiels à court et à long termes et l’élaboration de programmes de prévention de la pollution et de réduction des déchets.

Pour l’ingénieur, ces obligations à l’égard de l’environnement doivent être mises au premier plan, au même titre que celles liées à la santé et à la sécurité. L’évaluation du risque de contamination de l’eau de la nappe phréatique repose sur la modélisation mathématique et la simulation de l’écoulement des eaux souterraines. Des logiciels spécialisés existent pour effectuer ces simulations, mais la qualité des prédictions dépend de l’adéquation de la modélisation qui a été menée pour décrire le milieu à l’étude et de l’applicabilité pour ce milieu des méthodes et hypothèses utilisées par le logiciel.

La loi de Darcy

 

Parmi les éléments de modélisation mis à contribution dans ces logiciels, on retrouve la célèbre « loi de Darcy » qui permet de lier le débit dans un aquifère (couche du sous-sol saturée en eau) à la différence de pression entre deux points de l’aquifère. Cette différence de pression, directement liée à la différence de hauteurs (H1 - H2) qu’atteindrait l’eau dans un puits ouvert en ces points, engendre un mouvement de l’eau dans l’aquifère dont on peut mesurer le débit Q.

Ce principe s’applique aussi bien à grande échelle qu’à petite échelle. Pour isoler les paramètres et variables qui interviennent dans la modélisation du débit, considérons ce qui se passe au niveau d’un très petit volume de contrôle (entre les points a et b) le long d’une ligne de courant. Les différentes grandeurs susceptibles d’affecter le débit sont :

  • la différence de pression Pb – Pa,
  • la longueur L du volume de contrôle,
  • la section A du volume de contrôle,
  • la perméabilité K du milieu.

Comme l’a établi dans ses expériences Henry Darcy, si l’on fait varier une seule de ces grandeurs à la fois, on constate que le débit Q est :
a. proportionnel à la différence de pression ∆P = Pb - Pa,
b. proportionnel à la perméabilité K du milieu,
c. proportionnel à la surface A de la section de passage de l’eau dans le volume,
d. inversement proportionnel à la longueur L du volume

Darcy a donc postulé la relation suivante entre toutes ces variables :

Q = -KA (Pb - Pa)/L

Le signe négatif permet de rendre compte de la direction de l’écoulement : un fluide s’écoule toujours d’un point de haute pression (Pa) vers un point de basse pression (Pb). C’est d’ailleurs ce qui se passe avec les systèmes météorologiques.

On peut établir un lien entre cette équation et l’équation ∆V = RI qui lie, en électricité, un courant I à une différence de tension V et une résistance R. En considérant le milieu, caractérisé par le rapport L/(KA) comme une résistance à l’écoulement de l’eau souterraine, on a bien, au signe près, ∆P = RQ.

Lorsqu’on réduit la taille du volume de contrôle (section et longueur) et qu’on passe à la limite, on peut réécrire l’équation de Darcy sous forme différentielle et caractériser ainsi la vitesse de filtration u le long de la direction principale :

u = -K dP/dx

avec u = limA→0Q/A et dP/dx = limL→0(Pb - Pa) /L

Un milieu saturé

Du fait de leur âge respectable, bien des aquifères constituent un milieu saturé, qui ne peut absorber davantage d’eau : cette conservation de la masse d’eau fait en sorte qu’en tout point du milieu, s’il n’y a pas d’autres sources ou puits impliqués, la quantité d’eau qui entre doit être égale à celle qui sort.

Cette propriété équivaut à une variation nulle de la vitesse

du/dx = 0 et donc d (K dP/dx) / dx = 0.

Si la composition d’un aquifère ne varie pas dans l’espace, le milieu est dit homogène et K est constant à travers le domaine. On peut alors procéder à une simplification de l’équation ci-dessus (d2P/dx2 =0) et obtenir par conséquent une pression et une hauteur toutes deux linéaires en fonction de x. Dans le cas d’un écoulement unidirectionnel dans un aquifère hétérogène, K dépend de x et les fonctions de pression et de hauteur ne seront linéaires que par morceaux.

Passage au 3D

Le modèle précédent supposait qu’on pouvait ramener l’écoulement à une seule dimension, dans la direction de la pente de l’aquifère. Comme ça n’est pas toujours le cas, il convient de généraliser le modèle à l’espace tridimensionnel.

En considérant maintenant la vitesse de filtration (ou le flux) comme un champ de vecteurs u dans l’espace, la loi de Darcy en 3D fait intervenir le gradient de la fonction pression, noté ∇P, et le vecteur associé à la gravité ρ g (avec g = (0, 0, -g), ρ étant la masse volumique du fluide). La loi de Darcy s'écrit alors :

u = –K (∇P – ρ g)

Le flux va dans la direction où la chute de pression est la plus grande, c’est-à-dire celle opposée au gradient de pression.

Notons de plus que si le milieu n’est pas homogène, K devient ici une matrice. La conservation de la masse peut ici s’écrire en faisant appel au concept de divergence relativement au champ de vecteurs associé au flux.

div u = ∇.u = 0.

En résolvant ce système d’équations, on trouve le champ de vecteurs-vitesse u qui détermine les lignes de courant et donc la trajectoire de l’eau souterraine.

Comme, en général, il n’existe pas de solution analytique, on fait appel à des méthodes numériques, programmées dans les logiciels spécialisés. La qualité des résultats dépend de la finesse avec laquelle on a caractérisé le milieu et de l’adéquation des méthodes utilisées à reproduire dans le transfert de l’information numérique le flux physique en jeu.

En complément, voici un exemple d'images de simulation (calculées par Jocelyne Erhel).

  
 
Image de gauche : la perméabilité est homogène, égale à 1 dans tout le carré. On voit le maillage en petits carrés utilisé pour faire les simulations. On impose une pression égale à 1 à gauche, égale à 0 à droite ; on impose aussi un flux nul en haut et en bas (l'eau ne peut pas sortir). On s'attend donc à ce que l'eau coule de la forte pression vers la plus faible pression, c'est-à-dire de gauche à droite, avec une vitesse horizontale et constante. On sait aussi que la pression va diminuer régulièrement de gauche à droite, de 1 à 0, en étant constante sur les lignes verticales.
 
Image de droite : résultat de la simulation avec une perméabilité homogène. On obtient bien le comportement attendu. Cette simulation permet de vérifier que dans un cas simple, dont on connaît la solution, la simulation informatique est correcte. C'est une première étape de validation.
 
  
 
Image de gauche : On utilise le même maillage et les mêmes conditions aux bords, mais la perméabilité est hétérogène. Dans le rectangle jaune, la perméabilité vaut 100, donc le milieu est ici plus perméable et on s'attend à ce que l'eau coule plus vite dans ce rectangle. Dans le carré rouge, la perméabilité vaut 0,1 ; on s'attend à ce que l'eau y coule moins vite, et contourne ce carré, comme un obstacle.
 
Image de droite : résultat de la simulation avec cette perméabilité hétérogène. On voit que la vitesse a maintenant une composante verticale et qu'elle est plus grande dans le rectangle jaune, plus petite dans le carré rouge. On voit aussi que la charge n'est plus constante verticalement. On obtient bien un comportement qualitatif conforme à la prévision, mais il n'était pas possible de prévoir les valeurs quantitatives sans effectuer les calculs. Pour valider, on vérifie des invariants comme la conservation globale de la masse (il n'y a pas de fuite d'eau). Ici, le flux entrant à gauche est égal au flux sortant à droite. À l'échelle d'un petit carré du maillage, on a aussi conservation, dite locale, de la masse.
 

Une première version de ce document est parue dans la revue Accromath réalisée par l'Institut des sciences mathématiques et le Centre de recherches mathématiques du Québec, Volume 6.1, Hiver-printemps 2011.

Attention, avec Internet Explorer, le rendu de certaines formules mathématiques pourrait être altéré.