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.

">
Les Newsletters Interstices
    Niveau intermédiaire
    Niveau 2 : Intermédiaire

    L’eau sous nos pieds

    Environnement & Planète
    Modélisation & Simulation

    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.

    Entre 2008 et 2010, le projet d’un site d’enfouissement dans le comté de Simcoe (Ontario) s’est heurté à une vigoureuse opposition de citoyens inquiets de la menace que faisait planer ce site sur la qualité de l’eau potable de la région. Une étude réalisée par une firme d’ingénieurs à partir de simulations réalisées avec le logiciel MODFLOW (logiciel de modélisation des aquifères) soutenait que le site ne représentait pas un danger pour la nappe phréatique. S’appuyant sur la loi sur l’accès à l’information, le regroupement de citoyens demanda à connaître le modèle et les données utilisés par la firme pour les simulations. N’étant pas en mesure de se conformer à cette demande, le comté dut reculer dans son projet.

    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.


    Né à Dijon, Henry Darcy perdit son père à l’âge de 14 ans. Grâce à la détermination de sa mère, il put fréquenter les meilleures écoles. Remarqué pour sa facilité à concevoir et son assiduité au travail, il obtint à 23 ans son diplôme d’ingénieur civil.

    En 1834, il élabora un ambitieux plan d’alimentation en eau pour la ville de Dijon, dont l’état d’insalubrité préoccupait l’administration publique depuis trois siècles. Ce plan comprenait la dérivation d’une source distante de 12 km et la distribution de l’eau à travers un réseau d’aqueducs, de conduites et de fontaines. Le plan fut approuvé, et Darcy put en diriger l’exécution. En 1847, Dijon devint ainsi l’une des villes les mieux approvisionnées en eau potable.

    En 1856, Henry Darcy publie Les Fontaines publiques de la ville de Dijon, où il décrit la réalisation du projet et expose les principes qui doivent guider les ingénieurs dans ce type de travail. C’est aussi dans cet ouvrage qu’il fait le récit des expériences qu’il a conduites pour déterminer les lois de l’écoulement de l’eau à travers les sables. C’est là qu’il établit l’équation qui deviendra connue sous le nom de « loi de Darcy ».

    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).

    Le gradient en un point d’une fonction scalaire à plusieurs variables est un vecteur qui représente la variation locale de la fonction en ce point par rapport à la variation de ses différentes variables. Pour une fonction f à trois variables définie dans l’espace cartésien, le gradient ∇f au point (x, y, z) s’exprime ainsi :

    ∇f(x,y,z)= (∂f/∂x(x,y,z), ∂f/∂y(x,y,z), ∂f/∂z(x,y,z))

    La direction du gradient en un point correspond à celle où la fonction augmente le plus rapidement à partir de ce point, et sa grandeur témoigne de la rapidité de cette augmentation.

    L’ensemble des vecteurs gradients associés aux différents points du domaine de la fonction constitue un champ de vecteurs.

    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.

    La divergence d’un champ de vecteurs est une mesure de la dispersion du flux. Pour un champ de vecteurs V =(u, v, w) défini dans l’espace cartésien, la divergence au point (x, y, z) s’exprime ainsi :

    ∇.V(x,y,z) = ∂u/∂x(x,y,z) + ∂v/∂y(x,y,z) + ∂w/∂z(x,y,z)

    Une valeur positive de la divergence (i) correspond à une expansion du flux, alors qu’une valeur négative (ii) correspond à une compression du flux. Une divergence nulle en tout point témoigne d’un flux conservatif, en masse et en volume.

      i) V = (x, y)
    ∇.V = 2

      ii) V = (-x, -y)
    ∇.V = -2

    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é.

    Newsletter

    Recevez chaque mois une sélection d'articles

    Niveau de lecture

    Aidez-nous à évaluer le niveau de lecture de ce document.

    Si vous souhaitez expliquer votre choix, vous pouvez ajouter un commentaire (Il ne sera pas publié).

    Votre choix a été pris en compte. Merci d'avoir estimé le niveau de ce document !

    France Caron

    Chercheuse au Département de Didactique de l'Université de Montréal.
    Voir le profil

    André Garon

    Professeur à l'École Polytechnique de Montréal.
    Voir le profil