Des pollutions suivies à la trace

Les eaux souterraines représentent environ 23 % des réserves en eau douce de la planète, contre 76 % pour les glaciers et seulement 0,3 % pour les eaux de surface. En France, elles couvrent près des deux tiers des besoins en eau potable et une partie des besoins de l’agriculture et de l’industrie. Mais leur surexploitation hypothèque gravement la pérennité de ces réserves, sans compter le risque corollaire d’une intrusion d’eau de mer dans les nappes littorales. Quant à la qualité des eaux, elle souffre d’une multiplicité de pollutions. Comment prévoir le transport de ces polluants dans des aquifères dont la structure est généralement très hétérogène ? La compréhension des mécanismes physico-chimiques de dispersion des polluants suppose au préalable de connaître la structure de ces aquifères et par voie de conséquence la structure des couches géologiques au sein desquelles l’eau circule. Un grand nombre de mesures de terrain sont certes disponibles. Celles-ci restent toutefois parcellaires et sont donc insuffisantes pour fournir une image détaillée de ces géométries complexes, qu’il s’agisse de strates sédimentaires ou de réseaux de fractures dans des roches. Il faut donc développer des modèles numériques, lesquels seront ensuite « calés » sur les données recueillies in situ. Chercheurs en modélisation et simulation pour les sciences de l’environnement, nous concevons des modèles de nature stochastique (probabiliste) : la stratégie consiste en effet à pallier le manque de mesures en utilisant des données de type probabiliste dont les caractéristiques (moyenne, variance, loi de distribution) sont choisies de façon réaliste au regard des mesures disponibles.

L’étape suivante consiste à modéliser la circulation de l’eau, autrement dit sa vitesse dans les aquifères étudiés et sa charge hydraulique, ce dernier paramètre étant relié à la fois à la hauteur d’eau et à sa pression : l’écoulement se fait dans le sens des fortes charges vers les faibles charges. On rencontre alors une seconde difficulté, liée aux fortes hétérogénéités du milieu. C’est en particulier le cas des zones sédimentaires, les plus perméables étant constituées d’une juxtaposition d’argiles étanches et de sables perméables, eux-mêmes plus ou moins fins ou compactés. C’est aussi le cas des zones de roches fracturées, les fractures étant d’ailleurs susceptibles de former un réseau connecté, avec des écoulements assez importants pour offrir une ressource exploitable. Il faut donc recourir à des modèles capables de prendre en compte ces hétérogénéités. Or celles-ci couvrent une large gamme d’échelles, ce qui impose de découper la zone étudiée, généralement étendue (la taille d’un aquifère), en un très grand nombre de mailles de petite taille (de plusieurs millions à l’ordre du milliard).

La circulation de l’eau est régie par deux équations de base : la loi dite de Darcy et la loi de conservation de la masse. La première exprime que la vitesse est proportionnelle au gradient de charge hydraulique, le coefficient de proportionnalité – la perméabilité – étant caractéristique du milieu géologique. La charge varie elle-même dans le temps en fonction de la porosité du milieu, de la densité de l’eau (plus ou moins salée par exemple) et de la saturation en eau (mélangée ou non à de l’air ou à un autre fluide). Quant à la seconde équation, elle exprime que la masse d’eau en circulation reste constante tant qu’il n’y a ni recharge ni pompage. En l’absence de données géologiques suffisantes – perméabilité, porosité – en tout point de la zone étudiée et d’une connaissance précise des lois de comportement – densité, saturation –, ces différents paramètres sont eux aussi représentés par des variables aléatoires. Nous sommes ainsi conduits à résoudre un système d’équations stochastiques, et ce au sein de chaque maille. Pour ce faire, notre stratégie consiste à estimer une valeur moyenne de la circulation de l’eau en nous fondant sur des simulations associées à un tirage au hasard des valeurs des paramètres (d’une simulation à plusieurs centaines). Autant dire qu’il s’agit d’opérations très gourmandes en puissance de calcul. Le nombre de mailles peut en effet atteindre plusieurs millions, voire plusieurs milliards. D’où l’utilisation des superordinateurs ou des grappes de processeurs connectés entre eux par des réseaux ultrarapides.

Une fois la circulation de l’eau modélisée, il s’agit de déterminer les concentrations d’espèces chimiques présentes. Ces espèces sont soit en solution, soit sous forme de précipité, soit encore fixées aux grains de la roche. Le problème est de déterminer la répartition spatiale et temporelle de chacune d’elles, autrement dit l’évolution des concentrations des polluants dans telle ou telle maille, l’une des questions étant de savoir si un équilibre finit par s’instaurer. Nous sommes alors confrontés au couplage de deux modèles stochastiques, celui de la circulation de l’eau et celui représentant la dispersion des polluants.

À un instant donné, seules les espèces en solution dans l’eau sont mobiles sous l’effet de la dispersion. Elles réagissent aussi entre elles et avec les espèces fixes, ce qui modifie la composition chimique de l’eau. La dispersion des polluants est régie par trois mécanismes : un mécanisme d’advection (transport par le flux d’eau), un mécanisme de dispersion cinématique (passage d’une ligne de courant à une autre, du fait des hétérogénéités essentiellement), un mécanisme de diffusion moléculaire (de type mouvement brownien). À cela viennent s’ajouter les réactions chimiques éventuelles des espèces entre elles et avec la roche. C’est notamment sur la modélisation de ces mécanismes que portent nos contributions au projet MICAS de l’Agence nationale de la recherche (ANR) et au Groupement de recherche MOMAS.

Nous avons par exemple testé notre modèle sur un cas simple dans le cadre de MICAS : une espèce chimique unique et inerte chimiquement (sans aucune réaction avec le milieu environnant). La vidéo ci-dessus visualise les résultats de simulations de la dispersion dans un aquifère plus ou moins hétérogène. Si la perméabilité du milieu était homogène et s’il n’y avait pas de diffusion moléculaire, les trajectoires seraient identiques et horizontales. Dans l’exemple de cette vidéo en revanche, le mouvement des particules reflète la structure complexe du champ de vitesse et l’effet du mécanisme de diffusion. Au bout d’un temps assez long, on observerait une stabilisation du panache de pollution.

Un autre exemple intéressant est celui, purement théorique, que nous avons traité dans le contexte du groupement de recherche MOMAS. Nous avons simulé la dispersion au cours du temps, dans un aquifère, de treize espèces chimiques réagissant entre elles. Les caractéristiques géologiques ont été fixées a priori et nous avons attribué deux valeurs possibles à la perméabilité, dans deux zones distinctes. Le modèle de circulation de l’eau et de dispersion des polluants est donc dans ce cas déterministe. La principale difficulté numérique vient alors du couplage non linéaire entre les phénomènes d’advection-dispersion-diffusion, d’un côté, et les réactions chimiques, de l’autre. Cette étude peut s’appliquer à la pollution d’un site par des déchets non radioactifs.

Cette série de résultats de simulation est réalisée sur la base d’un modèle géologique déterministe et d’un système chimique comportant 13 espèces.
De t = 0 à t = 5000, on injecte des produits par la gauche ; de t = 5000 à t = 6000 on arrête l’injection. Les images successives montrent la dispersion dans l’aquifère d’une des espèces au cours du temps. En rouge sont représentées les régions où sa concentration est la plus élevée. L’échelle de couleurs est modifiée d’une image à l’autre pour que les différences soient bien visibles.

Cette image montre la vitesse de l’eau, qui transporte les produits chimiques ; au milieu, on voit une sorte de barrière beaucoup moins perméable. Dans cette simulation, la vitesse n’évolue pas au cours du temps.
Images © MOMAS, Inria, Andra.

Reste que l’un des défis à relever à court terme, dans toutes ces recherches, est de parvenir à intégrer un nombre croissant de processus hydrogéologiques et physicochimiques, afin de nous rapprocher de situations réelles. Cela suppose de mettre l’accent sur une optimisation des algorithmes, compte tenu de la puissance de calcul et du temps pour résoudre de tels systèmes fortement non linéaires. De nouvelles pistes commencent également à être explorées, en particulier la détermination des caractéristiques hydrauliques à partir des mesures (résolution du problème inverse). Enfin, les méthodes que nous mettons au point pourront s’appliquer à l’étude du comportement d’autres fluides souterrains comme le pétrole. Elles constituent en outre un outil précieux pour les études de sécurité des futurs stockages géologiques de dioxyde de carbone.

Une première version de cet article est parue dans les Cahiers de l’Inria en partenariat avec la revue La Recherche, n°430, mai 2009.