UNIVERSITE DE NEUCHATEL INSTITUT DE GEOLOGIE FACULTEDESSCIENCES CENTRE D'HYDROGEOLOGIE Une approche deterministe des distributions des temps de transit de l'eau souterraine par la theorie des reservoirs THESE Présentée à la Faculté des Sciences de l'Université de Neuchâtel pour obtenir Ie titre de DOCTEUR ES SCIENCES par David Etcheverry Soutenue le 3 juillet 2001 devant la Commission d'Examen : Professeur Pierre PERROCHET, Université de Neuchâtel, Suisse Directeur de thèse Professeur François ZWAHLEN, Université de Neuchâtel, Suisse Examinateur Professeur Piotr MALOSZEWSKI, Gesellschaft für Strahlenforschung, Allemagne Examinateur Dr. Yvan Rossier, ATE-GEOCLEAN, France Examinateur IMPRIMATUR POUR LA THESE Une approche déterministe des distributions des temps de transit de l'eau souterraine par la théorie des réservoirs de M. David Etcheverry UNIVERSITE DE NEUCHATEL FACULTEDES SCIENCES La Faculté des sciences de l'Université de Neuchâtel sur le rapport des membres du jury, MM. P. Perrochet (directeur de thèse), F. Zwahlen, P. Maloszewski (Neuherberg, D) et Y. Rossier (Meyzieu, F) autorise l'impression de la présente thèse. Neuchâtel le 8 octobre 2001 Le doyen: J.-P. Derendinger Remerciements à Pierre Perrochet sans qui ce travail n'aurait jamais vu le jour. C'est en réponse à quelques innocentes questions sur la simulation de l'âge de l'eau, pendant la pause café, qu'il m'a ouvert la porte vers un projet de recherche qui durera plusieurs années, et qui mènera finalement à ma thèse de Doctorat. Dire qu'il n'a été que mon directeur de thèse serait bien peu. Par sa complicité, par son soutien, par sa grande compétence, par son enthousiasme, et par la confiance qu'il m'a témoignée, il m'a donné les moyens d'aboutir. à Claudia, à Léo et à Maya. Suffit-il de vous remercier alors que c'est par ma faute que vous avez subi un mari grincheux, un papa absent ? Des excuses conviennent mieux. Alors : merci de tout mon cœur à vous trois et pardonnez-moi cette longue période difficile. Je vous dois beaucoup. Ne vous inquiétez pas, je ne ferai pas une deuxième thèse ! à Piotr Maloszewski de la Gesellschaft für Strahlenforschung, pour l'intérêt dont il a témoigné pour ce travail, pour l'accueil chaleureux qu'il m'a réservé à Munich, et pour avoir accepté de faire partie du jury. à Yvan Rossier de la société ATE-Geoclean, membre du jury, qui m'a introduit dans le monde de la modélisation numérique des sites contaminés. à François Zwahlen, Directeur du Centre d'Hydrogéologie de l'Université de Neuchâtel, membre du jury, pour l'intérêt qu'il a témoigné pour ce projet à Jean-Christophe Maréchal du Indo-French Center for Groundwater Research (BRGM) pour sa lecture attentive du manuscript par 40 degrés à l'ombre, et pour m'avoir demandé il y a quelques années : « pourquoi tu simules pas les âges surFeflow ? » à Michael Campana de Mexico University pour ses importants commentaires sur les solutions analytiques. à Làszlò Kiràly de VUniversité de Neuchâtel, pour son modèle hypothétique régional. 1 à Daniel J. Goode du US Geological Survey et à Jesus Carrera de la Technical University of Catalonia pour les conversations par e-mails sur la dispersion hydrodynamique de Tage de Peau. à Véronique Maître de VEcole Polytechnique Fédérale de Lausanne et à Fabien Cornaton de V Université de Neuchâtel pour leur aide précieuse. et à ma famille, à mes amis, à mes collègues, et aux autres qui ont eu la gentillesse de me donner un petit coup de main quand j'en avais besoin. 2 Table des matières 1 Introduction_______________________________________________11 1.1 L'âge de l'eau en hydrogéologie opérationnelle 11 1.2 Définitions 12 1.3 Problématique des méthodes d'investigation 13 2 Contribution du présent travail______________________________Ig 3 Théorie des réservoirs______________________________________12 3.1 Définitions 17 3.2 Distribution des âges internes 19 3.3 Distribution des temps de transit 21 3.4 Relation entre distribution des âges internes et distribution des temps de transit 22 3.5 Distribution des âges internes et distribution des espérances de vie 24 3.6 Allure des fonctions W(x) et <ï>(x) 27 3.7 Domaine de validité de la théorie des réservoirs 32 4 Dérivation analytique des distributions des temps de transit 35 4.1 Ecoulement confiné 37 4.1.1 Ecoulement confiné unidimensionnel 37 4.1.1.1 Aquifère homogène à épaisseur constante 37 4.1.1.2 Aquifère hétérogène à épaisseur constante 38 4.1.1.3 Aquifère homogène à épaisseur variable 39 4.1.2 Ecoulement tridimensionnel dans un aquifère en forme de couronne 39 4.1.3 Ecoulement horizontal confiné entre deux canaux parallèles 42 4.1.4 Ecoulement vertical vers une galerie de drainage 42 4.2 Ecoulement semi-confiné sous l'influence des précipitations 43 4.2.1 Aquifère à épaisseur variable et recharge uniforme 43 4.2.2 Paramètres de !'aquifère constants et recharge linéairement variable 49 4.2.3 Aquifère hétérogène unidimensionnel 51 4.2.4 Ecoulement radial vers un puits a pénétration totale 54 4.2.5 Influence de la dispersion longitudinale sur un écoulement semi-confiné 56 4.2.6 Aquifère semi-confiné drainé par un exutoire linéaire 58 4.2.7 Influence de la zone non-saturée 61 4.3 Conclusion 61 5 Simulation numérique des distributions des temps de transit______& 5.1 Introduction 63 5.2 Méthode classique de simulation par transport transitoire 64 5.2.1 Simulation des temps de transit 64 5.2.2 Procédure de calcul 65 5.2.3 Influence de la dispersion 66 5.3 Développement d'une méthode de simulation à l'état permanent par la théorie des réservoirs 67 5.3.1 Simulation des âges internes par l'équation de Goode 67 5.3.2 Dérivation des temps de transit par la théorie des réservoirs 69 5.3.3 Matériel informatique utilisé 69 5.3.4 Procédure de calcul 70 5.3.5 Vérification du schéma de calcul en condition de transport purement advectif 71 5.3.6 Influence du mélange numérique 72 5.3.7 Atténuation des effets du mélange par inversion des flux 74 5.4 Simulation d'un système hétérogène à cinq exutoires 78 4 5.4.1 Configuration hydrogéologique 78 5.4.2 Ecoulements 79 5.4.3 Décomposition du domaine en systèmes d'écoulement 80 5.4.4 Ages internes et temps de transit 82 5.4.4.1 Systèmes locaux 83 5.4.4.2 Systèmes intermédiaires 87 5.4.4.3 Système régional 91 5.4.5 Conclusion sur la méthode développée 93 6 CONCLUSION________________________________________________________2Z 7 ANNEXES__________________________________________________________IfiZ 5 Liste des figures Figure 1 Age de l'eau dans un système d'écoulement hypothétique à deux entrées et une sortie. L'eau entre dans le système avec un âge nul. L'âge augmente vers Pexutoire. Le volume VSï est déterminé par l'eau d'âge inférieur ou égal à x - 10 (isochrone Sx). 22 Figure 2 Espérance de vie de l'eau dans un système d'écoulement hypothétique à deux entrées et une sortie. Les conditions aux limites et les paramètres d'écoulement sont identiques à ceux de la Figure 1. Dans cet exemple, le volume vs correspond à l'eau d'espérance de vie inférieure ou égale à 10 (isochrone Sx). 25 Figure 3 Distributions d'âges internes et des temps de transit pour un cas hypothétique. 28 Figure 4 Ecoulement confiné unidimensionnel. Géométrie et conditions aux limites. 37 Figure 5 Ecoulement bidimensionnel dans un aquifère semi-circulaire confiné en forme de couronne ; a angle d'ouverture, R rayon externe, r0 rayon interne. 39 Figure 6 Distribution des temps de transit pour un écoulement contine semi-circulaire en forme de couronne. 41 Figure 7 Ecoulement confiné entre un canal et un cours d'eau. 42 Figure 8 Ecoulement vertical vers une galerie de drainage. 43 Figure 9 Ecoulement unidimensionnel semi-confiné dans un aquifère à épaisseur variable. L'infiltration est uniforme sur toute la longueur de Taquifere. 44 Figure 10 Distribution des temps de transit pour différentes geometries d'un aquifère homogène à épaisseur variable sous infiltration uniforme. 45 Figure 11 Distribution des temps de transit pour un aquifère à épaisseur et porosité constantes alimenté par une infiltration linéairement variable (i0 pour x = 0, iL à l'exutoire). 50 7 Figure 12 Ecoulement semi-confiné sans infiltration à proximité de l'exutoire. 51 Figure 13 Aquifère unidimensionnel hétérogène segmenté en N compartiments. 52 Figure 14 Distribution des temps de transit pour un aquifère unidimensionnel homogène alimenté par une infiltration non-uniforme i. L'infiltration défini trois zones d'égale longueur. L'exutoire est placé sur la gauche de la figure pour faciliter l'interprétation des résultats. L'infiltration est proportionnelle à la hauteur des rectangles gris. 53 Figure 15 Ecoulement radial semi-confiné sous l'influence des précipitations. 54 Figure 16 Champ captant dans un aquifère semi-confiné sous l'influence des précipitations. Ligne interrompue : limite des deux systèmes d'écoulement. Chaque système est décomposé en N secteurs Aj de surface et de longueur variables. Les flèches symbolisent les lignes de courant. 55 Figure 17 Influence de la dispersion sur un écoulement semi-confiné. 57 Figure 18 Ecoulement semi-confiné dans un aquifère alimenté par les précipitations et drainé par un exutoire linéaire. 58 Figure 19 Distribution des temps de transit dans un exutoire linéaire drainant un aquifère semi-confiné sous l'influence des précipitations. 60 Figure 20 Aquifère semi-confiné à zone non-saturée. 61 Figure 21 Influence des deux composantes de la dispersion sur la distribution des temps de transit. A gauche composante longitudinale, à droite composante transversale. 66 Figure 22 Algorithme de calcul du volume poreux d'âge inférieur ou égal à x dans un élément triangulaire linéaire. 70 Figure 23 Vérification du schéma de calcul par comparaison des solutions analytiques et numériques d'un écoulement circulaire. La géométrie du système est donnée en 3.1.2. 71 8 Figure 24 Influence du mélange sur les distributions obtenues par dérivation du champ permanent d'âge moyen pour un système d'écoulement hétérogène. 73 Figure 25 Age interne et espérance de vie pour un exutoire drainant deux systèmes d'âges différents. De forts gradients apparaissent pour les âges internes alors que la répartition de l'espérance de vie est plus régulière. 75 Figure 26 Atténuation du mélange numérique par inversion des flux. Le modèle exponentiel est montré à titre indicatif. L'exutoire est figuré à gauche du modèle conceptuel pour faciliter l'interprétation des résultats. 76 Figure 27 Géométrie et paramètres hydrogéologiques du système régional hypothétique. Les exutoires principaux sont numérotés de 1 à 5 ; les valeurs entre parenthèses représentent respectivement la perméabilité [m-s1] et la porosité [-] ; Les potentiels sont imposés égaux à l'altitude sur tous les noeuds de la limite supérieure. 78 Figure 28 Potentiels hydrauliques [m] et fonction de courant [m3*d-1]. En haut en gris : niveaux aquifères (K=IO^ m-s"1) ; flèches : zones d'exfiltration ; entre parenthèses : débit sortant approximatif [m3*d1]. 79 Figure 29 Isolement des systèmes d'écoulement de chaque exutoire par particle tracking inverse (noir) ; en gris : niveaux aquifères (K=IO"* m-s"'). 1 et 5 : systèmes locaux, 2 : système régional, 3 et 4 !systèmes intermédiaires. 81 Figure 30 Exutoire 1. a) potentiels hydrauliques ; en gris : niveau aquifere (K=IO"6 m-s"1), b) âges internes ; c) distribution des temps de transit. 84 Figure 31 Exutoire 5. a) potentiels hydrauliques ; en gris : niveau aquifere (K=IO"6 m-s"1), b) âges internes ; c) distribution des temps de transit. 86 Figure 32 Exutoire 3. a) potentiels hydrauliques ; en gris : niveau aquifere (K=IO"6 m-s'1), b) âges internes ; c) distribution des temps de transit. 88 Figure 33 Exutoire 4. a) potentiels hydrauliques ; en gris : niveau aquifere (K=IO"6 m-s"1), b) âges internes ; c) distribution des temps de transit. 90 9 Figure 34 Exutoire 2. a) potentiels hydrauliques ; en gris : niveau aquifère (K=IO"6 m-s*1), b) âges internes ; c) distribution des temps de transit. 92 Figure 35 Classes d'âges internes [a] correspondant aux 4 pics de la distribution des temps de transit : 7000, 14000, 22000 et 32000 a (trait gras) ; en pointillé : fonction de courant 93 10 1 Introduction Alors que la Terre va compter neuf milliards d'êtres humains au milieu du siècle, de nombreux pays seront bientôt en situation de « stress hydrique ». La quantité d'eau disponible, la répartition équitable de la production, mais aussi la qualité de l'eau, seront les enjeux de ce siècle pour les pays qui ne pratiqueront pas une gestion durable de leurs ressources hydrogéologiques. Marq De Villiers nous invite à réfléchir à ces problèmes dans son essai intitulé « L'eau » et paru aux éditions Actes Sud en 2000. Dans le vaste champ de connaissances multidisciplinaires que constitue l'hydrogéologie actuelle, Page de l'eau est omniprésent. 1.1 L'âge de Peau en hydrogéologie opérationnelle Ressources en eau La production d'eau potable doit être en équilibre avec la recharge de l'aquifere exploité. Une eau jeune indique un renouvellement des réserves. L'exploitation est envisageable à condition que la production ne dépasse pas les capacités de renouvellement de l'eau. En revanche, une eau ancienne indique soit un long parcours confiné depuis la zone d'alimentation, soit une faible recharge. Le premier cas est très favorable à l'exploitation durable des ressources alors que le deuxième ne l'est manifestement pas. Le monitoring de l'âge de l'eau en cours d'exploitation permet de détecter des apports d'eau récente ou ancienne, qui indiquent une modification des écoulements, et par conséquent une surexploitation de l'aquifere (voir par exemple Metcalfe et al., 1996 ; Blavoux, 1995, Blavoux & Letolle, 1995). Sites contaminés L'hydrogéologue concerné par l'étude des sites contaminés se pose souvent les questions : A quel moment la pollution va atteindre le champ de captage d'eau potable ? Avec quelle concentration ? Combien de temps reste-t'il pour mettre en place des mesures de remédiation ? A l'instant où la pollution atteindra la cible, comment va évoluer la concentration en fonction du temps ; diminuer, augmenter ou rester stable ? Après la mise en place du dispositif de 11 remediation, combien de temps faudra-t'il attendre pour que la production au captage puisse reprendre ? Mais toute Ia pollution n'arrive pas au même instant La concentration évolue en fonction du temps, commence par augmenter, puis diminue ou se stabilise, suivant que la pollution est ponctuelle dans le temps ou chronique. L'âge de l'eau fourni une estimation du temps de transit minimal d'un polluant conservatif. C'est donc le scénario le plus pessimiste de propagation de la pollution : déplacement rapide de la pollution, sans atténuation de la concentration par dégradation ou adsorption. Zones de protection La délimitation des zones de protection des captages en eau potable est essentiellement basée sur le temps de transit de l'eau depuis la zone d'alimentation (Blau, 1990). Les limites des zones sont fixées sur la base du temps nécessaire aux polluants pour se dégrader. Selon les directives fédérales suisses, le temps de transfert minimal d'un polluant, du captage à la limite « zone rapprochée-zone éloignée » est de dix jours. Stockage de déchets Qu'il s'agisse de déchets nucléaires ou de déchets spéciaux, les sites de stockage de déchets ont pour objet de ralentir suffisamment le déplacement des substances, pour permettre leur dégradation avant d'atteindre des cibles vulnérables. Le temps est le principal facteur qui influe sur les divers processus physico-chimiques ou bactériologiques qui mènent à rauto-épuration de l'eau souterraine. C'est pourquoi le dimensionnement des ouvrages de confinement est basé essentiellement sur les temps de transit des substances stockées. Tous ces aspects de l'hydrogéologie opérationnelle sont directement liés aux concepts théoriques d'âge et de distribution des temps de transit de l'eau ou d'une substance transportée par l'eau. 1.2 Définitions L'âge de l'eau souterraine est défini en hydrogéologie comme le temps qui s'est écoulé depuis son entrée dans un système d'écoulement. Il ne s'agit donc pas de l'âge absolu qui présente peu d'intérêt dans la pratique. L'âge de l'eau quand elle atteint l'exutoire est appelé temps de transit, pour signifier qu'elle s'est déplacée à travers tout le système, de la zone d'alimentation à 12 l'exutoire. L'âge de l'eau en un point dépend du chemin parcouru et de la vitesse à laquelle elle s'est écoulée. A l'exutoire, comme en tout point vers lequel les lignes d'écoulement convergent, l'âge n'est donc pas unique. Il doit être décrit par une distribution appelée « la distribution des temps de transit », dont le premier moment est « le temps de transit moyen ». Cette valeur moyenne est représentative de la répartition des âges pour des distributions symétriques et peu étalées. Mais comme nous allons le voir, les systèmes hydrogéologiques réels sont souvent caractérisés par des distributions des temps de transit multimodales et asymétriques, pour lesquelles le premier moment donne peu d'informations. 1.3 Problématique des méthodes d'investigation Méthodes de traçage L'âge de Peau souterraine est accessible sur le terrain par diverses méthodes de traçage. La méthode de datation la plus couramment employée est sans doute rapproche par les méthodes isotopiques. L'approche quantitative par modèles de mélange (en anglais lumped parameter models) est utilisée en hydrogéologie isotopique depuis des décennies. Dans cette approche, le système est assimilé à une fonction de transfert qui transforme le signal isotopique d'entrée en un signal de sortie. Toutes les propriétés du système sont agrégées dans un ou deux paramètres de la fonction. Le choix de la fonction de transfert est souvent justifié par les bons résultats qu'elle donne, et non par sa compatibilité avec la configuration hydrogéologique du système. Mais quelle valeur attacher à un temps de transit moyen obtenu avec un modèle qui ne correspond pas conceptuellement au problème à résoudre ? Si l'applicabilité de la fonction de transfert n'est pas vérifiée au préalable, rien ne permet de dire si les résultats de datation sont valides ou non. En revanche, si la fonction de transfert est une représentation satisfaisante du système étudié, alors la démarche mène à une datation effective de l'eau du système. Bien plus que le temps de transit moyen, c'est toute la distribution des âges qui est alors déterminée. Simulation numérique Les outils de simulation numérique permettent une approche déterministe des distributions des temps de transit. Une première estimation de l'âge peut être obtenue pai particle tracking. Mais cette méthode est mal adaptée au calcul des distributions des temps de transit, car elle livre avant 13 tout des résultats sous forme graphique. De plus, en milieu hétérogène, les trajectoires des particules sont très sensibles à la position du point de départ ; un faible écart dans la position de départ des particules peut donner des trajectoires tout à fait différentes. Les résultats sont donc difficilement reproductibles. Dans certaines configurations de la zone source de pollution et de la cible par rapport aux écoulements, les particules transportées par le flux purement advectif n'atteignent pas la cible, alors qu'elle l'atteignent si la dispersion transverse est prise en compte. La simulation de l'âge en conditions de transport dispersif est possible par simulation transitoire de la réponse du système à une injection instantanée de soluté. L'évolution de la concentration en tout point du système correspond par définition à la distribution des temps de transit pour une concentration d'injection unitaire. Cette approche nécessite des simulations sur de longues périodes, l'âge maximal n'étant pas connu à priori, mais fixé par l'utilisateur. Goode (1996) a présenté récemment une méthode de calcul des âges moyens en tout point du système, par résolution d'une équation de transport permanent Dans le cas d'exutoires ponctuels, la distribution des temps de transit se réduit par cette approche à une valeur unique, le temps de transit moyen. La bibliographie met en évidence à la fois une intense activité en datation des eaux souterraines, et un nombre limité de travaux sur la prise en compte de toute la distribution des âges. La majorité des études focalisent sur l'âge moyen de l'eau, sans que le caractère distribué de l'âge ne soit évoqué. En simulation numérique, le particle tracking est l'outil le plus utilisé dans Ia pratique pour simuler l'âge de l'eau, alors qu'il existe depuis longtemps d'autres méthodes qui permettent de simuler des distributions d'âge. Finalement, l'utilisation des fonctions de transfert en hydrologie isotopique sans justifier de leur applicabilité, a relégué les modèles de mélange au rang des « boites noires », alors qu'elles sont d'un point de vue mathématique rigoureusement identiques à la distribution des temps de transit. Le problème est récurrent : la distribution des âges de l'eau des systèmes réels n'est pas directement accessible par les méthodes de terrain ; l'influence de la géologie sur la distribution de l'âge est considérée comme bien documentée, alors que la grande majorité des publications ne traite que de l'âge moyen ; et finalement, les méthodes de simulation numérique sont mal connues et nécessitent des calculs transitoires. Il y a donc un besoin d'approfondissement des connaissances et de développement des méthodes d'investigation. 14 2 Contribution du présent travail Dans ce contexte, cette étude explore les potentialités de Ia théorie des réservoirs en vue d'une approche systémique et donc déterministe des distributions des temps de transit. L'objectif est de déterminer la distribution des temps de transit quand tous les paramètres hydrogéologiques du système sont connus (paramètres d'écoulement, géométrie et conditions aux limites). Ainsi, l'incertitude sur la forme des distributions des temps de transit est levée. Tout ce travail est basé sur la théorie des réservoirs, introduite en hydrogéologie principalement par Eriksson (1961, 1971) et par Bolin & Rodhe (1973). Dans la première partie de ce travail, la théorie des réservoirs est présentée d'un point de vue probabiliste bien adapté à l'étude macroscopique des âges en milieu poreux. Sa démonstration est formulée sur la base du théorème de la divergence de Gauss pour mettre en évidence la relation qui lie le champ vectoriel du flux d'eau aux âges à Texutoire. La relation qui lie la distribution des temps de transit au concept d'espérance de vie de l'eau est démontrée. La théorie des réservoirs est le fondement de tout ce travail, que ce soit pour les solutions analytiques ou pour le développement de la méthode numérique. La théorie des réservoirs est ensuite appliquée à des problèmes conceptuels simples pour constituer un catalogue de solutions analytiques à l'usage des praticiens. La distribution des temps de transit est résolue pour des écoulements permanents et des conditions de transport advectives. Certaines configurations ont pour solution des modèles de mélange bien connus, tels que le modèle exponentiel, le modèle piston, le modèle linéaire, le modèle exponentiel piston, ce qui valide leur applicabilité par une approche systémique. De nouveaux modèles sont introduits, en particulier une solution générale pour un aquifère à épaisseur linéairement variable. Le cas particulier d'un aquifère dont l'épaisseur diminue jusqu'à zéro vers Texutoire a pour solution le modèle linéaire, modèle considéré comme peu réaliste par les hydrogéologues isotopistes du fait de sa simplicité. Plusieurs solutions purement advectives montrent que le modèle exponentiel n'est pas la conséquence d'un mélange parfait de l'eau au sein de l'aquifère, ni même à l'exutoire, mais qu'il correspond à un rapport infîltrâtion/épaisseur/porosité constant sur tout le domaine. Ainsi, le modèle exponentiel est bien adapté pour des aquiferes homogènes, rechargés 15 uniformément et à épaisseur constante. Une version dispersive du modèle exponentiel est également introduite. Cette nouvelle solution montre une sensibilité à la dispersion longitudinale pour des nombres de Peclet inférieurs à 1. La liste des problèmes conceptuels résolus dans ce chapitre n'est pas exhaustive. L'approche proposée permet aux praticiens de résoudre d'autres problèmes simples, sans passer par les équations différentielles de transport de soluté. Dans la dernière partie de ce travail, une méthode de simulation numérique des distributions des temps de transit est développée, en vue de son application à des problèmes réalistes, tenant compte de la configuration géologique et des conditions d'écoulement. La distribution des temps de transit est calculée par application de la théorie des réservoirs au champ d'âge interne obtenu par résolution de l'équation de Goode à l'état permanent. Après avoir simulé les écoulements, l'âge moyen est calculé en tout point du modèle en conditions de transport advectif. L'intégration du volume d'eau en fonction des âges donne la distribution cumulée des âges. Deux dérivations par rapport au temps donnent successivement la distribution des âges internes puis la distribution des temps de transit, conformément à la théorie des réservoirs. Avec cette méthode, l'âge maximal à l'exutoire est toujours égal à l'âge maximal dans le système, alors que par les approches transitoires, c'est l'utilisateur qui fixe l'âge maximal (puisqu'il ne le connaît pas à priori). Les simulations à l'état permanent entraînent un gain de temps considérable par rapport aux approches traditionnelles transitoires. Des simulations nécessitant plusieurs minutes en transitoire sont ainsi effectuées en quelques secondes. Ce gain de temps permet d'aborder des problèmes beaucoup plus complexes, et ouvre la porte à des simulations tridimensionnelles. Le champ d'âge est obtenu sur tout le domaine, ce qui permet une interprétation systémique de la distribution des temps de transit à l'exutoire. La méthode est illustrée par la résolution d'un système d'écoulement régional hétérogène à cinq exutoires. Après isolation des systèmes d'écoulement, la distribution des temps de transit à chaque exutoire est dérivée du champ d'âge interne correspondant Les premiers développements de ce travail ont été présentés au Symposium de l'Agence Internationale de l'Energie Atomique à Vienne en mai 1999 et publiés dans Etcheverry & Perrochet (1999). La méthode de simulation numérique est présentée dans Etcheverry & Perrochet (2000). 16 3 Théorie des réservoirs La théorie des réservoirs {« reservoir theory » en anglais) a pour origine la théorie des temps de résidence (« residence time theory ») énoncée et popularisée par Danckwerts (1953, 1958) et par Zwietering (1959). La théorie des temps de résidence traite de particules qui entrent dans un système, qui transitent à travers lui, et qui en sortent. Chaque particule a un attribut que l'on appelle Tage et qui est le temps qu'elle a passé dans le système. Ce traitement est de première importance pour les réactions chimiques, puisqu'il permet de savoir combien de temps des substances sont restées en contact à l'intérieur d'un réacteur chimique. La théorie des réservoirs énoncée par Eriksson (1961 et 1971) et résumée dans un langage mathématique plus rigoureux par Bolin & Rodhe (1973), est l'application de la théorie des temps de résidence aux substances transportées dans des réservoirs naturels, qu'il s'agisse de l'océan, de l'atmosphère, d'un lac, ou dans notre cas d'un aquifère. Dans ce chapitre, la théorie des réservoirs est reformulée d'un point de vue probabiliste bien adapté à l'étude macroscopique des âges en milieu poreux. La démonstration est basée sur le théorème de la divergence de Gauss pour mettre en évidence la relation qui lie le champ vectoriel du flux aux âges à l'exutoire. Dans ce premier travail d'investigation, les écoulements sont considérés permanents et la dispersion est négligée. 3.1 Definitions Les définitions des termes relatifs au concept d'âge des eaux souterraines varient suivant les publications, qu'elles soient anglophones ou francophones. On retrouve dans la littérature anglophone « residence time », « transit time », et « exit time », alors que les auteurs francophones utilisent « âge », « temps de résidence », « temps de transit », et parfois « temps de séjour ». L'étymologie de « temps de résidence » (latin residere = rester) n'est pas appropriée à une substance qui quitte le système. Temps de sortie serait mieux adapté, mais nous allons nous limiter à l'expression «temps de transit» qui caractérise de façon non-ambiguë la durée du déplacement d'une substance entre deux points d'un domaine. Voici donc les définitions des expressions utilisées dans ce travail : 17 • Domaine d'écoulement : C'est un volume délimité à l'intérieur d'un milieu poreux dans lequel de l'eau circule. Cette définition est beaucoup plus vague que celle des systèmes d'écoulement de Toth (1963), et correspond mieux au champ de validité de la théorie des réservoirs. Un simple aquifère avec une zone d'alimentation et un exutoire, une zone délimitée au sein d'un système d'écoulement, ou encore un système d'écoulement composite à plusieurs exutoires satisfont tous trois à cette définition. • Particule d'eau souterraine : C'est un volume d'eau suffisamment petit pour être considéré comme une entité conservative (qui reste constituée des mêmes éléments au cours du temps). • Age de l'eau : C'est pour une particule d'eau, le temps qui s'est écoulé depuis son entrée dans le domaine d'écoulement. L'âge est donc positif ou nul, mais jamais négatif par définition. L'âge de Peau à son entrée dans le domaine est zéro. En conditions de transport purement advectives, conditions que l'on ne rencontre rigoureusement qu'en simulation mathématique, l'âge de l'eau souterraine en un point est une grandeur scalaire qui peut être décrite par une seule valeur. Dans les milieux poreux naturels, l'âge en un point n'est pas unique. Tout volume élémentaire représentatif est constitué d'eaux ayant circulé à différentes vitesses et ayant suivi des chemins différents dans la porosité de !'aquifère. Ceci est particulièrement vrai à l'exutoire vers lequel toutes les lignes d'écoulement convergent. L'âge de l'eau souterraine doit donc être décrit par une distribution. • Age interne : Le qualificatif « interne » sert à préciser qu'il s'agit de l'âge de l'eau qui n'a pas encore atteint l'exutoire. C'est l'âge de l'eau à l'intérieur du domaine d'écoulement, quelle que soit sa position. • Temps de transit : C'est le temps que met une particule d'eau pour traverser tout le domaine d'écoulement C'est donc l'âge de l'eau à l'instant précis où elle sort du domaine d'écoulement C'est la limite de l'âge interne quand la position tend vers Pexutoire. • Espérance de vie : C'est un concept peu utilisé dans Ia littérature. Cependant, il nous servira à résoudre des problèmes d'instabilité de calcul numérique quand des eaux d'âges très différents se rencontrent Vespérance de vie d'une particule d'eau est le temps qui lui reste à passer dans le système avant d'en sortir (Villermaux, 1985). 18 Soient T Tage d'une particule à un endroit donné, et "K le temps qui lui reste à passer dans le système. Si x, est le temps de transit de cette particule à la sortie du système, on a Eq. 1 xt=>x+k 3.2 Distribution des âges internes Considérons un domaine à travers lequel de l'eau s'écoule à l'état permanent donc sans perte ni gain de matière au cours du temps. L'âge de l'eau, x, a la dimension du temps et peut prendre toutes les valeurs d'un intervalle réel positif. C'est donc une variable continue. Les particules d'eau à 1*intérieur d'un domaine d'écoulement n'ont évidemment pas toutes le même âge. A un instant donné, certaines particules entrent dans le domaine, d'autres ont déjà séjourné plus ou moins longtemps en suivant les lignes de courant diverses, d'autres sortent du domaine avec un âge égal au temps de transit. A notre échelle d'observation, il est impossible de déterminer exactement l'âge de toutes les particules d'eau, x est donc une variable aléatoire dont la probabilité d'occurrence dépend de tout le contexte hydrogéologique. Kreiszic (1993) précise que la probabilité P qu'une variable aléatoire X prenne une valeur donnée à l'intérieur d'un échantillon doit être définie de façon univoque. Il en va de même pour la probabilité que X prenne toutes les valeurs d'un intervalle I. La probabilité que l'âge d'une particule d'eau soit inférieur ou égal à x est définie par Eq. 2 m(x) = P(Xs!T) m(x) est la fonction de répartition ou distribution cumulée des âges internes. Elle est par définition positive et strictement croissante. De cette première définition, on calcule Eq.3 P(x, (t) la densité de probabilité des temps de transit, ou plus simplement, la distribution des temps de transit, sera donnée par la relation Eq. 14 O(t) = —*-* =-------*-* dr F0 dt 21 3.4 Relation entre distribution des âges internes et distribution des temps de transit Soit q(x,y,z) le champ vectoriel du flux de particule défini sur le domaine d'écoulement, et soit la surface isochrone ST correspondant à la position de l'eau d'âge x (Figure 1). L'eau entre dans le domaine par Sjn la réunion des limites à flux entrant, et en sort par Sou, la réunion des limites à flux sortant. On défini pour chacune des limites au domaine d'écoulement une normale n(x,y,z) orientée positivement vers l'extérieur. normale aux surfaces S* Représentation schématique du flux aux entrées et sorties Figure 1 Age de l'eau dans un système d'écoulement hypothétique à deux entrées et une sortie. L'eau entre dans le système avec un âge nul. L'âge augmente vers l'exutoire. Le volume VSt est déterminé par l'eau d'âge inférieur ou égal à T - 10 (isochrone ST). L'eau entre nécessairement dans VSi, le volume d'eau d'âge inférieur à t, par toute la surface Sin car l'âge de l'eau à l'entrée dans le système est nul. Pour un écoulement permanent, la forme et la position des isochrones ST ne varie par en fonction du temps. En d'autres termes, t(x,y,z) est 22 une constante. Donc l'eau quitte Vx en traversant ST et éventuellement une partie de Sout que Ton notera S„„ . OUIt La loi de conservation de la matière permet de dire que le bilan des flux sur le sous-domaine d'écoulement V^ est égal à la différence entre le flux total entrant et le flux qui passe à travers ST et une portion de S0UI, ce qui s'écrit mathématiquement Eq. 15 JjQ'VqdV=jTqndS+/rq ndS + JJqndS Vs1 Sin Soot, S, A l'état permanent, il n'y a ni gain ni perte de matière à l'intérieur de V^, donc le terme de gauche est nul ( V • q « 0). Le premier terme de droite est égal au flux total entrant, soit Eq.16 JfTq-DdS=F0 siD grandeur négative du fait de l'orientation des normales. Le deuxième terme de droite représente l'eau qui sort du domaine d'écoulement avec un temps de transit inférieur ou égal à t, on a donc Eq. 17 JJq ndS=+F(T) Sou, Le troisième terme de droite représente l'eau qui sort du domaine V^, mais qui reste dans le domaine d'écoulement C'est la quantité d'eau qui passe en un temps dt d'un âge x à un âge x + dr tout en restant dans le domaine d'écoulement. En termes de probabilités on aura donc Eq. 18 chjTq-ndS = M0P(T(x) est la conséquence de la répartition des âges internes, qui elle même résulte du champ de vitesse q. Finalement, la fonction de transfert n'est influencée que par la répartition des paramètres hydrogéologiques dans le système et par les conditions aux limites. Il faut cependant noter que dans le cas d'aquifères à plusieurs exutoires, la differentiation de M(x) conduit à une distribution des temps de transit composite englobant tous les exutoires. Lorsque Ton s'intéresse à un exutoire particulier, il est nécessaire d'en isoler le système d'écoulement propre pour calculer la distribution des temps de transit sur la base du champ d'âge interne qui s'y rattache, comme le montre l'exemple hétérogène à cinq exutoires traité au chapitre 5.4. 3.5 Distribution des âges internes et distribution des espérances de vie La Figure 2 illustre le concept d'espérance de vie pour le système d'écoulement de la Figure 1. L'eau sort du système avec une espérance de vie nulle. L'espérance de vie augmente en amont de l'exutoire et est bien supérieure à zéro pour l'eau qui entre dans le système. 24 n normale aux surfaces Représentation schématique du flux aux entrees et sorties Figure 2 Espérance de vie de Peau dans un système d'écoulement hypothétique à deux entrées et une sortie. Les conditions aux limites et les paramètres d'écoulement sont identiques à ceux de la Figure 1. Dans cet exemple, le volume V3 correspond à l'eau d'espérance de vie inférieure ou égale à 10 (isochrone Sx). On peut définir pour l'espérance de vie les mêmes fonctions de distribution que pour l'âge inteme. M(X) est la masse d'eau d'espérance de vie inférieure ou égale à X, F(X) est le débit d'eau entrant dans le système avec une espérance de vie inférieure ou égale à X. De ces distributions cumulatives, on passe aux fonctions densité par simple differentiation en suivant le même raisonnement qu'en 3.2 et 3.4. Eq. 22 w Mn dX w E âk 1V(X) et O(X-) deviennent respectivement la distribution interne des espérance de vie et la distribution des espérances de vie à l'entrée dans le système. Toujours à l'aide du théorème de la divergence (voir 3.4), nous allons démontrer qu'il existe une relation simple tout d'abord entre ty(X) et ®{X), puis entre &(X) et <ê(t). Soit V8^ le volume d'eau en aval de l'isochrone d'espérance de vie Sx et q(x,y,z) le champ vectoriel du flux des 25 particules. Toute l'eau de Vs sort par l'exutoire Sou[ avec un débit égal au débit total F0. On a donc Eq.23 JTq ndS = +F0 Soul L'eau qui entre dans Vs avec une espérance de vie inférieure ou égale à X peut provenir en partie de la zone d'infiltration. Le débit d'eau qui entre dans le système avec une espérance de vie inférieure ou égale à X est par définition F(X), donc Eq. 24 JTq-IIdS =-F(X) Si-, quantité négative du fait de l'orientation des normales. Finalement, une partie de l'eau qui alimente Vs peut provenir d'eau ayant déjà séjourné un certain temps dans le système. Cette eau entre dans le volume Vs en traversant Sx et correspond à la perte de volume entre deux isochrones pour un intervalle de temps dX, ce qui s'écrit Eq.25 ffq.ndS = ^ JJ dX La somme de ces trois termes donne Eq.26 *-F(M-^ Examinons chaque terme des équations Eq. 26 et Eq. 20. M(t) et M(X) ne désignent à l'évidence pas la même quantité. Il suffit de constater la différence entre le champ d'âge interne et le champ d'espérance de vie pour s'en rendre compte (Figure 1 et Figure 2). Qu'en est-il de F(x) et de F(X) ? L'espérance de vie de l'eau quittant le système est nulle. Donc, d'après TEq. 1 on peut écrire pour l'exutoire Eq. 27 FW = F(xt) 26 Par définition, l'eau entre dans le système avec un âge interne nul. Donc, toujours d'après TEq. 1, l'espérance de vie de l'eau qui entre dans le système est égale au temps qu'elle mettra pour traverser le système. On a donc Eq. 28 F(X)=F(T.) et par conséquent F(t) « F(X). En combinant 1'Eq. 20 et PEq. 26 on démontre que Eq.29 dMÇr) = dMk)=dMfc] dr dX dxt M0 et F0 étant des constantes caractéristiques du système, on a aussi Eq. 30 W(%) » V(X) et 0>(t) « 0>(X) On peut expliquer ce résultat en remarquant que l'espérance de vie est équivalente à l'âge d'une eau qui s'écoulerait dans le même système mais en sens inverse. Quel que soit le sens de l'écoulement, le temps de transit et le débit de chaque tube de courant restent inchangés. Nous appliquerons ce résultat théorique en simulation numérique pour minimiser les effets du mélange quand des lignes de courant d'âges très différents convergent (voir 5.3.7). Le mélange d'eaux d'âges très différents induit de forts gradients d'âge interne, en particulier en conditions purement advectives, ce qui rend difficile le calcul de la distribution des temps de transit pour des raisons d'instabilités numériques. Sur le champ d'espérance de vie, l'interface entre les eaux d'âges différents défini des isochrones plus ou moins circulaires autour de l'exutoire sans que de forts gradients n'apparaissent. Il est ainsi possible d'éviter de sérieux problèmes d'instabilité. Le calcul de la distribution des temps de transit sur la base du champ d'espérance de vie est alors beaucoup plus stable qu'avec les âges internes. 3.6 Allure des fonctions 1P(T) et 0(x) Bien qu'il existe une infinité de cas de figures, les définitions précédentes des fonctions permettent de dégager leurs caractéristiques principales. La Figure 3 donne les différentes distributions obtenues pour un cas hypothétique. 27 tfWetdFWa0 eh dx f(0)-F(0)-0 W)-Y J U 0(t)äO 0(t) Figure 3 Distributions d'âges internes et des temps de transit pour un cas hypothétique. M(t), F(t) et les fonctions de répartition correspondantes ont pour ordonnée à l'origine zéro, car bien qu'il y ait de l'eau qui entre dans le domaine avec un âge nul, la masse représentative d'une classe d'âge infiniment petite est nulle. Ces fonctions sont également non décroissantes car cumulatives par définition. Il est important de noter que si M(x) est strictement croissante, F(x) peut en revanche ne pas croître. C'est le cas par exemple pour un aquifère dont l*exutoire est séparé de la zone d'alimentation, ou quand un aquifère comporte deux exutoires distincts. D'après TEq. 21 *P(x) est positive ou nulle. On sait également qu'elle est non-croissante et qu'elle peut aussi ne pas décroître. En effet, l'eau qui entre dans un aquifère peut soit rester dans raquifère auquel cas *P(t) ne varie pas (car il n'y a pas d'apports d'eau d'âge non-nul), soit 28 sortir partiellement par un exutoire, auquel cas V(x) diminue. En revanche, si de l'eau est sortie de l'aquifère avant d'avoir atteint Tage x,il y aura moins d'eau d'âge x que d'eau d'âge nul. Bolin & Rodhe (1973) montrent que l'ordonnée à l'origine de V(x) est une valeur remarquable. Nous venons de voir que F(x) est nul pour x = 0. Or, 1'Eq. 20 permet de calculer Eq. 31 Km^-F0 En insérant ce résultat dans TEq Eq. 32 l'inverse du temps de renouvellement de l'eau dans un aquifère. Bien que 4>(x) soit également une fonction de densité de probabilité, 1'Eq. 21 montre qu'elle peut être quelconque tout en restant positive ou nulle. Tous les cas de figure sont possibles en fonction de la complexité du contexte hydrogéologique. Nous verrons qu'il existe des fonctions simples, comme la fonction exponentielle ou la fonction Dirac, mais la simulation numérique en milieu hétérogène montre que les distributions sont souvent beaucoup plus complexes, multimodales et asymétriques. Le Tableau 1 résume les observations précédentes sous la forme d'un tableau de variation pour les différentes fonctions introduites. X-O+ 0 < x < +co T-* +00 m(x) M(x) dm(x) et dM(x) V(x) f(x) F(x) df(x) etdF(x) «M 0 0 >0 1 T0 0 0 aO aO OOOO OOOO AAA A AAAtAI 1 M0 0 0 1 Fo 0 aO Tableau 1 Tableau de variation des fonctions de répartition et de distribution. 10 on obtient F, Ii2V(T)-V(O) _. *-° M0 29 Nauman (1981) a mis en évidence la nécessité d'utiliser les moments pour décrire l'allure de ^(x) et de O(x). Les moments ^n d'une fonction A(x) sont donnés par la relation Eq. 33 H0=jVA(r)dt. o Par définition des probabilités (Eq. 4 et Eq. 12), les moments d'ordre zéro de *P(x) et O(x) sont égaux à l'unité. Les moments d'ordre 1 appliqués à *P(x) et à (x) donnent respectivement Eq. 34 =¾-/VF(T)Ch l'âge interne moyen qui est d'usage peu fréquent dans la pratique, et «o Eq. 35 X001 -J-XO(X) dt o le temps de transit moyen, grandeur de première importance. L'étalement et l'asymétrie d'une fonction A(x) * 3>(x) ou A(x) - 1P(X) autour de la moyenne sont donnés par les moments centrés GD Eq. 36 lï^fix-njAWèt o d'ordre n=l et n=2 qui sont la variance et le coefficient d'asymétrie. Ces moments ont été utilisés par Etcheverry & Perrochet (2000) pour décrire de façon objective la sensibilité de 4>(x) à la profondeur d'un niveau peu perméable. Les hydrogéologues utilisent couramment le temps de renouvellement T0 {turnovertimé). Pour un écoulement permanent on a Eq.37 T0-^l 30 Bolin & Rodhe (1973) mettent en évidence la relation qui existe entré le temps de renouvellement et le temps de transit moyen. En combinant 1'Eq. 35 à TEq. 21 on obtient Eq. 38 ^._^LjTdV Ma H0 I) En remarquant que l'intégrale est égale à l'opposé de l'intégrale des *P(t) sur t, Ia relation devient Eq.39 Îob = J^ ™(T)^ = Mt0T0 Le temps de transit moyen est donc rigoureusement égal au temps de renouvellement. Sur la base des premiers moments des distributions, Bolin & Rodhe (1973) distinguent trois cas que nous allons interpréter d'un point de vue hydrogéologique. a) Tjn, < T0n, : <ï>(t) est faible pour des eaux jeunes. Une grande partie de l'eau reste un certain temps dans le domaine d'écoulement avant d'en sortir. Ce cas correspond par exemple à un aquifère dont la zone de recharge est éloignée de l'exutoire. Une galerie drainant l'eau d'un massif correspond bien à cette configuration. On rencontre également cette situation dans la plupart des écoulements confinés. b) T « T : d'après Eq. 34 et Eq. 35, cette condition est satisfaite si A cft-0 Eq. 40 fTV(T)àz-fx(%)&z>*f%y(x) et - xl -T0—~ o En effectuant les intégrales et en simplifiant par t, cette relation mène à Eq-41 V(t)- -To^t& dr équation différentielle bien connue, qui s'applique à la désintégration radioactive ou aux réactions chimiques de premier ordre. La solution générale de cette équation est 31 Eq. 42 V(T)-C1C1' C1 étant une constante d'intégration. Nous pouvons utiliser comme condition aux limites 1'Eq. 32. Pour *P(t =0)« 1/T0, la constante d'intégration prend alors pour valeur C1^lIT0. Cette solution correspond à l'un des modèles les plus utilisés en hydrogéologie isotopique, le modèle exponentiel. En génie chimique, ce modèle s'applique à un réacteur à l'intérieur duquel un agitateur crée une uniformité parfaite dans la distribution des âges des réactants. Dans un tel système, la distribution des temps de transit est identique à la distribution des âges en tout point du réacteur. On ne peut bien sûr pas parler de mélange parfait pour un écoulement en milieu poreux. Zuber (1986) fait remarquer que le mélange a lieu au point de prélèvement, sans aucun mélange dans l'aquifère. Il cite Eriksson (1958) et Bredenkamp & Vogel (1970) qui interprètent ces distributions exponentielles des âges par une décroissance exponentielle de la perméabilité et de la porosité avec la profondeur. Ces interprétations sont cependant trop arbitraires pour être considérées comme des règles générales. Les solutions analytiques du chapitre suivant, obtenues pour des conditions de transport purement advectif, montrent qu'une distribution exponentielle des temps de transit est avant tout la conséquence d'un rapport infiltration/épaisseur/porosité uniforme. Nous verrons également qu'une simple variation linéaire d'épaisseur a une influence notable sur l'allure de la distribution. Le modèle exponentiel ne peut donc pas être appliqué à tous les systèmes, comme on a malheureusement tendance à le faire en hydrogéologie isotopique. c) x-lBt > T01- ¦ : On arrive à cette situation si l'exutoire et la zone d'alimentation sont assez proches pour autoriser de l'eau très jeune à sortir sans avoir parcouru un trop long chemin, et si la vitesse de l'eau qui reste le plus longtemps dans le domaine n'est pas trop forte par rapport à l'eau qui sort jeune. Cette configuration est typique des systèmes karstiques par exemple dans lesquels de l'eau très jeune court-circuite l'eau ancienne du massif. 3.7 Domaine de validité de la théorie des réservoirs Nauman & Buffham (1983) énoncent les principes de base de la théorie des temps de résidence de la façon suivante : 32 1) Le fluide qui s'écoule transporte des particules (atomes, molécules, etc.), ou tout autre élément de fluide, qui sont conservés durant le passage à travers le système. 2) Chaque particule entre à un moment donné et ressort de façon définitive à un autre instant. Une particule peut quitter temporairement le système et y revenir, mais aucune n'est dans le système depuis une durée infinie, et aucune ne peut rester indéfiniment dans le système. 3) Le système consiste en un volume fini bien délimité dans l'espace. On sait à chaque instant si une particule donnée se trouve dans le système ou non. Ces axiomes nécessitent quelques commentaires en vue de leur application aux eaux souterraines. Le premier axiome indique que Ia substance concernée peut être de l'eau, mais peu aussi être un traceur. Donc, toute l'approche développée ici s'applique également à un traceur transporté par l'eau, à condition de prendre en compte, le cas échéant, les phénomènes de sorption et de dégradation. Le deuxième axiome précise que toute l'eau doit être en mouvement Les eaux immobiles ne peuvent donc pas être traitées par cette approche, même si elles sont connectées à l'eau en mouvement On peut rencontrer de telles situations dans Ia nature dans des « zones mortes » du champ de vitesse, dans la porosité non connectée d'un aquifère, dans des poches sédimentaires de matériaux perméables, ou encore dans des réseaux karstiques. Le troisième axiome précise que Ton ne peut appliquer la théorie qu'à un domaine d'écoulement parfaitement connu. Sa géométrie doit être définie sans ambiguïté. Si l'on ne s'intéresse qu'à un seul puits d'un champ captant par exemple, il y aura lieu d'isoler le système d'écoulement du puits pour effectuer les calculs. Tout au long de ce travail, que ce soit pour les résolutions analytiques ou pour les simulations numériques, l'écoulement sera permanent, et les conditions de transport advectives. Il sera donc nécessaire dans le futur d'étendre les investigations à des écoulements transitoires et à des conditions de transport dispersif. 33 4 Dérivation analytique des distributions des temps de transit Les modèles de mélange (« lumped parameter models » en anglais) sont utilisés en hydrogéologie depuis des décennies, particulièrement en hydrogéologie isotopique. Dans cette approche, une fonction mathématique simule l'influence du système sur un traceur qui le traverse. La concentration d'entrée du traceur est transformée par la fonction de transfert en une fonction de sortie. Tous les paramètres de l'aquifère et de l'écoulement sont agrégés, rassemblés (« lumped ») dans un ou deux paramètres de la fonction de transfert. Le problème inverse est résolu par ajustement des données simulées sur les données de terrain. Cette méthode simple est bien adaptée aux systèmes mal connus, peu documentés. On trouvera une description détaillée de l'approche par modèles de mélange dans Amin & Campana (1996), Maloszewski & Zuber (1982, 1992,1993 et 1996), et dans Zuber (1986). D'après Maloszewski & Zuber (1992), la procédure à suivre pour développer un modèle est la suivante : 1) Le système doit être représenté sous la forme d'un modèle conceptuel, après quoi 2) une représentation mathématique du modèle est développée. 3) Le modèle est ensuite vérifié. La vérification est obtenue quand on a montré que le modèle se comporte comme il est supposé le faire, c'est à dire qu'il donne des résultats cohérents avec le problème fixé. Finalement, 4) le modèle est calibré pour ajuster des données observées sur Ie terrain. Dans l'approche par modèles de mélange, la compatibilité entre Ie modèle conceptuel et le problème posé n'est pas un critère de choix de la fonction de transfert. Les hydrogéologues justifient souvent leur choix par les bons résultats obtenus, ou en faisant référence à des travaux dans des conditions similaires où la fonction de transfert a donné de bons résultats. Havelle (1992) rappelle que "!'IAEA (1982), la US Nuclear Regulatory Commission (Silling, 1983), et le Swedish Nuclear Power Inspectorate (1986, 1987) font tous explicitement référence à la nécessité de démontrer qu'un modèle est une représentation « bonne », « correcte », ou « suffisante » de la réalité. Il est évident que l'utilisation d'une fonction de transfert, dont on n'a pas vérifié la compatibilité avec le système réel, est une approche stérile qui ne peut donner aucun résultat fiable. Par contre, la résolution du 35 problème inverse par une fonction de transfert conceptuellement équivalente au système réel mène à la détermination des paramètres hydrogéologiques de ce système. Les vérifications de l'applicabilité des modèles de mélange en hydrogéologie sont rares mais existent. Eriksson (1958) a introduit le modèle exponentiel (EM) en hydrogéologie en supposant une décroissance de la perméabilité avec la profondeur. Bredenkamp & Vogel (1970) montrent que EM est applicable quand la perméabilité et la porosité diminuent avec la profondeur. Vogel (1970) et Gelhar & Wilson (1974) démontrent également que EM est une solution pour des aquifères homogènes. Plus récemment, Beltman et al. (1996) présentent une solution pour un puits à pénétration totale qui se réduit à EM. Amin & Campana (1996) proposent de simuler les effets du mélange par une fonction gamma à trois paramètres. La fonction permet de modéliser tous les types de mélange du mélange parfait (EM) à l'advection pure (modèle piston). L'étude la plus complète et la plus rigoureuse du modèle exponentiel est sans doute celle d'Haitjema (1995). Il dérive notamment une solution pour une classe d'écoulement semi-confiné dans des aquifères stratifiés. Luther & Haitjema (1998) montrent la stabilité d'EM pour une variété de configurations hydrogéologiques. La signification hydrogéologique du modèle piston (PM) est bien décrite par Zuber (1986). Eriksson (1958) justifie intuitivement l'applicabilité du modèle linéaire (LM) par la proportionnalité entre le temps de transit et la longueur des lignes d'écoulement. Zuber (1986) affirme sans aucune discussion que «l'allure de cette fonction semble tellement peu naturelle qu'il n'est pas étonnant qu'aucune application pratique de ce modèle ne soit connue ». Le modèle exponentiel-piston (EPM) dérivé par WoIf & Resnick (1963) pour de systèmes chimiques réels, a été introduit en hydrogéologie pour des systèmes hétérogènes (Zuber, 1986). Il y a donc non seulement une pauvreté de la littérature en ce qui concerne les verifications de l'applicabilité des modèles de mélange, mais aussi et surtout un manque d'unité dans les approches qui permettent de résoudre ce problème. Dans une première publication (Etcheverry & Perrochet, 1999), nous avons montré que la théorie des réservoirs permet de résoudre des distributions des temps de transit dans un cadre de compréhension adapté au problème de l'âge des eaux souterraines. Dans le présent chapitre, de nombreux problèmes conceptuels simples sont résolus à l'état permanent pour des conditions de transport purement advectives (excepté un cas à dispersion longitudinale). L'applicabilité des principaux modèles de mélange est démontrée. De 36 nombreux nouveaux modèles sont introduits. Ce chapitre constitue un catalogue de solutions à l'usage du praticien. 4.1 Ecoulement confiné 4.1.1 Ecoulement confiné unidimensionnel La Figure 4 montre un système d'écoulement unidimensionnel confiné. Cette configuration est une bonne approximation pour un écoulement vertical vers une galerie, un écoulement à travers la zone non-saturée, ou un écoulement entre une rivière et un canal de drainage par exemple. La distribution des temps de transit de ce simple cas est évidemment du type Dirac. La démonstration analytique de cette hypothèse est loin d'être évidente, c'est pourquoi elle est développée en Annexe 1. Ce premier exemple est le seul de ce chapitre qui est basé sur le calcul des âges internes. Il illustre l'approche utilisée dans le chapitre suivant pour la simulation numérique des distributions des temps de transit sur la base du champ d'âge interne. 1 I=O I=L ' s s s : limite imperméable Figure 4 Ecoulement confiné unidimensionnel. Géométrie et conditions aux limites. 4J.Ì.1 Aquifère homogène à épaisseur constante M(t) est dérivée en Annexe 1 sur la base des isochrones pour un aquifère homogène à épaisseur constante. Deux dérivations successives mènent à la distribution des temps de transit sous la forme d'une fonction Dirac. Eq* 43 Tt 1 In '. POUriLfc STSt1 vmax \ *min i .2.2 -----^- et KAH 2tj2 Oa2R KAH avec AH«=H2-Hi le gradient hydraulique, r0 et R les rayons intérieur et extérieur de la couronne, Tn^n et T1n^ le temps de transit minimal et maximal. La Figure 6 montre l'allure de la distribution des temps de transit pour différentes valeurs du rapport r0/R. 40 1=50 Figure 6 Distribution des temps de transit pour un écoulement confiné semi-circulaire en forme de couronne. La distribution normalisée ne dépend pas de l'angle d'ouverture a. A l'exutoire, l'eau jeune est toujours prépondérante car la vitesse de l'eau diminue avec Ie rayon de la couronne. La contribution de l'eau ancienne augmente avec l'épaisseur de l'aquifère. L'étalement de la courbe reste important pour des aquifères peu épais. Le modèle exponentiel-piston pourrait être une bonne approximation des résultats obtenus pour r0/R > 0.80, mais le paramètre T| de ce modèle, défini comme le rapport du volume total sur le volume de l'écoulement exponentiel, n'a pas de signification pour ce type d'écoulement confiné. Une distribution Dirac des temps de transit est une bonne approximation pour les aquifères de rapport r0/R supérieur à 0.80. 41 4.1.3 Ecoulement horizontal confiné entre deux canaux parallèles La Figure 7 montre un écoulement confiné horizontal entre deux canaux parallèles. La différence de potentiel entre les deux canaux est constante. Figure 7 Ecoulement confiné entre un canal et un cours d'eau. Si Ton s'intéresse au temps de transit au point de prélèvement, l'âge de l'eau à l'entrée dans l'aquifère doit être fixé à zéro par définition. L'écoulement dans l'aquifère est de type piston puisque les deux canaux sont supposés parallèles et que le gradient hydraulique est constant Soit Tq le temps de transit dans l'aquifère. Le canal de droite draine l'aquifère de façon homogène sur toute sa longueur. Nous verrons en 4.2 que la distribution des temps de transit pour une recharge uniforme correspond au modèle exponentiel EM. Soit Tq le temps de transit moyen dans le canal de droite. L'Annexe 6 montre que pour un prélèvement ponctuel dans le canal exutoire, la distribution des temps de transit correspond au modèle exponentiel-piston (EPM) que l'on note ^ePM (T)- ^ solution analytique peut être obtenue par les formules de la théorie des circuits linéaires (modèles piston et exponentiel en série), mais il est plus élégant de la dériver par Ia théorie des réservoirs. Dans la plupart des cas, Tg le temps de transit dans Pexutoire peut être négligé par rapport au temps de transit dans l'aquifère, et Ia solution se réduit à c.1d(t). 4.1.4 Ecoulement vertical vers une galerie de drainage L'eau qui s'infiltre perpendiculairement à une galerie (Figure 8) donne la plus petite approximation du temps de transit entre Ia surface et la galerie. 42 lèvement Figure 8 Ecoulement vertical vers une galerie de drainage. Si l'on néglige la différence de longueur entre les lignes d'écoulement, la distribution des temps de transit entre la surface et la galerie peut être assimilée à un modèle piston de temps de transit moyen Tp. La galerie est rechargée uniformément sur toute sa longueur par l'eau du massif. La distribution des temps de transit dans la galerie est donc de type exponentiel (voir 4.2.1) de temps de transit moyen T^. Ici encore, l'Annexe 6 montre que la distribution des temps de transit entre la surface et un point de prélèvement dans la galerie correspond au modèle exponentiel-piston ^ePM(x)- On P^ négliger Tq le temps de transit dans la galerie, et la solution devient Oc_1D(x) le modèle piston, avec pour temps de transit moyen T^. 4.2 Ecoulement semi-confiné sous l'influence des précipitations 4.2.1 Aquifère à épaisseur variable et recharge uniforme Considérons un écoulement saturé dans un aquifère homogène dont l'épaisseur b(x) varie de b0 pour x = 0 à bL pour x = L (Figure 9). 43 F(x) F(X) infiltration = i -t> ^> -> -> H> Fo x=0 x=L Figure 9 Ecoulement unidimensionnel semi-confiné dans un aquifère à épaisseur variable. L'infiltration est uniforme sur toute la longueur de l'aquifère. b0 et bL peuvent être égales. La porosité et la perméabilité sont uniformes. La recharge est constante sur toute la zone d'alimentation. L'épaisseur moyenne est faible par rapport à la longueur de !'aquifère. On peut donc supposer la vitesse horizontale et constante sur une même verticale. La solution dérivée en Annexe 3 fait intervenir la fonction W de Lambert qui correspond à la solution en x de l'équation xe* « y (Corless et al., 1996). Pour -l/e sx < 0 Ia fonction prend simultanément deux valeurs négatives. Pour x a 0 il n'y a qu'une branche réelle positive. Seule cette dernière a un sens physique pour notre problème. On la notera LambertW. Elle nous permet de calculer la distribution des temps de transit pour un écoulement semi-confiné, soit LambertW Eq. 51 ^ScM0, bL) \ u ~Tio,t>t--bo\ 0(bL-b0) 1 + LambertW pL ~Doceb„ bo 44 La longueur de l'aquifère n'influe pas sur ^sc(T»^o» ^l) 0^ 'a fonction est normalisée par le flux total. La Figure 10 montre la distribution des temps de transit pour différentes geometries caractérisées par le rapport b0/bL. To'*_________________________________________________________________ 0 OJ 1 1.5 2 2.5 3 t/TO Figure 10 Distribution des temps de transit pour différentes geometries d'un aquifère homogène à épaisseur variable sous infiltration uniforme. La distribution oscille entre une fonction en escalier (pour b0 =0) et une fonction exponentielle décroissante (pour une épaisseur nulle à l'exutoire). La distribution est moins sensible à la géométrie de l'aquifère quand l'épaisseur diminue vers l'exutoire (b0/bL > 1) que quand l'épaisseur augmente vers l'exutoire. L'allure de la distribution est la conséquence du rapport vitesse de pore / longueur des lignes d'écoulement. Quand l'épaisseur augmente vers l'exutoire (b0/bL > 1) la vitesse de pore est faible à proximité de l'exutoire. L'eau qui entre dans le système loin de Texutoire circule plus vite que l'eau qui s'infiltre à proximité de l'exutoire. C'est ce qui 45 mène à la forme en escalier de la distribution des temps de transit. Au contraire, les aquifères dont l'épaisseur diminue vers l'exutoire (b0/bL <1) sont caractérisés par de fortes vitesses à l'exutoire. Ces fortes vitesses associées à lignes d'écoulement très courtes donnent à l'exutoire systématiquement plus d'eau jeune que d'eau vieille. C'est ce qui justifie la forme exponentielle de la distribution des temps de transit Considérons trois configurations typiques (traits gras sur la Figure 10). b0 e=bL: L'expansion de premier ordre en série de ¢$0(1,^,^) pour bL par rapport à b0 permet de calculer la limite de 1 'Eq. 51 pour une épaisseur constante. lim bo-bi, *0 Cette distribution correspond au modèle exponentiel (EM) utilisé depuis des décennies par les chimistes et les hydrogéologues isotopistes. Eriksson ( 1958) introduisit ce modèle en hydrogéologie en supposant une perméabilité décroissante avec la profondeur. Bredenkamp & Vogel (1970) assument une décroissance exponentielle de la perméabilité et de la porosité avec la profondeur. Vogel (1970) considère des âges décroissant exponentiellement avec la profondeur. Ces vérifications basées sur des suppositions plus ou moins arbitraires sont responsables de la grande confusion qui règne en ce qui concerne l'applicabilité du modèle EM aux systèmes hydrogéologiques. Certains auteurs affirment que le mélange parfait doit avoir lieu au sein même de l'aquifère, les autres répondent qu'il suffit que le mélange se fasse à l'exutoire. En fait, ce modèle est utilisé depuis des décennies dans la plupart des études isotopiques sans que son applicabilité aux systèmes hydrogéologiques soit bien comprise. Son utilisation est souvent justifiée par les bons résultats obtenus, mais pas par sa compatibilité avec le problème conceptuel étudié. Récemment, Haitjema (1995) a déterminé la solution analytique pour une classe d'aquifères à épaisseur constante. La distribution des temps de transit étant une fonction simple du temps de transit moyen, Luther & Haitjema (1998) en déduisent que la forme de la distribution est la conséquence d'un rapport b00/i0 constant. Pour vérifier leur hypothèse, nous allons dériver une solution plus générale pour une infiltration i(x) et une épaisseur b(x) variables, telles que le 46 rapport 9b(x)/i(x) reste constant. b(x) doit être proportionnel à i{x). Pour simplifier le calcul, les paramètres sont supposés varier linéairement avec x. On a alors Eq. 53 i(x) = rb(x)= r(b0 + -^j) r étant une constante de proportionnalité de dimension [T1]. F(x) devient Eq. 54 F(x)-r(bL~b°)x2 + rb0x .IbJLr Après avoir calculé la vitesse de pore et intégré l'inverse de x jusqu'à l'exutoire, on obtient / Eq. 55 x « - In r L2+r(b0 + bL) Ux2(^-bo)+2xboL)i dont les racines en x donnent des valeurs identiques de F(x) après substitution dans 1'Eq. 54. Après dérivation et normalisation par le flux total, on obtient ¢^(1). Ceci montre que si infiltration et épaisseur varient linéairement tout en restant proportionnelles, alors la distribution des temps de transit correspond au modèle exponentiel O80(T). Ce résultat n'est pas systématiquement vrai pour tout r = i(x)/ b(x), mais renforce la conclusion de Luther & Haitjema (1998). b0 »0: La limite de ^sc^^o»^) en ^o "® est indéterminée car le signe de -xi0 + 0bL n'est pas connu. Cependant, si l'on élimine b0 de (A. 18), la vitesse de pore devient constante et égale à Eq. 56 VW = S DLÜ Le temps de transit de x à l'exutoire est donc Eq. 57 t(x)«-^(L-x) et s'étend de zéro à bL0/i0 = 2T0. En substituant x(x) de 1'Eq. 57 dans F(x) on obtient 47 Eq. 58 FW-ff Après differentiation et normalisation par le flux total on obtient t \ Heaviside(2Tfl-T:) b,9 Eq.59 <*>SCML) =---------.!} ° l avecT0 -•£- Cette fonction en échelon correspond à la limite de TEq. 51 à b0 <=0. Elle est identique au modèle linéaire des hydrogéologues isotopistes. Zuber (1986) affirme que «l'allure de cette fonction semble tellement peu naturelle qu'il n'est pas étonnant qu'aucune application pratique de ce modèle ne soit connue ». La résolution analytique présentée ci-dessus démontre pourtant la validité de LM pour des écoulements dans des aquifères très réalistes dont l'épaisseur augmente de zéro en x = 0 à une valeur quelconque à I'exutoire. b0/bL = oo; La solution est obtenue par simple annulation de bL dans TEq. 51. Eq- 60 **&•**> ~ïb J -e LambertWl-e*0 1+LambertW / -Ti0 t\\ _eebo Cette distribution correspond au maximum d'influence de l'eau rapide infiltrée à proximité de I'exutoire. Elle est clairement distincte de 4>sc(t). Ce cas d'un écoulement à épaisseur variable résolu en conditions de transport purement advectives montre bien que le degré de mélange n'est pas responsable de la distribution des temps de transit Si c'était Ie cas, le modèle exponentiel devrait être obtenu pour b0/bL = oo où la convergence des lignes d'écoulement est plus forte à I'exutoire que pour une épaisseur constante. Mais ce n'est pas le cas. L'interprétation en termes de mélange est donc mal posée. Une distribution exponentielle des temps de transit doit être considérée comme la conséquence d'un rapport 6(x)b(xVi(x) constant sur toute la longueur de l'aquifère. 48 4.2.2 Paramètres de l'aquifère constants et recharge linéairement variable Des variations horizontales de lithologie comme dans les terrains alluviaux ou dans les propriétés de la couverture peuvent introduire des variations de taux d'infiltration. Considérons une infiltration linéairement variable sur un aquifère semi-confiné à épaisseur et porosité constantes. La solution est obtenue par la même procédure de calcul que dans l'Annexe 3. Pour une infiltration variant linéairement de i0 pour x « 0 à iL on obtient la distribution des temps de transit Eq. 61 ^8cU'o4L)- qk V / (iL+i0)eb°e + (iL-i0)e lu M Ob, Ol+io) ebo6-iL + i0 Cette distribution ne dépend pas de la longueur de l'aquifère car elle est normalisée par le flux total. En annulant iL dans 1'Eq. 61 on obtient Eq. 62 2*1 JSSU L eb°9-eM *scfao) tt Iimfcscfao.>l) ° 4TfRT^----\* eM+l iL-*0 b0e la distribution des temps de transit pour une infiltration décroissant de i0 en x « 0 à zéro à l'exutoire. La distribution des temps de transit pour une infiltration augmentant de zéro en x « 0 à iL à l'exutoire, est obtenue par expansion en série puis annulation de i0 Eq. 63 i b 2O2 ^scfaJ - Hm O80Ki0, iL) - 8 L ff , a (2b09+xiL) io—0 On vérifie par substitution que Eq. 64 limOsctx.io.iJ-fcscto 10-lL L'allure de la courbe est montrée à la Figure 11 pour différentes valeurs du rapport i0/iL 49 T0Mr) O---------1—I—i—I—[—I—I—I—r~i—i—i—i—i—I—i—i—i—i—j—i—i—i—i—i—i—r~~i—r O 0.5 1 1.5 2 2.5 3 t/TO Figure 11 Distribution des temps de transit pour un aquifère à épaisseur et porosité constantes alimenté par une infiltration linéairement variable (i0 pour X = 0, iL à l'exutoire). La distribution oscille autour de fc^x). Pour une infiltration augmentant vers l'exutoire (i0/iL si), l'eau qui entre dans le système à proximité de l'exutoire circule plus vite que l'eau qui s'infiltre à l'opposé (x = 0). C'est pourquoi il y a systématiquement plus d'eau récente à l'exutoire que d'eau ancienne. Pour les systèmes où l'infiltration décroît vers l'exutoire, l'eau qui entre à proximité de l'exutoire se déplace plus lentement que l'eau infiltrée plus en amont. Quand le rapport i0/iL dépasse une valeur d'environ 2, la proportion d'eau récente à l'exutoire diminue, 3>sc(T*io) donnant la plus grande quantité d'eau ancienne. Ces résultats obtenus en conditions de transport purement advectif montrent, comme dans l'exemple à épaisseur variable, qu'une distribution exponentielle des temps de transit est la conséquence d'un rapport b(x)8(xyi(x) constant sur tout le domaine. Le degré de mélange n'intervient pas dans ce cas purement advectif. 50 4.2.3 Aquifère hétérogène unidimensionnel La Figure 12 montre le cas simple d'un aquifère semi-confiné avec une couverture imperméable à proximité de l'exutoiie. - t / / / / / / / S///////. Figure 12 Ecoulement semi-confiné sans infiltration à proximité de l'exutoire. Ce problème est bien connu en ingénierie chimique (WoIf & Resnick, 1963). Il correspond à un modèle exponentiel (uniformément rechargé par la limite supérieure) et un modèle piston (près de l'exutoire) en série. On peut résoudre facilement ce système par la théorie des circuits linéaires en multipliant les transformées de Laplace des deux fonctions de transfert. Cependant, l'Annexe 6 montre que la même solution est obtenue par application de la théorie des réservoirs. On peut étendre cette solution à un système hétérogène à N compartiments homogènes d'épaisseur b^, porosité 6V et longueur Lj (Figure 13). Chaque compartiment est rechargé sur limite supérieure par une infiltration i. et collecte l'eau provenant des compartiments plus en amont (à l'exception du compartiment le plus à gauche). L'eau quitte chaque compartiment par la limite verticale aval. Toutes les autres limites sont imperméables. 51 Fj(X) X(T) FjW <7 <7 ^ ^ ij 1N-I 1N <7 <7 <; W 4* 6J 'N-I 'N bJ 0N-I *N ry ? / / / s / s / s / // s s f / / / ? j /// / s s / /'/ / / / / / / / / / j / s s / / /Ys —O H hi-i hi x=0 X=ZL. Figure 13 Aquifère uni di men si on nel hétérogène segmenté en N compartiments. Cette configuration n'est pas équivalente à plusieurs compartiments en série ou en parallèle car chaque compartiment est rechargé simultanément par de l'eau d'âge nul et par de l'eau d'âge non-nul. La théorie des circuits linéaires ne peut donc pas s'appliquer simplement. La solution est dérivée en Annexe 4 pour N compartiments selon le même principe que pour un aquifère homogène. On obtient la fonction de concaténation Eq. 65 *(tJb,i,L.el)-2- Oq Ll_ï3 bJflJ *>A N N U avec 2)tnHCinsT<2)tinaxil D=J+I n=j où t^,, est le temps de transit à travers un compartiment n (voir définition en Annexe 4). On vérifie par substitution que si la porosité, l'épaisseur et l'infiltration sont égales pour tous les compartiments, alors la solution se réduit au modèle exponentiel. Ceci confirme que la distribution exponentielle des temps de transit est la conséquence d'un rapport b(x)6(x]/i(x) constant. Des hétérogénéités horizontales dans la nature des terrains de couverture ou dans la perméabilité de raquifère peuvent entraîner des variations discontinues d'infiltration. L'Eq. 65 a été résolue pour un aquifère homogène à épaisseur constante alimenté par une infiltration non-uniforme (Figure 14). 52 T0O(T) 1.6 . - i gi i s B -\ I \ N s \ I I I 1 I -----------I I I I I I" I 1 I I I I I I I I I T I I I 1 I I I I I I I I I I I I 1 I 1 I I 0.5 1 I I I 1 I I I I F' I i i i i i 1.5 2 2.5 Figure 14 Distribution des temps de transit pour un aquifère unidimensionnel homogène alimenté par une infiltration non-uniforme i. L'infiltration défini trois zones d'égale longueur. L ' ex u toi re est placé sur la gauche de Ia figure pour faciliter l'interprétation des résultats. L'infiltration est proportionnelle à la hauteur des rectangles gris. La comparaison avec la solution obtenue pour une infiltration uniforme (EM), facilite l'interprétation. Les variations de la distribution des temps de transit sont discontinues comme l'infiltration. La quantité d'eau d'un certain âge à ï'exutoire est proportionnelle à l'infiltration. Quand l'infiltration dépasse l'infiltration moyenne, les âges correspondant sont sur-représentés à Texutoire par rapport au modèle exponentiel. Et inversement quand l'infiltration est inférieure à l'infiltration moyenne. Des résultats similaires seraient obtenus pour des variations discontinues d'épaisseur ou de porosité. 53 4.2.4 Ecoulement radial vers un puits a pénétration totale La distribution des temps de transit d'un écoulement radial semi-confiné vers un puits à pénétration totale (Figure 15) est dérivée en suivant le même schéma de calcul que pour un écoulement à épaisseur, porosité et infiltration constante (voir Annexe 3). infiltration = i Figure 15 Ecoulement radial semi-confiné sous l'influence des précipitations. Le flux F(r) à travers une section verticale à la distance r du puits est obtenu par la relation Eq. 66 F(r) = ai(R2-r2) avec R le rayon maximal et a l'angle d'ouverture du système d'écoulement. Le temps de transit de l'eau Eq. 67 R2 *(r) = TolnR2_r2 à une distance r du puits ne dépend pas de l'angle d'ouverture du système d'écoulement. En substituant r(x) de 1'Eq. 67 dans PEq. 66 on obtient la distribution des temps de transit Eq. 68 *sc-,M = ^r qui est identique à Osc(t). La distribution ne dépend ni de l'angle d'ouverture, ni de la longueur du système d'écoulement. Ce résultat peut être étendu à plusieurs puits dans un aquifère semi- 54 confiné à épaisseur, et porosité constante, alimenté par une recharge uniforme. La zone de capture d'un puits peut être segmentée en une infinité de secteurs de surface A, de longueur et angle d'ouverture variable, comme le montre la Figure 16. Figure 16 Champ captant dans un aquifère semi-confiné sous l'influence des précipitations. Ligne interrompue : limite des deux systèmes d'écoulement. Chaque système est décomposé en N secteurs Aj de surface et de longueur variables. Les flèches symbolisent les lignes de courant. La surface totale de Ia zone d'alimentation d'un puits est A0. La démonstration précédente montre que la distribution des temps de transit à ï'exutoire de chaque secteur est du type <ï>sc-r(T) pour une recharge uniforme. Elle ne dépend pas de la longueur, ni de Tangle d'ouverture des secteurs. On peut donc calculer la distribution des temps de transit à la sortie d'un puits en sommant les contributions de chaque secteur. On obtient à nouveau -X eT° Eq.69 -J-------a--£- A0I T0 le modèle exponentiel. Ceci montre que la distribution des temps de transit d'un puits dans un champ captant ne dépend ni de la géométrie des lignes d'écoulement, ni des dimensions du champ de captage, si l'infiltration est uniforme, et si l'épaisseur et la porosité de l'aquifère sont constantes. Ce résultat s'explique si l'on considère chaque ligne d'écoulement séparément. Une ligne d'écoulement de ce système peut être considérée comme un système unidimensionnel, 55 rechargé uniformément. La distribution des temps de transit à l'exutoice de cette ligne d'écoulement correspond donc au modèle exponentiel et ne dépend pas de la géométrie de récoulement La rencontre de toutes ces lignes d'écoulement à Texutoire donne le modèle exponentiel. 4.2.5 Influence de la dispersion longitudinale sur un écoulement semi-confiné Quelle est l'influence de la dispersion sur un écoulement semi-confiné unidimensionnel ? En ce qui concerne l'influence de la dispersion sur la distribution des temps de transit, la littérature ne fait référence qu'au modèle dispersif (voir par exemple Lenda & Zuber 1970, Zuber, 1986) qui est la solution unidimensionnelle d'une injection instantanée. Considérons un aquifère d'épaisseur et porosité constantes, rechargé sur sa limite supérieure par une infiltration uniforme. Pour simplifier la résolution du problème, la dispersion transverse et la diffusion moléculaire sont négligées. Seule la composante longitudinale de la dispersion est considérée. Le système d'écoulement peut être vu comme une infinité de tubes de courant parallèles. Chaque tube constitue un écoulement unidimensionnel de temps de transit moyen t0. La distribution des temps de transit à la sortie de chaque tube, est notée sc. L'influence de Ia dispersion devient sensible pour des systèmes intermédiaires quand les effets de Tadvection et de Sl la dispersion sont du même ordre de grandeur (Pe = 0.1 à 1). La contribution des temps de transit inférieurs à T0 augmente avec la dispersion de 76% pour les systèmes intermédiaires à 98% pour les systèmes dominés par la dispersion. Les systèmes purement advectifs drainent à l'exutoire 64% d'eau d'âge inférieur à T0. Le modèle exponentiel ^sc(t) doit être utilisé pour les systèmes dominés par l'advection. Pour les systèmes très lents, (Pe < 1), la distribution des temps de transit dispersive ¢^.^(1) ne peut être approximée ni par DM, ni par EM. 4.2.6 Aquifère semi-confiné drainé par un exutoire linéaire Luther & Haitjema (1998) donnent la solution analytique de la distribution des temps de transit pour différentes largeurs d'un exutoire linéaire. Us montrent que l'exutoire n'influe pas notablement sur la distribution des temps de transit tant que sa largeur ne dépasse pas 10% de la largeur de l'aquifère. Hs en concluent que le retard pris par l'eau dans l'exutoire peut être généralement négligé. La Figure 18 montre un aquifère semi-confiné drainé par un exutoire linéaire. Contrairement à Luther & Haitjema qui s'intéressent à l'écoulement perpendiculaire à l'axe de l'exutoire, nous allons nous intéresser au temps de transit de l'eau entre le point d'exfiltration et le point de prélèvement dans l'exutoire. Figure IS Ecoulement semi-confiné dans un aquifère alimenté par les précipitations et drainé par un exutoire linéaire. 58 L'eau est échantillonnée sur une section verticale de l'exutoire à la position (X,Y). La distribution des temps de transit de l'eau dans cet échantillon est dérivée en Annexe 7. On obtient — ^. jit ja 1O ~~ 1O avec Tgq et T™ les temps de transit moyens de l'aquirere et de l'exutoire respectivement (voir définition en Annexe 7). On peut vérifier ce résultat en assimilant le système à deux systèmes exponentiels en série. En effet I'aquifère est rechargé uniformément sur toute sa limite supérieure. Nous avons vu qu'à cette configuration correspond une distribution exponentielle des temps de transit (voir 4.2.1). De la même façon, l'exutoire est un système unidimensionnel rechargé uniformément sur toute sa longueur du fait de la géométrie de I'aquifère. La distribution des temps de transit de l'exutoire est donc également de type exponentiel. On peut alors appliquer la théorie des circuits linéaires (Eq. A.39) à Vl(t) - (j/T0"l)e-t''17 et à 10, la distribution est pratiquement confondue à «!^(t) . L'influence de Texutoire commence à être visible pour des rapports b9X/S inférieurs à 10. Si Ton considère en première approximation que l'exutoire a une section rectangulaire égale au produit de la hauteur par la largeur (xb), on obtient bÖX/S«9X/x. La porosité étant de l'ordre de grandeur de 0.1, on en déduit que l'influence de la dispersion est visible quand la largeur de Pexutoire est environ 100 fois inférieure à celle de Taquifère. On peut retrouver cette situation pour certains systèmes de très petite taille. Mais dans la plupart des cas, cette solution exacte montre que le temps de transit dans Texutoire peut être négligé et que <&sc(x) est une bonne approximation de la distribution des temps de transit. 60 4.2.7 Influence de la zone non-saturée L'influence de la zone non-saturée a été négligée dans les exemples qui précèdent. La Figure 20 montre un aquifère semi-confiné à géométrie radiale avec une zone non-saturée en surface. E prélèvement 0 Figure 20 Aquifère semi-confiné à zone non-saturée. L'infiltration est homogène sur toute la limite supérieure du système. Les paramètres d'écoulement de la zone non-saturée et de la zone saturée sont homogènes. Si la dispersion est négligée, la zone non-saturée peut être considérée comme un système d'écoulement unidimensionnel à distribution des temps de transit de type Dirac (voir 4.1.1.1). La zone saturée est rechargée uniformément par de l'eau d'âge égal au temps de transit moyen de la zone non- saturée. La distribution des temps de transit est donc une fonction exponentielle retardée. La solution de ce problème est donnée en Annexe 6 et correspond au modèle exponentiel-piston ^epmW- 4.3 Conclusion De nombreux modèles conceptuels sont résolus analytiquement pour la distribution des temps de transit sur la base des paramètres de l'aquifère et de l'écoulement. Certaines configurations hydrogéologiques ont pour solution des modèles de mélange utilisés en hydrogéologie isotopique. Plusieurs nouvelles solutions sont introduites. 61 Le modèle piston (PM) est la solution pour un écoulement unidimensionnel confiné. Le modèle exponentiel-piston (EPM) est la solution pour 1) deux canaux parallèles communiquant par Vintermédiaire d'un aquifère confiné, 2) une galerie de drainage alimentée uniformément sur toute sa longueur, 3) un aquifère semi-confiné avec zone non-saturée, et pour 4) un aquifère semi-confiné sans zone non-saturée, mais avec une lacune d'infiltration à proximité de l'exutoire. Dans les cas 1) et 2), la partie exponentielle de la distribution des temps de transit peut être négligée, et la solution se réduit au modèle piston. Une solution générale est dérivée pour une classe d'écoulements semi-confinés dans des aquifères à épaisseur variable. La distribution des temps de transit fait intervenir la fonction W de Lambert. Le cas particulier d'une épaisseur constante se réduit au modèle exponentiel (EM). Pour les systèmes où l'épaisseur diminue jusqu'à zéro à Fexutoire, la distribution des temps de transit se simplifie et on obtient le modèle linéaire (LM) qui est considéré comme peu réaliste par les hydrogéologues isotopiste. Une solution générale est développée pour des systèmes d'écoulement à épaisseur constante mais alimentés par une infiltration linéairement variable. Un cas à hétérogénéité discontinue est également résolu. EM est la solution exacte pour un écoulement radiai semi-confiné vers un puits à pénétration totale. Une version dispersive du modèle exponentiel prenant en compte la dispersivité longitudinale est introduite. Finalement, la distribution des temps de transit est calculée en un point d'un exutoire drainant un aquifère semi-confiné. Plusieurs solutions développées ici en conditions de transport purement advectif montrent que le mélange de l'eau au sein du système ou à l'exutoire ne peuvent pas être responsables d'une distribution exponentielle des temps de transit. EM doit être considéré comme la conséquence d'une infiltration uniforme sur un aquifère semi-confiné à épaisseur et porosités constantes. Les solutions obtenues ici par résolution de problèmes conceptuels montrent que les modèles de mélange ne doivent pas être considérés comme les fonctions de transfert de boites noires, mais comme les solutions analytiques de problèmes hydrogéologiques bien définis. La méthode de résolution analytique des distribution des temps de transit développée dans ce chapitre permettra aux praticiens de dériver leurs propres solutions pour des modèles conceptuels adaptés à la configuration hydrogéologique du problème à résoudre. 62 5 Simulation numérique des distributions des temps de transit 5.1 Introduction Plusieurs méthodes sont adaptées à la simulation numérique de Tage de l'eau souterraine. La plus répandue est probablement le « particle tracking » (voir par exemple Toth, 1996 et Oliveira & Baptista, 1997). Cette méthode basée sur le champ de vitesse de pore ne donne que des informations ponctuelles sur Tage de l'eau. Il est possible de calculer la distribution des temps de transit par un classement statistique des temps d'arrivée de nombreuses trajectoires. C'est la méthode utilisée par la Nagra sur le Wellenberg pour déterminer Ia distribution des temps de transit de particules quittant les sites de stockage de déchets (Nationale Genossenschaft für die Lagerung Radioaktiver Abfälle, 1997). Mais les résultats sont très sensibles aux points de départ des trajectoires, particulièrement en milieu hétérogène comme l'ont montré Varni & Carrera (1998). Une autre méthode est basée sur la fonction de courant Le volume de chaque tube de courant est divisé par le débit qui le traverse, ce qui donne le temps de transit moyen de chaque tube. Après réarrangement, la distribution des temps de transit est obtenue pour tout le système (Deutscher Verband für Wasserwirtschaft und Kulturbau, 1995). Le calcul de la fonction de courant limite cette méthode à des domaines bidimensionnels. L'approche la plus récente et la plus efficace fait appel aux équations de transport. La distribution des âges en tout point d'un système est déterminée par la réponse du système à l'injection d'un soluté. La fonction de densité de probabilité correspond à une injection impulsionnelle (Kinzelbach, 1992 par exemple), alors que la distribution cumulée des âges est obtenue par une injection constante (Varni & Carrera, 1998 par exemple). Goode (1996) étend Ie principe à la génération de l'âge moyen en tout point du système par la résolution d'une équation de transport permanent E>es âges moyens à l'exutoire, on peut en déduire directement la distribution des temps de transit si le nombre de points représentant l'exutoire est suffisant. Toutefois, en simulation numérique, le maillage à l'exutoire ne peut pas être densifié indéfiniment pour les exutoires ponctuels comme des sources ou des cours d'eau en 63 2D. De plus, le mélange numérique, inévitable quand les lignes d'écoulement convergent, a pour effet de moyenner les âges. Ainsi, le temps de transit maximal est souvent inférieur à l'âge interne maximal. Vami & Carrera généralisent les résultats de Goode à la distribution des âges en tout point, par la résolution d'une équation de transport transitoire. Mais ces calculs transitoires sont gourmands en temps et en puissance de calcul, particulièrement pour des modèles tridimensionnels. De plus, pour les cas purement advectifs, il est nécessaire de raffiner le modèle dans le temps (pas de temps) et dans l'espace (densité du maillage) pour capturer le signal de concentration à travers le domaine d'écoulement. L'approche proposée dans ce chapitre consiste en un couplage des travaux de Goode (19%) et de la théorie des réservoirs. Le champ d'âge moyen est généré par résolution de l'équation de transport permanent, duquel on déduit la fonction cumulée des âges internes. Deux dérivations numériques conduisent successivement à la distribution des âges internes puis à la distribution des temps de transit. Des résultats que Ton obtient après des dizaines de minutes de calcul transitoire ne demandent que quelques secondes par cette approche. Cette rapidité d'exécution n'a pas pour seul avantage un gain de temps. Elle permet d'aborder des problèmes beaucoup plus complexes qu'en transitoire, et ouvre la porte à des simulations tridimensionnelles. Les premiers résultats obtenus sur des modèles conceptuels régionaux ont été publiés dans Etcheverry & Perrochet (2000). Ils laissent entrevoir la variété des distributions des temps de transit des systèmes réels. A la fin de ce chapitre, la fonctionnalité de la méthode est démontrée sur un système hypothétique régional, hétérogène, faille, à deux niveaux aquifères et cinq exutoires. 5.2 Méthode classique de simulation par transport transitoire 5.2.1 Simulation des temps de transit Danckwerts (1953) décrit la méthode expérimentale de détermination de la distribution des âges d'un traceur. Pour un traceur non réactif, on sait que l'âge du traceur est égal à l'âge de l'eau qui le transporte (Zuber, 1986). Soit c(t,x) la concentration observée au point x =(x,y,z) et au temps t en réponse à une injection instantanée de Ia forme 64 Eq. 74 c(t,x) = ô(t-t0)pourxer" avec t0 la date d'injection, et T" l'union de toutes les limites à flux entrant, et ö la fonction Dirac. c(t, x) est par définition proportionnelle au nombre de molécules qui ont atteint le point x à Tage x = t-t0. A l'exutoire, la distribution des âges en tout point $(t,x) est donc égale à c(t, x). La somme pondérée par les débits donne la distribution des temps de transit sur l'exutoire Jq*nc(t,x) Eq. 75 n--------------&(%) Fo avec T* l'union de toutes les limites à flux sortant, q le débit, n la normale à Texutoire, et F0 le débit total. Varni & Carrera (1998) donnent le même résultat pour une injection de type échelon de la forme Eq. 76 c(t,x) = Heaviside(t-t0) pourxer" La fonction Heaviside étant l'intégrale du signal Dirac, la réponse du système à une injection en échelon correspondra à l'intégrale de TEq. 75, soit à Texutoire Jq-nc(t,x) Eq. 77 ^-------~F(t) 5.2.2 Procédure de calcul Tous les codes de transport sont adaptés à la simulation des distributions des temps de transit par transport transitoire. Les deux types d'injection sont utilisés dans la pratique. Cependant, il est plus facile d'imposer un signal en échelon qu'un signal impulsionnel. Nous utiliserons dans ce qui suit un signal en échelon. Une concentration unitaire est imposée sur toutes les limites à flux entrants. La réponse à ce signal est observée sur tous les points de l'exutoire et multipliée par le débit local. Une dérivation par rapport au temps donne la distribution des temps de transit. Cette méthode a l'avantage de rendre compte des effets du mélange, mais les calculs transitoires en font une approche gourmande en puissance et en temps de calcul. 65 5.2.3 Influence de la dispersion L'influence de la dispersion est testée sur le système d'écoulement semi-circulaire de la Figure 5 pour lequel une solution analytique a été développée en 4.1.2. Ce système d'écoulement particulier aux lignes d'écoulement parallèles permet de dissocier les effets des composantes longitudinales et transversales de la dispersion. La Figure 21 montre que la dispersion longitudinale a pour effet d'étaler les extrémités de la distribution des temps de transit. L'âge maximal augmente, ainsi que la proportion d'eau jeune. La dispersion transversale induit l'effet inverse. Le mélange entre des lignes de courants adjacentes a tendance à moyenner les âges. A l'extérieur de la couronne, les eaux vieilles se mélangent avec les eaux jeunes ce qui diminue l'âge maximal. De la même façon, sur le bord intérieur de la couronne, les eaux jeunes se mélangent avec des eaux plus vieilles, ce qui a tendance a augmenter l'âge minimal. Quand la dispersion transverse augmente, la distribution des temps de transit se rétracte sur elle-même autour du temps de transit moyen. 0.003 CLp=constante (1 m) =1000 a. =Oj=0 (solution analytique) T1T I 1 [' I I I I I I l"l I 1J T1I 1T- I' I I1I eu:^constante (1 m) eu =Op=0 (solution analytique) / oty=10m cy=100m -0.003 7O.OO25 7O.OO2 L0.0015 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 x[d] x[d] r0 = 250m R=IOOOm K-IO-4HiS"1 6 = 0.2 AH=IOOm F0=SSImV1 Figure 21 Influence des deux composantes de ta dispersion sur la distribution des temps de transit A gauche composante longitudinale, à droite composante transversale. 66 Les deux composantes de la dispersion ont donc des effets opposés. L'influence de la dispersion ne devient perceptible qu'à partir de valeurs très élevées (aL « 500 m et 04. « 100 m) pour un système d'écoulement de 2350 m de longueur moyenne et de 750 m de large. Sur des systèmes réels, les effets de Ia dispersion se limiteront donc à un lissage des distributions d'âge. Dans les systèmes réels, ce sont les longueurs des lignes d'écoulement et le champ de vitesse qui influeront le plus sur la distribution des temps de transit. 5.3 Développement d'une méthode de simulation à l'état permanent par Ia théorie des réservoirs L'allure générale de la distribution des temps de transit des systèmes réels est avant tout fonction de la géométrie de l'écoulement et du champ de vitesse. Les effets de la dispersion se limitent à un lissage de la distribution des temps de transit, comme nous venons de le voir. Travailler en dispersif ne se justifie donc que pour des systèmes très lents. Dans cette section, une méthode de simulation en conditions purement advectives à l'état permanent est proposée. Le gain en puissance et en temps de calcul est considérable. La méthode est basée sur une application de la théorie des réservoirs au champ d'âge moyen généré par la méthode de Goode (1996). 5.3.1 Simulation des âges internes par l'équation de Goode Nous venons de voir que la réponse à une injection impulsionnelle est équivalente à la distribution des âges en tout point. On peut donc calculer l'âge moyen ponctuel A(x) en appliquant l'opérateur moment à la concentration, ce qui donne pour une injection au temps t«0 œ /tc(t,x)dt Eq. 78 A(x)--£---------- J"c(t,x)dt 0 Le transport de masse dans un écoulement permanent satisfait en tout point l'équation Eq. 79 e-^^ - -q- Vc(t,x) + V- 0DVc(t,x) 67 avec 9 la porosité, D le tenseur de macro-dispersion, et q le flux permanent de Darcy tel que Eq. 80 Vq = O Goode reprend le schéma adopté par Spalding (1958) pour résoudre l'équation d'âge. En multipliant TEq. 79 par t, et en intégrant t de zéro à l'infini, on obtient OO /OD \ / <0 Eq.81 e(te(t,x)£ -8jc(t,x) dt- -q-V Jtc(t,x)dt +V-ODV J*tc(t,x)dt 0 \0 / \0 > Du fait des conditions initiales (injection instantanée), la concentration est nulle en tout point à t « ». Au temps t = 0, la multiplication par t sur tout le domaine annule Ie terme tc(t,x), bien 00 que la concentration soit égale à l'unité au point d'injection. De plus, le terme C c(t,x) dt est constant sur tout le domaine, ce qui permet d'écrire Eq. 82 e-q-v J*tc(t,x)dt Ji___________ /c(t,x)dt , 0 -V-GDV /tc(t,x)dt Jc(U) dt t 0 soit d'après 1'Eq. 78 Eq. 83 q-VA(x)-V'0DVA(x)-e En conditions de transport purement advectives, TEq. 83 devient Eq. 84 q-VA(x)«9 L'âge moyen de l'eau en tout point d'un écoulement permanent satisfait donc une équation de transport permanent avec un terme source égal à la porosité. Tous les codes de transport peuvent être modifiés par introduction de ce terme source et de la condition aux limites Eq. 85 A(x)=Opourxer 68 l'âge étant nul à rentrée dans le système. La distribution spatiale de l'âge moyen est obtenue directement sans aucune manipulation particulière. 5.3.2 Dérivation des temps de transit par la théorie des réservoirs On pourrait directement calculer la distribution des temps de transit en appliquant 1'Eq. 75 aux âges obtenus à l'exutoire. Mais l'âge moyen maximal à l'exutoire est souvent inférieur à l'âge interne maximal du fait de la dispersion transverse (voir 5.2.3). De plus, des exutoires ponctuels comme des sources ou des cours d'eau en 2D ne sont représentés que par un nombre limité de points, ce qui influe sur la définition des résultats. La théorie des réservoirs exprime la relation entre la distribution des âges internes et la distribution des temps de transit Pour un transport purement advectif, la distribution des âges en un point correspond à l'âge moyen calculé par 1'Eq. 84. Du champ d'âge moyen, on peut calculer Ia distribution cumulée des âges M(t) sur tout le domaine d'écoulement V. Celle-ci s'obtient par Eq. 86 M(x) « f8Heaviside(x - A(x)) dV v selon une procédure numérique décrite en 5.3.4. La largeur des classes d'âges est choisie par l'utilisateur. Pour chaque valeur de x, et pour chaque élément dont au moins un nœud a un âge A(x) inférieur ou égal à x, le volume d'eau d'âge inférieur ou égal à x est calculé par interpolation linéaire des âges entre les nœuds de l'élément. Deux dérivations successives mènent à la distribution des âges internes puis à la distribution des temps de transit. Ceci permet de calculer $(x) en résolvant une équation de transport permanent. La méthode est applicable quelle que soit la taille de l'exutoire. Le champ d'âge interne et la distribution des temps de transit sont obtenus pratiquement simultanément. Le temps de transit maximal est toujours égal à Tage interne maximal. 5.3.3 Matériel informatique utilisé Les équations de conservation du fluide et du soluté sont résolues par le simulateur à éléments finis Feflow® (Diersch, 1996), installé sur un ordinateur personnel compatible Windows (Pentium IH, 733 MHz, 64 Mo RAM). Ce matériel informatique permet de traiter en quelques secondes des 69 problèmes permanents à plusieurs dizaines de milliers d'éléments triangulaires linéaires. La résolution des problèmes transitoires nécessite plusieurs dizaines de minutes. L'extension de la méthode à des problèmes 3D ne pose pas de difficulté conceptuelle, mais nous nous sommes limités dans ce travail à des coupes verticales 2D. 5.3.4 Procédure de calcul Après résolution du problème d'écoulement permanent, l'âge est imposé à zéro sur tous les nœuds à flux entrant. Un terme source égal à la porosité est attribué à chaque élément. L'âge moyen de l'eau tel qu'il est défini en Eq. 78 est obtenu à chaque nœud par résolution des équations Eq. 83 ou Eq. 84. Le fichier de résultats fourni par Feflow®est traité par un programme Fortran 90 pour calculer Me(x) la masse d'eau d'âge inférieur ou égal à x dans chaque élément. L'algorithme de calcul Me(x) est détaillé à la Figure 22. Quatre cas de figure se présentent suivant le nombre de nœuds d'âge inférieur à l'isoligne x. L'isoligne est considérée linéaire à l'échelle de l'élément. Pour des éléments triangulaires à trois nœuds, l'âge en b et c est calculé par interpolation linéaire. I)M6(T)=O 2) M6(T)=BS1 3) M6(T)=BS6 - 9Sabc 4) M6(T)=BS6 Figure 22 Algorithme de calcul du volume poreux d'âge inférieur ou égal à X dans un élément triangulaire linéaire. 70 La somme de toutes les surfaces élémentaires donne pour chaque isoligne t Eq. 87 ^M6(T) = M(T) Les distributions W(x) et d>(t) sont ensuite obtenues par dérivation numérique à trois points de M(x) conformément aux Eq. 10 et Eq. 21. 5.3.5 Vérification du schéma de calcul en condition de transport purement advectif Une première vérification du schéma de calcul en conditions advectives a été publiée dans Etcheverry & Perrochet (2000) pour un cas homogène semi-circulaire. La géométrie du système est donnée à la Figure 5, et la solution analytique en Annexe 2. La dispersion est minimisée. La Figure 23 montre la parfaite concordance entre solution numérique et solution analytique. Cet exemple valide la méthode de calcul de la distribution des temps de transit par application de Ia théorie des réservoirs au champ d'âge moyen résolu par TEq. 84. r0 = 250m R = IOOOm K=IO-4InS-1 8 = 0.2 AH=IOOm F0 = 38ImV1 t . = 142.8 d min t =2284.6d 0 500 1000 1500 2000 2500 x[d] Figure 23 Vérification du schéma de calcul par comparaison des solutions analytiques et numériques d'un écoulement circulaire. La géométrie du système est donnée en 4.1.2. Dans cet exemple, les lignes d'écoulement sont parallèles et aucun échange de soluté ne se produit entre les tubes de courant. De telles conditions sont à priori rares dans la nature. II est \1.\AJ4D 0.002 0.0015 0.001 0.0005 O phi (numérique) — phi (analytique) 71 plus réaliste d'envisager des systèmes où les lignes d'écoulement convergent au niveau de l'exutoiie ou à l'intérieur du système. 5.3.6 Influence du mélange numérique Dans les systèmes réels, les exutoires sont souvent des zones où des eaux d'âges différents convergent et se mélangent. La Figure 24 montre l'influence du mélange sur la distribution des temps de transit calculée sur la base du champ permanent d'âge moyen. L'aquifère a une épaisseur constante. La recharge est uniforme. Le système est divisé en trois segments d'égale longueur. 72 T0O(T) 2 1.5- 1-: OJt Ot -0.5 t -It -1.5 0(T) 0.=0^=0.01 m 0(x) analytique 4>(x) numérique TTTTTTTT tttttttttttttttttttttttt 01L=0T= =1 m *(t) analytique 4>(x) numérique I I I J [1T' I I I I I M I I I I I I I I I I I MIfIIII 4 5 <Ï>(t) 0-=0-=0.1 m 0(t) analytique 4(X.) obtenue numériquement sur la base des espérance de vie. T0-Ct(T)1T0O(X) T0-«(tX T0-G(A) t/ro,?yro m i |'i i -i i i i i i i T/ro,x/ro = a-p = 30 m) M0 = 0.3-IO6 m3 F0 = HMmV1 T0 = 633.1 a wstt|,,oï| 2000 2500 3000 "T-----1------1-----1-----1-----1-----1-----1------1-----T- 3500 4000 4500 NN 15000 T[a] 5000 10000 Figure 30 Exutoire 1. a) potentiels hydrauliques ; en gris : niveau aquifère (K=IO"6 ms"1), b) âges internes ; c) distribution des temps de transit Le système d'écoulement de Texutoire 5 (Figure 31) est recoupé par deux niveaux aquifères de perméabilité égale à 10"6 nrs"1 (en gris sur la figure). Le niveau peu perméable situé sous et à droite de l'exutoire a une perméabilité de 10"8 m-s*1. L'extrémité droite du système a une perméabilité plus faible de 10"7 m-s"1. Le débit total est nettement supérieur à celui du système 1 84 du fait du volume des niveaux perméables. Ceci explique que le temps de transit moyen soit inférieur à celui du système 1 bien que les volumes soient très différents. Le niveau aquifère à gauche de l'exutoire draine 70 % du débit total. L'âge de l'eau de cette partie du système est globalement inférieur à 500 a. C'est donc ce niveau aquifère qui est responsable de la forme de la distribution des temps de transit entre 0 et 500 a. L'influence de la partie droite du système n'est visible que par l'étalement de Ia distribution des temps de transit. L'âge interne maximal est du même ordre de grandeur que pour le système de l'exutoire 1. La forme générale de la courbe est parfaitement rendue par l'approche transitoire, avec toujours un étalement du fait du mélange. 85 4Kt) [a-1] 0.009- permauent(a» = 10m, a™-= 1 m) transitoire (ai = a™. = 30 m) M0=LO-IO6In3 F0 »5 355.7 mV1 T0= 185.9 a 1000 I]III T-J-I- I" 1500 2000 ¦' I ' 2500 3000 3500 4000 4500 5000 10000 15000 T [a] Figure 31 Exutoire 5. a) potentiels hydrauliques ; en gris : niveau aquifère (K=IO"6 nvs'1), b) âges internes ; c) distribution des temps de transit. 86 5.4.4.2 Systèmes intermédiaires L'eau des deux systèmes intermédiaires (3 et 4) s'écoule jusqu'à environ 4000 m de profondeur. Leur zone d'alimentation est discontinue. La majeure partie de l'infiltration se produit aux abords de l'exutoire, mais une faible proportion de l'eau provient de zones éloignées par l'intermédiaire d'un mince filet d'eau qui circule dans le niveau aquifère. La Figure 32 montre les résultats obtenus pour le système d'écoulement de l'exutoire 3. La quasi totalité du système a une perméabilité de 10"8 m-s"1. L'eau qui circule en profondeur provient d'un fin tube de courant qui s'étend vers la droite à la faveur de niveaux perméables de perméabilité 10"6 et 10"7 m-s"1 (voir Figure 27 et Figure 29). Le débit total du système (138.9 m^d1) est faible par rapport à celui des systèmes locaux du fait des faibles perméabilités. Le système d'écoulement de Texutoire 3 draine trois sous-systèmes d'âges très différents, ce qui crée de forts gradients d'âge à l'exutoire. La distribution des temps de transit défini deux pics. Le deuxième pic n'apparaît pas sur les résultats transitoires puisque la concentration à l'exutoire est déjà égale à 1 quand le panache de soluté de l'eau profonde arrive en surface. Le premier pic (entre 0 et 10 000 a) est dû à des écoulements rapides d'eau infiltrée à proximité de Pexutoire. Il représente environ 80 % du débit total. Le deuxième pic correspond à de l'eau qui quitte le système avec un âge compris entre 24 000 et 28 000 a. La Figure 32-b montre que des pertes importantes de volumes se produisent entre 24 000 et 26 000 a pour l'eau ayant circulé en profondeur. Ce deuxième pic qui représente 20 % du débit total correspond donc à la sortie de l'eau profonde. 87 c ¦ 125000 10000 120000 13000Ot [a] Figure 33 Exutoire 4. a) potentiels hydrauliques ; en gris : niveau aquifère (K=IO"6 nvs"1), b) âges internes ; c) distribution des temps de transit 90 5 A A3 Système régional L'exutoire 2 draine un système d'écoulement régional dont les limites verticales et horizontale inférieure se confondent avec celles du modèle global (Figure 34-a). La faille recoupe les aquifères médian et profond. Le champ de potentiel défini des écoulements locaux et régionaux. Les premiers correspondent à de l'eau qui s'infiltre à proximité de l'exutoire et qui circule à relativement faible profondeur. Le débit total du système reste inférieur à celui de l'exutoire 4 pour un volume poreux trois fois supérieur. Il est également intéressant de noter que le temps de transit moyen (7566 a) est comparable à celui du système d'écoulement de l'exutoire 3 (7389 a), bien que le rapport des volumes poreux soit proche de 15. Ceci met en évidence le peu d'informations véhiculées par le temps de transit moyen. La recharge locale à gauche de la faille représente 3% du débit total et la partie gauche seulement 12%. Les écoulements profonds représentent donc 85% du débit total, 59% provenant de Ia partie à gauche de la faille, 26% de la partie à droite. Ce sont donc ces écoulements profonds qui influencent principalement la distribution des temps de transit Le champ d'âge interne (Figure 34-b) montre l'influence des deux niveaux aquifères. Les isochrones sont déformées dans Ie sens de l'écoulement dans les zones à forte vitesse. L'âge interne diminue dans le sens de l'écoulement dans la faille à proximité de l'exutoire du fait du mélange. Une zone de forts gradients d'âge est visible sur la partie à droite de la faille sur le flanc gauche du système d'écoulement de l'exutoire 3. Ces forts gradients d'âge sont la conséquence d'une « zone morte » dans le champ de vitesse. L'eau arrive des deux côtés de la faille avec des âges très différents, c'est pourquoi il a été nécessaire de calculer la distribution des temps de transit sur Ia base du champ d'espérance de vie par inversion des flux (voir 5.3.7). La distribution des temps de transit (Figure 34-c) a un caractère clairement quadrimodal, l'intensité des pics décroissant avec le temps de transit. La courbe obtenue en transitoire ne retranscrit que très vaguement cette multimodalité, du fait de l'important mélange dans Ia faille. 91 a ^¾¾^^¾^ ^¾ S*mMsÊ>^ ^S^S^li ^ / / n35S£ ^^*^fe^J^N\ M \ // r O • ^¾*¾ ^**^¾¾*¾¾¾*¾¾ CD \%. c *W[a_1] 0.0003 -0.00005- pcnnanent (oh = 10 m, Op = 0 m) flux inversés transitoire (a» = ay = 30 m) M0=ISJlO6In3 F0 a 2045.4013¾*1 T0 = 7566.5 a Tma. = 88.3103 a max -i—i—i—i—i—i—i—i—i—i—i—i—i—]—i—¦—i—i—i—i—i—i—i—i—i—i—i—i—i—I-I—i—i—i—i—i—i—i—j—i—i—i— 0 10000 20000 30000 40000 50000 60000 70000 80000 90000 t [a] Figure 34 Exutoire 2. a) potentiels hydrauliques ; en gris : niveau aquifère (K=IO* ms'1), b) âges internes ; c) distribution des temps de transit Une inspection détaillée du champ d'âge interne pour les classes d'âge correspondant aux quatre pics permet d'identifier l'origine de cette multimodalité (Figure 35). L'eau qui sort du système avec un âge inférieur à 7000 a provient principalement des écoulements dans les niveaux 92 aquifères à gauche de l'exutoire. Les écoulements locaux de part et d'autre de l'exutoire y contribuent aussi, mais dans une moindre proportion. L'eau de temps de transit compris entre 7000 et 14000 a provient principalement de la partie à droite de l'exutoire par l'intermédiaire de l'aquifère médian. Le pic entre 14000 et 22000 a est la conséquence des circulations dans l'aquifère profond à droite de l'exutoire. Finalement, la classe d'âges 22000 - 32000 a correspond vraisemblablement à l'eau qui entre dans la faille dans le fond du modèle. Ce dernier pic obtenu en flux inverses pourrait aussi être la conséquence du mélange numérique à la convergence des lignes d'écoulement à l'extrémité droite du système. Figure 35 Classes d'âges internes fa] correspondant aux 4 pics de la distribution des temps de transit : 7000,14000, 22000 et 32000 a (trait gras) ; en pointillé : fonction de courant. 5.4.5 Conclusion sur la méthode développée L'exemple traité dans cette section met en évidence les avantages et inconvénients de la simulation des distributions des temps de transit par transport permanent. Tout d'abord les avantages : • Les temps de calcul sont insignifiants par rapport à l'approche transitoire. Une centaine d'itérations nécessitent environ 10 min sur un Pentium® III à 733 MHz et 64 Mo de R.A.M., alors que la simulation à l'état permanent prend 10 fois moins de temps. Ces performances permettent d'aborder des problèmes complexes. Les tests de sensibilité et autres ajustements, fastidieux en transitoire, sont effectués confortablement par l'approche permanente. • L'âge maximal n'est évidemment pas connu avant d'effectuer les simulations transitoires. Il est donc nécessaire d'effectuer plusieurs essais pour déterminer le pas de temps et le nombre d'itérations idéaux. Ce problème n'intervient pas dans l'approche développée puisque l'âge 93 maximal à l'exutoire est toujours égal à Tage interne maximal calculé par résolution de l'équation appropriée. Il est fixé par la configuration hydrogéologique et non par l'utilisateur. * Même pour de très faibles valeurs du coefficient de dispersion, le mélange numérique à Pexutoire entraîne une perte d'information en transitoire, comme le montrent les résultats obtenus pour les exutoires 1, 3, et 2 notamment. La méthode permanente est beaucoup moins sensible au mélange numérique puisque la distribution des temps de transit est dérivée de tout le champ d'âge interne. Il serait donc judicieux d'inverser les flux quand les lignes d'écoulement convergent à l'exutoire. * Le champ d'âge interne, qui est la base du calcul de la distribution des temps de transit, permet l'interprétation des résultats par une approche systémique, c'est à dire sur la base de la répartition du paramètre « âge » au sein du système. L'approche transitoire ne donne elle aucune information sur les âges internes. L'interprétation des résultats ne peut se faire que par l'observation du déplacement du front de soluté après chaque itération. Le calcul des distributions des temps de transit par transport permanent présente cependant quelques inconvénients qu'il convient de signaler : * La méthode ne permet pas de simuler les effets de Ia dispersion sur la base d'un champ permanent d'âge moyen interne. Des conditions de transport dispersif nécessiteraient de calculer la distribution des âges en chaque point du système et d'intégrer la masse correspondant à chaque classe d'âge sur tout le domaine et sur tous les âges. * Pour les systèmes composites à plusieurs exutoires, il est nécessaire d'isoler les systèmes d'écoulement étudiés, ce qui peut demander un certain temps selon la complexité des cas. * La double dérivation numérique entraîne de fortes oscillations sur la distribution des temps de transit. U est cependant possible de les lisser par différentes méthodes (dérivation sur un large intervalle, moyenne mobile ou encore ajustement d'une fonction polynomiale). Toutefois, des recherches actuellement en cours au Centre d'hydrogéologie de l'Université de Neuchâtel, indiquent que les inconvénients de ces deux derniers points pourront être résolus dans un proche futur (Comaton & Perrochet, 2001, communication personnelle). D'une part, les 94 systèmes composites pourront en effet être décomposés au moyen du champ de capture établi par écoulement inverse, et sur lequel la théorie des réservoirs s'applique. D'autre part, en combinant les champs d'âge interne et d'espérance de vie, il est possible d'obtenir les fonctions de transfert, après cumul des masses concernées, au moyen d'une loi ad hoc ne nécessitant qu'une seule dérivation. 95 6 Conclusion Cette thèse explore les potentialités de la théorie des réservoirs dans le cadre d'une approche déterministe des distributions des temps de transit. Dans la première partie, la théorie des réservoirs est reformulée sur la base du champ de vitesse en appliquant le théorème de la divergence à un volume d'eau d'âge donné. Contrairement aux démonstrations simplifiées d'Eriksson (1961, 1971) et de Bolin & Rodhe (1973), l'approche par le théorème de la divergence met en évidence la relation qui existe entre le champ de vitesse à l'intérieur du système et les âges à Pexutoire. Ainsi, l'influence des paramètres d'écoulement et de la structure géologique sur les âges à l'exutoire apparaît clairement d'importance primordiale. La formulation adoptée met en évidence l'applicabilité de la théorie des réservoirs à la simulation numérique des âges de l'eau souterraine à l'état permanent. Le concept d'espérance de vie est défini et la relation qui le lie à la distribution des temps de transit est démontrée, en vue de l'appliquer en simulation numérique. La deuxième partie constitue un catalogue de solutions analytiques. Une méthode de résolution analytique des distribution des temps de transit est introduite. Elle permet de calculer des distributions des temps de transit en conditions advectives sans passer par la résolution d'équations différentielles de transport. Les solutions sont toutes obtenues en suivant le même schéma de calcul, à savoir la détermination du débit en fonction du temps de transit sur la zone d'alimentation. De nombreuses solutions sont obtenues à l'état permanent dans des conditions de transport purement advectives. Certains problèmes conceptuels ont pour solution des modèles de mélange bien connus des hydrogéologues isotopistes, (modèle exponentiel, modèle linéaire, modèle exponentiel-piston etc.), ce qui confirme leur applicabilité par une approche déterministe. Plusieurs nouvelles solutions sont introduites, en particulier une solution générale pour un écoulement semi-confiné dans un aquifère à épaisseur linéairement variable, alimenté par une infiltration homogène. Le cas particulier d'une épaisseur constante correspond au modèle exponentiel. Pour les systèmes où l'épaisseur tend vers zéro à l'exutoire, on retrouve le modèle linéaire, peu utilisé car considéré comme peu réaliste par les hydrogéologues isotopistes. Une classe d'aquifères à infiltration ou 97 épaisseur non uniformes est résolue. La solution générale permet de retrouver le modèle exponentiel-piston pour une infiltration nulle à proximité de l'exutoire. Un écoulement radial semi-confiné vers un puits complet donne à nouveau le modèle exponentiel pour une infiltration et une épaisseur uniformes. Cette solution appliquée dans les mêmes conditions à un champ captant permet de démontrer que la distribution des temps de transit de chaque puits correspond au modèle exponentiel, quelle que soit la géométrie des écoulements. Les solutions développées en conditions de transport purement advectif montrent que le mélange de l'eau au sein du système ou à l'exutoire ne justifient pas l'utilisation du modèle exponentiel. Une distribution exponentielle des temps de transit doit être considérée comme la conséquence d'un rapport infiltration/épaisseur/porosité constant sur tout le domaine. Un cas traite de l'influence de la dispersion longitudinale sur une distribution des temps de transit exponentielle. Cette version dispersive du modèle exponentiel montre que la dispersion augmente notablement la proportion d'eau jeune à l'exutoire par rapport à des conditions de transport advectives et pour des nombres de Peclet inférieurs à 1. Ces solutions analytiques incitent à ne plus considérer les modèles de mélange comme des fonctions de transfert de type « boîtes noires », mais comme des solutions analytiques de problèmes hydrogéologiques bien définis. La troisième partie explore les possibilités d'application de la théorie des réservoirs en simulation numérique. Une méthode de calcul de la distribution des temps de transit par résolution d'un problème de transport permanent est introduite en combinant la théorie des réservoirs aux récents travaux de Goode (1996). Après résolution des écoulements, l'âge moyen en chaque point du système est calculé par résolution de l'équation de Goode à l'état permanent en conditions de transport advectif. La distribution des temps de transit à l'exutoire est directement dérivée du champ d'âge interne par application de la théorie des réservoirs. Par cette méthode, Ie temps de transit maximal est toujours égal à l'âge interne maximal. Dans les approches transitoires, l'âge maximal n'est pas accessible puisqu'il est fixé par l'utilisateur. Le mélange numérique à l'exutoire entraîne une perte d'information par l'approche transitoire, même pour de très faibles valeurs de dispersion. La méthode permanente est beaucoup moins sensible au mélange puisque la distribution des temps de transit est dérivée du champ d'âge 98 inteme sur tout le domaine d'écoulement. De plus, on démontre que les calculs effectués sur l'espérance de vie après inversion des flux, permettent de minimiser les effets de la dispersion et les aberrations numériques liées à l'interpolation des âges sous fort gradient. La distribution des temps de transit et le champ d'âge interne sont calculés simultanément, ce qui permet d'identifier les paramètres internes au système qui influent sur l'âge de l'eau à l'exutoire. La simulation à l'état permanent permet un gain de temps considérable par rapport à l'approche transitoire. Cette rapidité d'exécution n'a pas pour seul avantage un gain de temps ; elle permet d'aborder des problèmes complexes, d'effectuer des tests de sensibilité, et ouvre la porte à des simulations tridimensionnelles. L'approche numérique par Ia théorie des réservoirs permet de simuler rapidement la distribution des temps de transit sur des modèles hydrogéologiques réalistes. L'exemple de simulation en fin de chapitre démontre les performances de la méthode sur un système d'écoulement hétérogène régional à cinq exutoires. Une telle approche déterministe de l'âge des eaux souterraines, basée sur la structure géologique et sur les paramètres hydrogéologiques du système, met en évidence toute ia complexité des distributions des temps de transit des systèmes réels. La méthode numérique développée présente cependant quelques inconvénients qu'il convient de signaler. Dans le cas de systèmes composites à plusieurs exutoires, il est nécessaire d'isoler les systèmes d'écoulement Les logiciels de simulation récents permettent d'effectuer ces manipulations facilement, en superposant au modèle les trajectoires simulées par particle tracking inverse. Un autre inconvénient apparaît. Il s'agit des oscillations numériques engendrées par les deux dérivations numériques successives de la masse cumulée. Ces oscillations peuvent être lissées en augmentant la densité du maillage. Les résultats sont également améliorés en effectuant les dérivations sur plusieurs points. En dernier recours, les courbes peuvent être lissées, avant dérivation, par moyenne mobile ou par ajustement d'une fonction polinomiale par exemple. Cette approche des potentialités d'application de la théorie des réservoirs n'est évidemment pas exhaustive. Des développements ultérieurs sont nécessaires, en particulier en ce qui concerne l'approche numérique. L'extension de la méthode proposée à la troisième dimension ne pose à priori aucun problème conceptuel ou technique. En revanche, des développements théoriques 99 semblent nécessaires en vue de l'extension de la méthode aux écoulements transitoires. Les travaux de Nir & Lewis (1975) et de Niemi (1977) pourraient servir de point de départ. La prise en compte de la dispersion est également une nécessité pour de nombreux systèmes, comme en témoignent les récents travaux de Goode (1996) et de Varni & Carrera (1998). D'ores et déjà, il est théoriquement et techniquement envisageable d'appliquer la théorie des réservoirs non pas au champ d'âge moyen, mais à la distribution complète des âges en chaque point du système pour simuler des conditions dispersives de transport. Une telle approche nécessite apparemment des simulations transitoires et des temps de calcul importants. Le présent travail démontre qu'il est possible de calculer la distribution des temps de transit pour des problèmes conceptuels réalistes, que ce soit par une approche analytique ou numérique. De nombreux problèmes simples peuvent être résolus par les méthodes analytiques. Pour l'hydrogéologue isotopiste, les résultats introduits constituent un outil d'aide à la décision pour le choix de fonctions de transfert adaptées à la configuration hydrogéologique du système étudié. La simulation numérique contribue à l'étude de systèmes complexes ou à grande échelle dont la distribution des temps de transit n'est pas apprehensible par les seules méthodes isotopiques. La modélisation mathématique ne se substitue cependant pas aux méthodes de datation basées sur la prise d'échantillons d'eau sur le terrain. Les deux approches sont complémentaires. Les modèles numériques permettent une meilleure compréhension de l'âge de l'eau souterraine en simulant des conditions proches de la réalité, alors que les approches de terrain apportent les données indispensables à la validation des résultats numériques. 100 Références AMIN I.E., CAMPANA M.E, 1996. A general lumped parameter model for the interpretation of tracer data and transit time calculation in hydrologie systems. Journal of hydrology 179,1-21. BELTMAN W.H., BOESTEN J.J., VAN DER ZEE S.E., QUIST J.J., 1996. Analytical modeling of effects of application frequency on pesticide concentrations in wells. Ground Water 34(3), 470-479. BLAU R.V., 1990. Groundwater protection zones in Switzerland: state of the art. Mémoires of the 22nd Congress of IAH, Vol. XXH, 1077-1086. BLAVOUX B., 1995. Apport des techniques isotopiques à la connaissance des gisements d'eau minérale. La houille blanche 2(3), 51-58. BLAVOUX B. et LETOLLE R., 1995. Apport des techniques isotopiques à la connaissance des eaux souterraines. Géochroniques 54,12-15. BOLIN B.,RODHE H., 1973. A note on the concepts of age distribution and transit time in natural reservoirs, Tellus XXV(I), 58-62. BREDENKAMP D.B., VOGEL J.C., 1970. Study of a dolomitic aquifer with carbon-14 and tritium. In: International Atomic Energy Agency (Ed.), Proc. Isotope hydrology Symp 1970, Vienna, 349-372. CORLESS R.M., GÖNNET G.H., HARE D.E., JEFFREY DJ., KNUTH D.E., 1996. On the Lambert W Function. Advances in Computational Mathematics 5,329-359. DANCKWERTS P.V., 1953. Continuous flow systems: distribution of residence times. Chemical engineering science 2(1), 1-18. DANCKWERTS P.V., 1958. The effect of incomplete mixing on homogeneous reactions. Chemical engineering sciences 8,93-102. 101 DIERSCH HJ., 1996. Interactive, graphics-based finite dement simulation system FEFLOW for modeling groundwater flow, contaminant mass and heat transport processes. WASY GmbH, Berlin. DEUTSCHER VERBAND FÜR WASSERWIRTSCHAFT UND KULTURBAU, 1995. Speicher-Durchfluss-Modeile zur Bewertung des Stoffein- und Stoffaustrags in unterschiedlichen Grundwasser-Zirkulationssystemen, Deutscher Verband für Wasser- wirtschaft und Kulturbau e. V. (DVWK), no. 109,94 p. ERIKSSON E., 1958. The possible use of tritium for estimating groundwater storage. Tellus X(4), 472-478. ERIKSSON E., 1961. Natural reservoirs and their characteristics, Geofisica International 1(2), 27-43. ERIKSSON E., 1971. Compartment models and reservoir theory, Annual Review of Ecology and Systematics 2,67-84. ETCHEVERRY D., PERROCHET P., 1999. Reservoir theory groundwater transit-time distributions and lumped parameter models. In : Proc. Int. Symp. Isotope Techniques in Water Resources Development and Management 1999, 10-14 May 1999 Vienna. International Atomic Energy Agency (Ed.), (CD-ROM). ETCHEVERRY D., PERROCHET P., 2000. Direct simulation of transit-time distributions using the reservoir theory. Hydrogeology Journal 8,200-208. FLAVELLE P.A, 1992. A quantitative measure of model validation and its potential use for regulatory purposes. Advances in Water Resources 15,5-13. GOODE D.J., 1996. Direct simulation of groundwater age, Water Resources Research 32(2), 289-296. GELHAR L.W., WILSON J.L., 1974. Ground-water quality modeling. Ground Water 12(6), 399-408. 102 HAITJEMA H.M., 1995. On the residence time distribution in idealized groundwatersheds. J. HydroL 172,127-146. I.A.E.A., 1982. Radioactive water management glossary. IAEA-TECDOC-264. International Atomic Energy Agency, Vienna. KINZELBACH W., 1992. Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser, 2. Auflage, Oldenburg Verlag, München, 343 p. KREYSZIG E., 1993. Advanced engineering mathematics. John Wiley & Sons Inc. (Eds), 7th ed, 1271 p. KRESIC N., 1997. Quantitative solutions in hydrogeology and groundwater modeling. Bocca Raton: CRC Press LLC (Lewis Publishers), 1997,471 p. LENDA A., ZUBER A., 1970. Tracer dispersion in groundwater experiments. In Isotope hydrology 1970, International Atomic Energy Agency, 619-641. LUTHER K.H., HAITJEMA H.M., 1998. Numerical experiments on the residence time distributions of heterogeneous groundwatersheds. Journal of Hydrology 207,1-17. MALOSZEWSKI P., ZUBER A., 1982. Determining the turnover time of groundwater systems with the aid of environmental tracers, I. Models and their applicability. Journal of hydrology 57,207-231. MALOSZEWSKI P., ZUBER A., 1992. On the calibration and validation of mathematical models for the interpretation of tracer experiments in groundwater, Advances in Water Resources 15,47-62. MALOSZEWSKI P., ZUBER A., 1993. Principles and practice of calibration and validation of mathematical models for the interpretation of environmental tracer data in aquifers, Advances in Water Resources 16,173-190. MALOSZEWSKI P., ZUBER A., 1996. Lumped parameter models for interpretation of environmental tracer data. In Manual on Mathematical Models in Isotope Hydrogeology. IAEA-TECDOC 910, International Atomic Energy Agency, Vienna. 103 METCALFE R., HOOKER P.J., DARLING W.G., CRAWFORD MB., 1996. The geochemical estimation of solute residence times in some reference hydrogeological settings. Nuclear science and technology series (EUR 16844), Office for official publications of the european communities (Ed.). NAUMAN E. B., 1981. Residence time distributions and micromixing, Chem. Eng. Commun 8, 53-131. NAUMAN E.B., BUFFHAM, B.A., 1983. Mixing in continuous flow systems. John Wiley and Sons, Inc. (Eds.). NATIONALE GENOSSENSCHAFT FÜR DIE LAGERUNG RADIOAKTIVER ABFÄLLE, 1997. Geosynthese Wellenberg 1996 - Ergebnisse der Untersuchungsphasen I und II - Textband, NAGRA technischer Bericht 96-01. NEEMI A.J., 1977. Residence time distributions of variable flow processes. International journal of Applied Radiation and Isotopes 28,855-860. NIR A., LEWIS S., 1975. On tracer theory in geophysical systems in steady and non-steady state. Part I and H. Tettus XXVII (4), 372-383. OLIVEIRA A., BAFIÎSTA M., 1997. Diagnostic modeling of residence times in estuaries, Water Resources Research. 33(8), 1935-1946. SILLING S.A., 1983. Final technical position on documentation of computer codes for high-level waste management, US Nuclear Regulatory Commission, NUREG-0856. SPALDING D.B., 1958. A note on mean residence-times in steady flows of arbitrary complexity. Chemical Engineering Science 9,74-77. SWEDISH NUCLEAR POWER INSPECTORATE, 1986. INTRACOIN final report, Levels 2 and 3, model validation and uncertainty analysis. SKI 86: 2. SWEDISH NUCLEAR POWER INSPECTORATE, 1987. INTRAVAL project proposal. SKI 87:3. 104 TOTH J., 1962. A theory of groundwater motion in small drainage basins in central Alberta, J. Geophys. Res. 67,4375^387. TOTH J., 1963. A theoretical analysis of groundwater flow in small drainage basins, J. Geophys. Res. 68,4795-4812. TOTH J., 1995. Hydraulic continuity in large sedimentary basins, Hydrogeology Journal 3(4), 4- 16. TOTH J., 1996. Enhancing the safety of nuclear waste disposal by exploiting regional groundwater flow: The recharge area concept, Hydrogeology Journal 4(4), 4-25. VARNl M., CARRERA J., 1998. Simulation of groundwater age distribution. Water Resources Research 34(12), 3271-3281. VILLERMAUX J., 1985. Génie de la réaction chimique. Conception et fonctionnement des réacteurs. Technique et Documentation, Lavoisier (Ed.), 26"6 édition. VOGEL JC, 1970. Carbon-14 dating of groundwater. In Isotope hydrology 1970, International Atomic Energy Agency, Vienna, 225-240. WOLF D., RESNICK W.. Residence time distribution in real systems. Industrial and Engineering Chemical Fundamentals 2(4), pp 287-293. ZUBER A., 1986. Mathematical models for the interpretation of environmental radioisotopes in groundwater systems. In: Fritz P. & Fontes J.-Ch. (Eds), Handbook of Environmental Isotope Geochemistry, vol. 2,1-59. ZWIETERING T.N., 1959. The degree of mixing in continuous flow systems. Chemical engineering sciences 11,1-15. 105 7 Annexes Annexe 1. Ecoulement confiné unidimensionnel La Figure 4 montre un système d'écoulement unidimensionnel confiné de longueur L. Le gradient hydraulique est dh/dl. Pour un aquifère homogène, la vitesse de pore v(l) est donnée par la loi de Darcy r. At /a Kdh Eq. Al V(I) = — Le volume M(*c) d'eau d'âge inférieur ou égal à t est obtenu en multipliant la vitesse de pore par le temps et par l'épaisseur de l'aquiiere. Eq. A 2 M(T) = TbK- dl Le temps de transit maximal est obtenu quand M(x) est égal au volume total d'eau mobile dans le système, « a* r bL9 M, Eq. A 3 hm t « —-s- *= —ft M(T)-M0 |jK— F0 dl quantité égale au temps de transit moyen. En dérivant TEq. A 2 par rapport au temps et en normalisant par le volume poreux total, on obtient Eq. A 4 b 107 Eq.A6 ^c-idW=ô(t-T0) la fonction Dirac ô étant définie pour tout a et b par Eq. A 7 ô(a-b) » pour a = b 0 pour a * b et/ô(a-b)da = l Annexe 2. Ecoulement confiné semi-circulaire La géométrie du modèle conceptuel est détaillée à la Figure 5. La vitesse de Darcy Eq. A 8 q(0- — ar à une distance r du centre de la couronne est fonction du gradient et de a - r la longueur totale de la ligne d'écoulement, avec AH =H2 - H1 constant sur tout le domaine. Le débit total du système est obtenu par intégration de tous les débits élémentaires sur le rayon Eq. A 9 Fo =jq(r) ^ = MS1nJEj Le volume total est fonction des rayons interne et externe et de Tangle d'ouverture du système d'écoulement. Eq. A10 Mo=0^(R2-ro2) Le temps de transit pour chaque ligne d'écoulement est donné par le rapport de la longueur sur la vitesse Eq. A II w)-------- On déduit le temps de transit maximal et minimal des rayons internes et externes respectivement ^ A ,„ 60;2¾2 9a2R2 Eq.A12 X011n=-£etx KAH 0^ KAH 108 x(r) augmente de façon monotone avec r. Le flux cumulé F(t) à l'exutoire est donc Eq. A 13 w , r(x), ,_, KAH1 (rix) F(x)=/q(r)dr--------In-" En substituant r(t) de TEq. A 11 on obtient Eq. A14 P(^MH1n(KAH iji 2a [QaX) Finalement, on calcule la distribution des âges internes et la distribution des temps de transit de Eq.A9etdeEq.A13 Eq. A 15 V(X) et InP03M, pourO La recharge étant constante sur toute la surface supérieure, le flux en x est une simple fonction linéaire de l'infiltration. Eq. A17 F(x) = i0-x 109 On en déduit la vitesse de pore Eq. A 18 v(x) IaX 0 b0+*(b---b°)l 0 L Par intégration de l'inverse de v(x) entre x et l'exutoire on obtient Eq. A19 T(x) = Ty(L(bL - b0 +b0lnL) + (bL-b0)x-Lb0lnx) La solution en x de cette équation fait intervenir la fonction W de Lambert (Corless et al., 1996) qui est définie comme la solution en y de Eq. A 20 yeJ = x Cette équation a une infinité de solutions complexes, une seule étant réelle. Seule cette dernière que Ton note LambertW(x) a un sens physique pour notre problème. La solution de TEq. A 19 en x s'écrit Eq. A 21 x(%)~ Lb° LambertW bL-*>o V b0 En substituant x par x(t) dans TEq. A 17 et en soustrayant cette quantité du débit total, on obtient Eq. A 22 F(T)-I0L 1------*— LambertW bL-b0 Pl - Pq . 6b„ b0 Une dérivation par rapport à x donne la solution générale pour la distribution des temps de transit d'un écoulement semi-confiné à épaisseur variable 110 LambertW D1 -Dn flu ' Y la „^ Eq. A 23 scMoA) ö(bL-b0) 1+LambertW / -Ti0ibL-bn\\ DL~DQcflb0 bo /> Annexe 4. Aquifère hétérogène unidimensionnel La Figure 13 montre un aquifère dont les hétérogénéités définissent j segments. La somme de l'eau qui s'est infiltrée sur le segment j et de l'eau entrée en amont dans le système par les segments n *=1 à j -1, à un débit F0n « in -Ln, correspond à Eq. A 24 \ 11-1 / Q-I le débit à travers une section verticale x(*c). La vitesse de pore est obtenue en divisant le flux par l'épaisseur et par la porosité du segment, soit Eq. A 25 VjW Vj Pour simplifier la notation, on défini N n-J.1 / i N Eq. A 26 N J dx bß, "> ' «,(x) u'lu In 1J On J-' 2F* le temps nécessaire à l'eau pour parcourir L, la longueur totale d'un segment j, L0 étant la longueur totale de l'aquifère. Le temps de transit entre une position x et l'exutoire est obtenu en sommant les temps de parcours à travers le segment j et à travers tous les compartiments en aval. On a donc 111 i> E,.A» ^)-j^.|: bfl in li f j \ 2 F0n STi J-I f & Ì 2F0n + J*-2^ k«-' V n=l If N + 2*. De la même façon que pour le cas homogène, la substitution de Xj(t) de l'Eq. A 27 dans l'Eq. A 24 donne Eq. A 28 La concaténation pour tous les segments de la première dérivée de Fj(x) normalisée par le débit total donne la distribution des temps de transit Eq. A 29 ^T1[WMl)-V il On Vj N D-JtI bß F0 n=j+l n=j Annexe 5. Influence de la dispersion longitudinale sur un écoulement semi-confiné Un domaine d'écoulement semi-confiné, d'épaisseur et de porosité constantes, alimenté par une infiltration uniforme, est assimilé à une infinité de tubes de courants en parallèle. Chaque tube de courant peut être considéré comme un système d'écoulement unidimensionnel, de temps de transit moyen t0, drainant un volume d'eau mobile m0(t0) à un débit f0(t0). Le temps de transit moyen Eq. A 30 m M fo(to) de chaque tube de courant ne dépend pas de la dispersion, de la même façon que T0 le temps de transit moyen de tout le système, car la quantité d'eau drainée par chaque tube de courant ne dépend pas des paramètres de transport de masse. De TEq. 52 on obtient 112 eT" Eq-A 31 fo(t„) = F„— T0 .[,. iJ,cft e Ax äi /t„Pe La distribution des temps de transit à l'exutoire de chaque tube de courant pour un transport dispersif est obtenue par la solution unidimensionnelle à une injection instantanée (Lenda & Zuber, 1970 par exemple) Eq. A 32 c.^W - ~----L„ ,V e +4 tPe 113 la distribution des temps de transit pour tout le système. On peut vérifier que le premier moment <*e ^sc-dispW est ^S3I à l'unité, et que le second moment central est égal à T0. En conditions de transport purement advectives, le nombre de Peclet tend vers l'infini, et la distribution tend vers 4>sc(t) , comme on peut Ie vérifier par substitution. Annexe 6. Distribution exponentiel-piston Plusieurs cas (sections 4.1.3,4.1.4 et 4.2.7) se réduisent à un écoulement piston en série avec un écoulement exponentiel. Cette configuration correspond au modèle exponentiel-piston (EPM) des chimistes et des hydrogéologues isotopistes (Maloszewski & Zuber, 1982 ; Wolf & Resnick, 1963). Nauman (1981) donne la distribution des temps de transit à la sortie de deux réacteurs en série de distributions des temps de transit (t) à la sortie du système est obtenue par la convolution des deux fonctions T Eq. A 37 0(x) -JVp1(T -t)