I. — Introduction

Il est notoirement reconnu aujourd'hui que les hydrosystèmes souterrains sont, par définition, peu accessibles en dépit de nombreux efforts prometteurs venant de l'imagerie géophysique de sub-surface (par ex. Mari & Porel, 2008 ; Mari et al., 2009 ; Hinderer et al., 2009). Ces systèmes restent très sous-échantillonnés alors même que la problématique ressource en eau et gestion voudrait que l'on puisse avoir une bonne vision de la qualité et de la quantité d'eau, ainsi que leur évolution face à de multiples perturbations naturelles et anthropiques. Très vite dans l'histoire des sciences hydrologiques, les modèles analytiques, analogiques puis numériques ont constitué des substituts à l'absence d'observation afin de mieux comprendre la dynamique des hydrosystèmes. Ainsi, dès le début des années 80, apparaissent des modèles numériques distribués pour la simulation des écoulements souterrains, les efforts de recherche passant rapidement dans le domaine public pour la construction d›outils opérationnels (par ex. McDonald & Hargbaugh, 1980 ; Thiery, 1990). Ces outils servent, dans un premier temps, de schéma explicatif et interprétatif des données disponibles. Ensuite, ils deviendront des demandes formulées par les décideurs pour l'aide à la gestion des ressources en eau. Plus spécifiquement en France, le développement dans les années 90 des SAGE et SDAGE (Schéma (Directeur) d'Aménagement et Gestion des Eaux) incite les gestionnaires à développer la modélisation des grands hydrosystèmes régionaux (par ex. Amraoui et al., 1999 ; Mougin et al., 2008). Quelques études plus locales existent, majoritairement construites comme le raffinement spécifique d'approches à plus grande échelle (par ex. Stockle et al., 1994 ; Thiery, 1993 ; Douez et al., 2010). Dans tous les cas et jusqu'au début des années 2000, l'accent est mis sur le fonctionnement des écoulements souterrains (plus rarement sur le transfert de masse), les flux hydrométéorologiques, les écoulements de surface ou encore le fonctionnement de la zone vadose étant assimilés à de simples fonctions d'entrée-sortie de l'hydrosystème (par ex. Kolghi et al., 1996).

Certainement motivés par la recherche sur le climat qui, pour partie, essaye de décortiquer les grandes circulations fluides à l'échelle globale, les travaux sur le cycle continental de l'eau « redécouvrent » que les compartiments surface et sub-surface du cycle de l'eau sont intriqués (par ex. Winter et al., 1998 ; Fleckenstein et al., 2010). La physique du simple échange de volumes ou de flux d'eau entre les compartiments ne répond que très partiellement aux faits d'observation. Les couplages sont de fait très instationnaires, montrent de forts contrastes de temps caractéristiques et sont bien moins linéaires que les simples échanges mentionnés plus haut. Parallèlement, des évolutions des contraintes exercées sur les hydrosystèmes comme l'irrigation de l'agriculture intensive génèrent de nouvelles difficultés de gestion de la ressource. L'irrigation intensive pompe majoritairement des eaux de surface et des eaux de nappe d'accompagnement des rivières (Corpen, 2007). On s'aperçoit alors que de nombreux hydrosystèmes plus profonds sont rapidement en récession avec, en retour, des étiages sévères constatés dans les vallées de grands fleuves. La rotation accélérée de l'eau sur son cycle accroit également les contaminations souterraines et l'excès de minéralisation des eaux (par ex. Caswell - Zilber-man, 1986 ; Skaggs et al., 1994 ; Yadav et al., 2002). Dès lors, Les modèles hydrologiques de seconde génération apparaissent. Ils sont qualifiés de modèles « intégrés » en raison de leur capacité à traiter de l'ensemble des trois compartiments surface, zone vadose et nappe souterraine (par ex. Furman, 2008 ; Goderniaux et al., 2009 ; Frei et al., 2010 ; Weill et al., 2013). La partie hydrométéorologie de ces modèles reste encore aujourd'hui assez fruste, vraisemblablement parce que la prédiction des flux d'entrée et l'identification des paramètres météorologiques à l'échelle locale reste un domaine de recherche amont parsemé de questions encore sans réponse (bilan radiatif local, séparation évaporation-transpiration des végétaux, etc., Sophocleous, 2002 ; Mirus et al., 2011). Dans ces nouveaux modèles intégrés, les trois compartiments majeurs sont plus ou moins couplés, sachant que les approches les plus abouties traitent de la physique des écoulements dans un unique continuum surface – sol – aquifère.

Alors que les premiers modèles d'hydrosystèmes souterrains s'adressent très majoritairement aux grandes échelles et sont en partie homogénéisés, les modèles intégrés sont généralement plus locaux et de forte hétérogénéité. En premier lieu, la dynamique des flux hydriques est très contrastée entre surface, sol et nappe, grossièrement les temps caractéristiques de séjour de l'eau sont de l'ordre de quelques jours, quelques semaines à quelques mois, quelques mois à plusieurs années dans chacun des compartiments respectifs. En second lieu, de nombreuses applications ciblent la haute résolution sur des systèmes de petite et moyenne échelle (1 – 1000 km²⁾. La forte hétérogénéité de structure de nombreux réservoirs (milieux poreux fracturés, sols à porosités multiples, …) ne peut plus être occultée. L'ensemble des caractéristiques évoquées plus haut a rendu les modèles intégrés de plus en plus complexes et finalement de moins en moins bien conditionnés par les données disponibles, l'échantillonnage pratiqué étant le plus souvent insuffisant pour rendre compte de l'hétérogénéité des compartiments. Si ces données sont plus nombreuses aujourd'hui, c'est surtout parce que des variables locales continuent d'être suivies pour disposer de longues chroniques, comme : les débits des rivières, les hauteurs de nappe, la chimie des éléments majeurs. Excepté sur les cas particuliers de sites expérimentaux dévolus à la recherche, les sites d'exploitation pétrolière ou de stockage souterrain (par ex. Mari et al., 2009), on n'a guère amélioré l'analyse « en routine » de l'hétérogénéité spatiale des hydrosystèmes souterrains ou encore la mesure des flux d'échange nappe – rivière. Le bilan est d'avoir à disposition d'une part des modèles potentiellement précis, tridimensionnels et à multiples lois physiques et, d'autre part, de manquer d'information haute résolution pour documenter une application quelconque. Sans nier les avancées scientifiques et techniques des modèles intégrés, leur exploitation pratique doit passer par une phase de « réduction » qui ajuste la complexité du modèle aux capacités de conditionnement et également aux attendus souhaités.

Deux illustrations de cette réduction sont proposées dans ce travail. La première s'appuie sur le fait que les données d'observation usuelles ne distinguent pas forcément toute la géométrie des écoulements dans les hydrosystèmes. Ainsi, les mesures de débits en rivière ou de hauteurs d'eau fournissent des grandeurs moyennes qui distinguent mal les composantes tridimensionnelles de l'écoulement à surface libre d'un réseau drainant. La même constatation s'applique aux charges hydrauliques dans les nappes. Elles sont le plus souvent mesurées comme des hauteurs d'eau dans des puits complètement ouverts où le potentiel hydraulique s'équilibre à une valeur moyenne uniforme sur l'épaisseur mouillée de la nappe. Cette mesure ne discrimine pas l'hétérogénéité verticale de l'aquifère, pas plus que les composantes verticales de l'écoulement (Delay et al., 2011, 2012). Il est proposé de ramener le système drainant de surface (fossé, rivière…) à un réseau de liens monodimensionnels, chaque lien gardant néanmoins les paramètres d'une géométrie simplifiée de la section de drain modélisé. Les écoulements de sub-surface sont traités dans un continuum zone vadose – zone saturée mais en intégrant sur la verticale les propriétés hydrodynamiques du sol et de l'aquifère et les charges hydrauliques (et/ou pressions capillaires). Zone vadose et nappe sont donc modélisées comme une surface bidimensionnelle déformée, tenant ainsi compte des variations topographiques de la surface du sol et du mur (fond) de l'aquifère. La seconde illustration se focalise sur le compartiment souterrain de l'hydrosystème dont on a vu que l'hétérogénéité verticale et ses effets restaient peu visibles sur les mesures classiques de charge hydraulique. En présence de forte hétérogénéité, comme dans les milieux poreux fracturés (cas de nombreux aquifères carbonatés plus ou moins karstifiés), une approche classique associant descriptions géométriques des fractures et de la matrice poreuse se heurtera au manque d'information sur le réseau de fractures. A l'opposé, une approche simple continuum (type milieu poreux unique) devra être finement discrétisée et présenter de très forts contrastes de propriétés hydrodynamiques susceptibles de rendre compte de la chenalisation des écoulements. En gardant une approche bidimensionnelle de la nappe mais en rendant plus complexe la physique des écoulements, on peut « réduire » le modèle de nappe. L'hétérogénéité locale de l'écoulement est assimilée dans un double continuum qui décrit explicitement les flux dans le milieu ouvert du réseau de fractures et les flux dans la matrice. A l'échelle d'un site expérimental hydrogéologique, un exemple de résolution de problème inverse est traité sur des données d'interférences hydrauliques entre puits. Cet exemple montre qu'un aquifère simple milieu nécessite le choix d'une forte variabilité locale des propriétés hydrodynamiques et par conséquent une paramétrisation lourde du modèle. A l'opposé, l'accroissement de la complexité physique d'une approche double continuum décrit localement la forte hétérogénéité et n'a besoin que d'une résolution spatiale grossière de cette hétérogénéité, réduisant ainsi la paramétrisation du modèle.

II. — Réduction de dimension d’un modèle hydrologique intégré

1) Écoulement dans le réseau drainant

On limite ici le modèle hydrologique aux écoulements du réseau drainant de surface et aux écoulements de sub-surface de la zone vadose et de la zone saturée. Usuellement, les écoulements à surface libre dans le réseau drainant et les rivières sont traités par les équations de Saint-Venant (par ex. Panday & Huyakorn, 2004) que l'on va réécrire ici, après simplifications, sous la forme d'une équation d'onde diffusive monodimensionnelle (Govindaraju & Kavvas, 1990). On suppose que chaque tronçon de rivière est de section trapézoïdale définie par deux paramètres l [L] la largeur du fond du tronçon et α [-], la pente des berges (Fig. 1). On admet également que l'on peut définir un flux moyen unitaire (vitesse) u_x [LT^-1] dans une direction x normale à la section d'écoulement. La direction x n'est pas figée dans l'espace et suit la pente principale du lit de la rivière. Dans ces conditions, le système d'équations traduisant le principe de conservation des flux volumiques en eau (1a) et la conservation de la quantité de mouvement (1b) est écrit sous la forme :

Équations 1a et 1b

A [L2] est la surface en eau normale à la direction principale d'écoulement x dans la rivière, h [L] est la hauteur d'eau dans la rivière, l' [L] est la largeur de la surface d'eau libre dans la rivière, z [L] est côte topographique du fond de la rivière, q [LT^-1] caractérise l'ensemble des entrées-sorties d'eau « latérales » par les berges et le fond de la rivière, g [LT^-2] est la composante scalaire de l'accélération de pesanteur, s_f [-] est le vecteur correspondant à la pente de « friction » de la rivière.

On peut apporter une simplifi cation majeure à l'équation de conservation de la quantité de mouvement en supposant que les variations du régime d'écoulement sont lentes (écoulement pseudo-stationnaire) et rendent nul le terme d/dt =∂/∂t + u∂/∂x. Les flux latéraux q sont également considérés comme négligeables comparés au flux circulant dans le lit de la rivière. Dans ces conditions, on obtient pour l'équation (1b) :

Équation 2

Bien que la vitesse ait disparu de l'équation (2), la pente de friction réintroduit cette vitesse dans l'équation via la relation de Manning :

Équation 3

Figure 1

avec N_man [L^1/3 T^-1] le coefficient de Manning et R_h [L] le rayon hydraulique de la rivière.

Pour une section trapézoïdale de rivière, on identifie facilement la largeur en eau de la rivière l' et la section d'écoulement A en fonction de la hauteur d'eau h et des paramètres l et α (Fig. 1) :

Équation 4

Le rayon hydraulique R_h correspond au rapport entre la section A et le périmètre « mouillé » du lit la rivière :

Équation 5

L'introduction de s_f (équation 2) dans la relation de Manning (équation 3) permet d'identifier le flux moyen unitaire (vitesse) le long de la direction principale d'écoulement de la rivière :

Équation 6

En substituant U_x dans l'équation de conservation des flux (1a), on obtient :

Équation 7

Sachant que la dérivée de A par rapport au temps peut se réécrire :

Équation 8

l'équation de conservation des flux prend la forme d'une équation de diffusion non-linéaire de la hauteur d'eau h dans la rivière :

Équation 9

2) Écoulement en zone vadose et zone saturée

Le compartiment souterrain depuis la surface du sol jusqu'au fond de la nappe est assimilé à un unique continuum. Pour des questions de simplification, on occulte ici l'existence d'écoulements gravitaires en lames de surface (le ruissellement selon Horton) ou de sub-surface (écoulements hypodermiques à la faveur d'horizons tassés peu perméables du sol). Ils pourraient également être traités comme le réseau drainant (voir plus haut) mais cette fois dans une approche bidimensionnelle de l'onde de diffusion. Zone vadose et nappe sont régies par des écoulements laminaires dont la dynamique est représentée par l'équation de Darcy-Buckingham (par ex. Huyakorn & Pinder, 1983). Supposant que le milieu est peu déformable et que la masse volumique du fluide reste constante (compressibilité du fluide nulle et absence de fortes variations de concentration en solutés), on écrit :

Équation 10

ϴ [-] est la teneur en eau du milieu correspondant au ratio du volume en eau sur le volume total du milieu (pores + solide), Sw [-] est la saturation en eau correspondant au ratio du volume en eau sur le volume des pores, Ss [L^-1] est l'emmagasinement spécifique du milieu, K [LT^-1] est le tenseur de conductivité hydraulique du milieu, h [L] est la charge hydraulique (ou la hauteur capillaire), et q_w [T^-1] un terme puits-source (en volume d'eau entré ou extrait par unité de volume de milieu et unité de temps).

La réduction de la dimension d'écoulement est opérée en intégrant l'équation (10) selon la direction verticale z ou, en cas de forte pente du bassin versant, selon z pris comme la normale au mur de l'aquifère. Pour simplifier cette intégration (éviter en l'occurrence le développement d'une intégrale de Leibnitz), on suppose que les écoulements sont majoritairement parallèles au mur de l'aquifère (∇_zh=0), que les bornes d'intégration z_b, z_s, correspondant respectivement aux côtes altimétriques du mur de l'aquifère et de la surface du sol, sont de gradient local négligeable. Enfin, on sépare l'intégration selon la direction z en deux composantes en introduisant la borne d'intégration z_w correspondant à la côte de l'interface zone non-saturée/zone saturée (ou le toit d'une nappe semi-captive) :

Équation 11

La dichotomie permet de dissocier les divers comportements des paramètres hydrodynamiques entre zone vadose et zone saturée. Ainsi, les variations de teneur en eau θ sont évidemment nulles dans la zone saturée et la saturation en eau Sw vaut 1. Le terme relatif à l'emmagasinement spécifique Ss est négligeable comparé aux variations de teneur en eau dans la zone vadose, et enfin, la conductivité hydraulique K ne dépend pas de la teneur en eau dans la zone saturée. On peut donc réécrire l'équation (11) sous la forme :

Équation 12

Équation 13

Finalement l'intégration selon z permet d'établir une équation bidimensionnelle commune à la zone vadose et la zone saturée de la forme :

Équation 14

Les paramètres moyens S̅ et T̅ sont obtenus par intégration numérique, l'intégration dans la zone vadose entre z_w et z_S supposant que soient classiquement imposées des équations d'état liant conductivité K et teneur en eau θ, capacité capillaire C et charge hydraulique h. Ces mêmes équations d'état supposent également que l'on puisse établir un profil selon z de la hauteur capillaire alors que l'équation (14) n'est a priori construite que pour calculer une charge moyenne h unique le long de z. Il est choisi pour le calcul des paramètres S̅ et T̅ que le profil capillaire entre z_w et z_S soit hydrostatique, conduisant à une évolution linéaire de la hauteur capillaire entre zéro en z_w et -(z_S- z_w) en z_S. Il faut également noter que l'équation d'écoulement calcule une charge moyenne h positive dès que le milieu est localement saturé (présence d'une nappe). Une valeur négative de h (succion capillaire) ne peut être obtenue que pour un milieu complètement non saturé (absence de nappe). Usuellement, les flux latéraux (dans le plan normal à z) sont tels que : 1- dans la zone saturée le moteur des flux est le gradient des charges hydrauliques h, 2- dans la zone vadose, le moteur est le gradient des hauteurs capillaires (h négatif). Ne disposant que d'une seule valeur de h sur z, les flux latéraux totaux zone vadose et zone saturée sont en partie erronés ; par exemple, le flux dans la zone vadose peut être la conséquence du gradient d'une charge h positive correspondant à la hauteur de nappe. A l'usage, de nombreux tests sur différentes configurations d'hétérogénéité de milieu montrent que les erreurs induites par l'intégration sur z restent très faibles et non significatives à l'échelle d'un bassin versant complet (non détaillé ici).

Pour conclure sur le modèle, le couplage à l'échelle du bassin versant entre système souterrain et réseau de drainage est réalisé localement via des flux d'échange entre les deux systèmes passant par les termes puits-source de l'équation (9) et Q de l'équation (14). On considère que ces flux sont proportionnels aux différences entre hauteur d'eau dans la rivière et charge hydraulique dans la nappe (Gunduz & Aral, 2005). Cette hypothèse reste acceptable tant que la charge dans la nappe est supérieure à la côte topographique du fond de la rivière. Lorsque la nappe est plus basse, le flux d'échange devient simplement proportionnel à la hauteur d'eau dans la rivière. La réduction de dimension euclidienne du modèle intégré présenté ci-dessus permet bien évidemment une forte réduction des coûts calcul et simplifie également la linéarisation des équations pour la zone vadose et le système drainant de surface. L'objectif ici n'est pas d'entrer dans les détails de cette réduction des coûts en fonction du degré de complexité géométrique et d'hétérogénéité du système étudié. On peut néanmoins considérer sans difficulté qu'à réduire les dimensions euclidiennes d'un modèle, le temps calcul est a minima réduit au prorata du nombre d'inconnues éliminées. Par exemple, un modèle tridimensionnel à 10 « couches » de calcul (un minimum pour une vue 3D correcte) sera calculé au moins 10 fois plus vite dans son approche bidimensionnelle.

3) Exemples

Deux simulations numériques sur des objets simplifiés permettent d'illustrer le comportement du modèle intégré. Il ne sera pas procédé ici à une comparaison des résultats avec ceux d'autres approches. La géométrie simple du « livre ouvert » ou de la « vallée en Y » est paradoxalement un exemple de référence largement usité dans la littérature (par ex. Sulis et al., 2010 ; Maxwell et al., 2014). La géométrie simple du modèle ne permet pas de cacher les erreurs numériques ou de conception derrière les résultats nécessairement difficiles à appréhender sur un objet complexe. On reprendra ci-dessous la vallée en Y dont le fond est occupé par trois branches connectées d'un réseau drainant. Chaque branche est ici un fossé ou un ruisseau de 1 m de large, 1 m de profondeur et les berges sont verticales. La topographie de la vallée est telle que la pente moyenne est de 3% (50 m de dénivelé pour 1.5 km de long). Le compartiment souterrain est un aquifère d'accompagnement (en sus du sol) du réseau drainant, de faible épaisseur égale à 5 m sur l'ensemble du système. Le mur de l'aquifère suit la topographie de la surface. On teste le scénario d'un petit bassin-versant dont l'aquifère initialement à 1 m de profondeur sous la surface topographique se vidange naturellement par effet gravitaire. Les ruisseaux sont initialement secs. Au cours de cette vidange, le bassin est arrosé d'une forte pluie uniformément distribuée, d'intensité 5 mm.h^-1, qui se déclenche au temps t=48 h et dont la durée est de 4 h.

Figure 2

La Figure 2 représente les cartes de charges hydrauliques dans la nappe, représentées en profondeur de l'interface saturé/ non saturé sous la surface topographique. On observe clairement que la pluie générée à 48 h engendre en 4 h une forte augmentation piézométrique comparée à ce que serait le régime de vidange sans évènement pluvieux. Cette rapidité de réaction, difficile à simuler, est la conséquence de la faible profondeur de l'aquifère, de son caractère peu capacitif et d'un sol encore suffisamment humide (la nappe est initialement à 1 m de profondeur) pour permettre une infiltration rapide. Succédant l'évènement pluvieux, le fond du bassin garde un niveau piézométrique élevé pendant que la partie amont s'assèche rapidement (vidange de la nappe). Le caractère fortement transitoire du scénario simulé et la rapidité de réaction du système sont confirmés par l'analyse de la variation du débit à l'exutoire du Y formé par les ruisseaux. La courbe débit versus temps est globalement croissante convexe, conséquence de la vidange gravitaire de la nappe (Fig. 3). Elle devrait ensuite décroître avec l'assèchement de la nappe, mais l'évènement pluie à 48 h se manifeste par un pic de débit très asymétrique avec une croissance abrupte en raison du flux pluvieux immédiatement drainé par les ruisseaux. La relaxation nettement plus lente ramène ensuite le pic de débit vers un régime général de vidange progressive de la nappe.

Figure 3

Le modèle à dimension euclidienne réduite se comporte sans anomalie de calcul (oscillations, diffusion numérique, etc.) qui pourrait adoucir et amortir ses réponses. Ceci se vérifie tant pour un raffinement important du modèle (maille de 5 m pour un domaine de 500x1500 m, pas de temps de 4 s pour 120 h de temps simulé) que pour une approche plus grossière (maille de 20 m, pas de temps de 120 s). Les erreurs entre discrétisation fi ne et grossière sont négligeables. En pratique, l'intégration verticale des écoulements souterrains et la simplification du réseau de surface en chenaux sont des éléments qui rendent la résolution numérique plus robuste et notamment moins sensible aux hétérogénéités spatiales des paramètres ou encore leur évolution d'une itération temporelle à l'autre.

Le second exemple proposé est bâti sur une géométrie de bassin plus réaliste qui reprend celle d'un bassin-versant réel de moyenne montagne (Le Strengbach, Vosges). La géométrie est néanmoins légèrement retouchée afin d'accentuer, pour l'exercice, la réactivité du système. 72 segments de 1 m de large et 1 m de profondeur forment le réseau de drainage en surface. Le bassin-versant de 0,8 km² est à forte pente (20-30%) : de 1148 m au point haut jusqu'à 883 m au point bas sur une distance de moins d'un kilomètre de long. Cette caractéristique engendre des écoulements très majoritairement contrôlés par les effets gravitaires pour un aquifère d'accompagnement du réseau de surface et dont la faible épaisseur (8 m) est uniforme sur l'ensemble du bassin. En raison des fortes pentes locales, le maillage bidimensionnel du compartiment de sub-surface est raffiné aux abords du réseau drainant afin d'améliorer la simulation des échanges surface-souterrain. Néanmoins, la taille minimale des mailles ne descend jamais en dessous de 5 m. On simule sur une période de 120 h le drainage gravitaire rapide de l'aquifère pour une situation initiale où la hauteur de nappe initiale varie linéairement entre -8 m et -4 m sous la surface topographique entre le haut et le bas du bassin versant. On ne discute ici que de l'évolution du compartiment souterrain, évolution plus complexe que dans le premier exemple mentionné en raison d'une géométrie réaliste du bassin et de sa forte pente. Le drainage progressif de l'aquifère est décrit avec précision en espace et en temps. Les zones élevées s'assèchent pendant que les vallées se saturent, en particulier celles dont le fond est occupé par un drain de surface dont la « perméabilité » très forte facilite le transfert des flux. Globalement, pour tous les temps de simulation, les flux convergent (Fig. 4) vers le réseau drainant et la taille des zones drainées diminue avec le temps lorsque l'aquifère se vide. Les cartes piézométriques (Fig. 5) confirment l'évolution temporelle rapide du bassin et le comportement majoritairement gravitaire des écoulements en raison de la forte pente du système.

Figure 4

Figure 5

Une comparaison, non présentée ici, entre un modèle tridimensionnel et l'approche à dimension réduite ne montrerait pas de différence significative tant pour les flux du réseau drainant que les flux échangés entre compartiments et les flux souterrains. Notons néanmoins que l'exercice de simulation n'intègre pas de forte hétérogénéité du milieu souterrain dans la mesure où l'on regarde à petite échelle un aquifère peu épais accompagnant un petit réseau de surface. L'intégration d'une hétérogénéité conséquente du milieu souterrain doit être travaillée à plus grande échelle mais constitue une analyse longue au-delà des préoccupations principales de ce travail. D'ailleurs, s'il est question de traiter d'aquifères profonds à forte hétérogénéité de structure, la proposition d'une réduction de modèle via une physique plus fouillée de la dynamique des écoulements souterrains se révèle une perspective également intéressante.

III. — Réduction de modèles d'hydrosystèmes souterrains à forte hétérogénéité

Comme cela est évoqué dans l'introduction, les hydrosystèmes souterrains peuvent présenter des hétérogénéités de structure très peu « visibles » sur les données hydrauliques disponibles. C'est particulièrement le cas des aquifères carbonatés qui mixent des caractéristiques hydrodynamiques de milieux poreux et de milieux fissurés, voire karstifiés, sur des distances parfois très courtes. Le caractère diffusif des écoulements amortit très fortement les charges hydrauliques du système, passant d'une zone d'écoulement poreux à une zone d'écoulement préférentiel en drains et fractures sans observer de variation marquée des charges (par ex. Kaczmaryk & Delay, 2007a). Identifier les contrastes des propriétés de conduction et de stockage du milieu sur la base de valeurs locales de charges hydrauliques peut se révéler un exercice ardu et nécessiter un degré élevé de raffinement spatial et temporel du modèle. Il en ressort des difficultés de paramétrisation du modèle qui rendent l'exploitation de l'outil difficile tant dans sa phase de calage (inversion) que dans une phase ultérieure de prédiction et d'évaluation des incertitudes.

La plupart des modèles hydrogéologiques opérationnels a pour habitude de traiter des milieux complexes sous l'approche d'un simple continuum dans lequel on compte sur le raffinement de la discrétisation et une affectation maille par maille de propriétés hydrodynamiques contrastées pour simuler l'hétérogénéité de structure. Dans ces conditions, l'écoulement est vu au travers d'une équation unique de conservation des flux volumiques de l'équation de Darcy. En éliminant les conditions initiales et celles des limites, puis en occultant le temps et l'espace pour la variable d'état et les paramètres afin d'alléger l'écriture, le modèle est décrit par l'équation (par ex. de Marsily, 1996) :

Équation 15

avec h [L] la charge hydraulique dans la nappe, K [LT^-1] le tenseur de conductivité hydraulique, Ss [L^-1] l'emmagasinement spécifique, et q [T^-1] un terme puits-source par unité de volume de milieu. Ce sont les disparités dans l'espace (et éventuellement le temps) des valeurs des paramètres Ss et K qui décrivent l'hétérogénéité du milieu.

On peut également transcrire les contrastes de propriétés hydrodynamiques entre milieu poreux et milieu fissuré plus ouvert en supposant que le modèle est double. En chaque point de l'espace, il est susceptible de coexister un continuum fractures généralement très conducteur mais peu capacitif et un continuum matrice peu conducteur mais capacitif (par ex. Barrenblatt et al., 1960 ; Warren & Root, 1963 ; Delay et al., 2007). Les écoulements sont donc régis localement par deux équations :

Équation 16

La définition des grandeurs h, Ss, K et q est conforme à celle énoncée pour l'équation (15) sachant que, spécifiquement, les indices f et m se rapportent aux continuums respectifs fractures et matrice. a [L^-1T^-1] est un coefficient d'échange de flux entre les deux continuums. Plusieurs configurations du modèle sont possibles. Considérant que les deux milieux sont présents partout, on renseignera les équations (16) sur l'ensemble de leurs paramètres en chaque maille du domaine. On peut également supposer que seul un milieu peut exister localement, mettant alors les paramètres de l'autre milieu et l'échange a à la valeur zéro. Cette configuration reste néanmoins d'intérêt réduit puisque qu'elle est conforme à une approche simple continuum (équation 15) éventuellement bimodale sur la distribution de ses paramètres. Enfin, il est fréquemment admis que la matrice peu conductrice est de conductivité hydraulique K_m négligeable devant celle des fractures (Landereau et al., 2001). Le milieu devient stagnant, comparable à un réservoir statique qui alimente (est alimenté par) le milieu fracturé via une cinétique de premier ordre :

Équation 16

Cette cinétique a un comportement très similaire à celui d'une équation de diffusion (ici de pression) entre les deux milieux.

Les utilisateurs des modèles opérationnels de nappe rechignent généralement à s'appuyer sur une approche double milieu, prétextant : 1- que sa paramétrisation intrinsèque par maille est plus forte, 5 paramètres en excluant les termes puits au lieu de 2 pour un simple milieu, et 2- que les formes similaires des deux équations d'écoulement conduisent invariablement à des équifinalités sur les sorties du modèle sauf fort contraste entre les paramètres du milieu fracturé d'une part et les paramètres de matrice d'autre part. Il est vrai également, en cas d'inversion du modèle ou simplement d'ajustement « à la main » que les charges h_f, h_m sont de sensibilités très différentes vis-à-vis des paramètres (plusieurs ordres de grandeur). Il est, par exemple, très difficile d'ajuster le coefficient a si les autres paramètres ne sont pas déjà établis ou leur sensibilité artificiellement très atténuée (Delay et al., 2007 ; Kaczmaryk & Delay, 2007b). Pour autant, l'exemple proposé ci-dessous démontre que l'approche double continuum n'est pas sans avantage en termes de réduction de modèle, même dans l'exercice de paramétrisation optimale et d'inversion du modèle.

L'exemple proposé repose sur l'inversion d'un simple et d'un double continuum sur la base de tests d'interférences réalisés sur le Site Expérimental Hydrogéologique (SEH) de Poitiers. La description générale du site et des tests d'interférences est détaillée dans plusieurs références, notamment (Kaczmaryk & Delay, 2007b ; Audouin et al., 2008 ; Bodin et al., 2012). Le test d'interférence consiste à pomper un débit constant dans un puits et observer dans un ou plusieurs puits distants les rabattements de charge hydraulique en fonction du temps. Comme il est stipulé en introduction de cette contribution, les charges moyennes mesurées dans un puits ouvert sur l'ensemble de l'épaisseur mouillée d'un aquifère enregistrent mal les composantes verticales de l'écoulement. Par conséquent, les deux types de modélisation sont basés sur une simplification bidimensionnelle considérant des écoulements plans horizontaux. On admet également que la conductivité hydraulique des milieux est isotrope localement ce qui ramène l'identification du tenseur K à celle d'un scalaire K. L'approche double continuum considère que le milieu matrice est un réservoir statique où les écoulements sous gradient de charge sont négligeables (cf. plus haut). Spécifiquement dans le cas du SEH, l'aquifère carbonaté investigué est karstifié localement à la faveur de 4 horizons de 2-3 m d'épaisseur aux profondeurs de 30, 50, 80 et 110 m environ. Les écoulements préférentiels dans ces horizons engendrent des réponses très rapides aux tests d'interférence. Afin de produire une information spatialisée la plus exhaustive possible sur la dynamique des écoulements, on génère des jeux de données d'interférences qui, en chaque point d'observation, ajoutent séquentiellement dans le temps la réponse du puits observé à chacun des stress engendrés par chaque puits pompé. La réponse en chaque puits observé correspond donc aux rabattements de charge hydraulique générés pour n puits pompés individuellement en séquence. Chaque séquence est suivie d'une période de relaxation permettant à la nappe de revenir à sa hauteur initiale avant pompage (Fig. 6). Les jeux de données assemblés en chaque puits observé sont utilisés pour identifier par approche inverse les paramètres hydrodynamiques des deux modèles, simple et double continuums. On calcule une fonction objectif F (h(p),h*,p) correspondant à la somme du carré des résidus entre les sorties du modèle h(p) et les observations h*. Les paramètres p sont ajustés par un algorithme d'optimisation qui cherche itérativement les meilleures valeurs de p pour réduire F. Cette optimisation suppose que l'on puisse calculer les composantes dF / dp du gradient de la fonction F. Ces composantes sont obtenues par la méthode de l'état adjoint qu'un lecteur intéressé pourra retrouver dans (Ackerer & Delay, 2010 ; Ackerer et al., 2014).

Figure 6

Attendu que les deux modèles, simple et double continuums, sont programmés avec les mêmes approches numériques, le même niveau de discrétisation en espace et en temps sur le même domaine, la réduction de modèle passe ici par la paramétrisation du calcul. On cherche le nombre de paramètres spatialement distribués pour décrire l'hétérogénéité du milieu identifi able sur les données d'interférences hydrauliques. Afin de ne pas infl uer artifi ciellement sur le nombre, la distribution spatiale et la distribution statistique des paramètres à identifi er, on utilise une technique évolutive de paramétrisation (Ackerer & Delay, 2010 ; Trottier et al., 2014 ; Ackerer et al., 2014) qui positionne des points de calcul de valeurs de paramètres aux nœuds d'une grille grossière surimposée au domaine de calcul. Il n'y a pas de relation entre la grille des paramètres et la grille de résolution des équations du modèle. L'affectation des valeurs locales des paramètres sur la grille de calcul des équations s'effectue par simple interpolation linéaire des valeurs aux nœuds de la grille des paramètres. Lors des itérations de convergence minimisant la fonction objectif F, la grille des paramètres peut être affi née, augmentant ainsi le nombre de paramètres à ajuster pour caler le modèle (Fig. 7). Ce raffinement est local et différencié en fonction de la position sur la grille des paramètres. Il s'appuie sur des critères multiples comme celui de réduire localement les écarts locaux entre h(p) et h* ou encore, faire décroitre les composantes locales dF / dp . Au final, on obtient le nombre minimum et la valeur de paramètres locaux distribués sur le domaine modélisé pour ajuster les sorties du modèle aux données d'observation.

Figure 7

Reprenant l'exemple des tests d'interférence mentionnés plus haut, la Figure (8) propose deux cartes des conductivités hydrauliques K (simple continuum) et K_f (double continuum) obtenues sur la grille de calcul après inversion (Ackerer & Delay, 2010 ; Trottier et al., 2014). On rappelle qu'ayant annihilé les écoulements du milieu matrice du double continuum, la conductivité K_f du milieu fracturé devient comparable à la conductivité K d'un simple continuum. Alors que les ajustements des deux modèles sur les données d'interférence sont de qualité égale (en général, moins de 1 cm d'écart pour chaque valeur de charge entre modèle et données), les distributions spatiales de K et K_f diffèrent largement. A l'échelle de l'ensemble du domaine modélisé, les deux approches identifient sensiblement les mêmes gammes de variation des paramètres. La délinéation spatiale des zones de forte conductivité comparée aux faibles valeurs peut varier d'un modèle à l'autre. Cela étant, l'approche de paramétrisation décrite plus haut est typiquement stochastique et peut dupliquer le nombre de champs de paramètres équiprobables pour le même ajustement du modèle aux données. Regardant ces ensembles, simple et double continuum identifient les mêmes grandes zones de conductivité forte versus zones peu conductrices. La différence majeure entre approche simple et double continuum tient au degré de raffinement nécessaire pour obtenir un bon ajustement du modèle. Le simple continuum est obligé de décrire finement l'hétérogénéité locale pour rendre compte des données d'interférences. Ainsi la carte des conductivités hydrauliques K, peut localement présenter de fortes variations de valeurs locales pour rendre compte des contrastes de comportement de l'hydrosystème. En pratique, il faut un nombre important de paramètres locaux servant de point d'appui à l'interpolation pour produire le champ de conductivités. A l'opposé, le double continuum dispose intrinsèquement (dans ses équations) de la capacité à rendre compte de contrastes locaux de comportement hydraulique du milieu. La carte des conductivités ressemble à un assemblage simple de zones de valeurs uniformes plus ou moins étendues. Elle nécessite un faible nombre de paramètres (donc de points d'appui pour être établie). On rappelle que pour autant, l'ajustement des données est similaire entre les modèles. Un simple continuum à deux paramètres K et Ss nécessite entre 2x1036=2072 et 2x4108=8216 paramètres pour une inversion correcte. Le double continuum quant-à-lui est nettement moins gourmand avec pour les 4 paramètres K_f, Ss_f, ⍺ et Ss_m, un nombre de valeurs à identifier entre 4x31=124 et 4x89=356.

Quand la résolution numérique du double continuum est bien écrite, le coût calcul d'une simulation est à peine supérieur à celui affiché par un simple continuum. Bien que la procédure de raffinement de la grille des paramètres soit en partie dépendante de la « randomisation » qui permet d'obtenir plusieurs solutions équiprobables au problème inverse, on peut grossièrement considérer que le temps calcul total est proportionnel au nombre de paramètres à identifier. Dès lors et pour le cas présent, une approche double continuum s'avère entre 6 à 66 fois plus rapide à inverser qu'une approche simple continuum. La complexification apparente de la physique est à l'usage une réduction drastique du modèle qui nécessite alors bien moins d'effort d'ajustement.

Figure 8

IV. — Conclusion

Si le besoin d'une analyse plus fouillée de la dynamique du cycle de l'eau continentale s'impose dans de nombreux travaux de recherche et d'applications pratiques, ce besoin est également à l'origine d'une complexité accrue des modèles mécanistes de l'Hydrologie. Ainsi, les modèles actuels, dits « intégrés », appréhendent sous le même ensemble la complexité du fonctionnement d'un bassin versant depuis les données d'entrées hydrométéorologiques en passant par le drainage superficiel, l'infiltration dans les sols, jusqu'aux écoulements en nappe. La tridimensionnalité du système et la forte variabilité spatio-temporelle de comportement sont les éléments majeurs de complexité. Cela étant, un modèle ne reste qu'un outil générique de calcul qui doit être adapté puis renseigné sur le contexte local pour une application dédiée. De fait, le sous- échantillonnage chronique des hydrosystèmes, en particulier leur compartiment souterrain, n'aide pas à la documentation du modèle, lui-même d'autant plus demandeur d'information qu'il est complexe et finement résolu. Pire encore, une grande partie des variables observables se révèlent peu discriminantes du fonctionnement de l'hydrosystème. On a pu ainsi prendre l'exemple des flux et hauteurs d'eau en rivière ou encore des charges hydrauliques en aquifère et constater que ces variables moyennes (agrégées) restituent mal la tridimensionnalité des écoulements. Elles renseignent mal un modèle hydrologique intégré sauf à adopter la dimensionnalité du modèle à celle que la donnée est susceptible d'identifier.

Ceci nous a conduits à dégrader un modèle hydrologique sans pour autant amputer complètement la physique sous- jacente. Cette dernière est en grande partie préservée car elle ne subit qu'une phase d'agrégation permettant de réduire la dimension des écoulements du système. Il n'y a pas non plus d'incohérence conceptuelle entre réduction de dimension, variables manipulées et paramètres du modèle. Spécifiquement sur les paramètres, les valeurs moyennes utilisées ne sont que le strict reflet de l'agrégation et non pas des valeurs empiriques ou issues de connaissance a priori ou encore d'une manipulation algébrique sans relation avec la manière de réduire le modèle. Comme l'idée de réduction du modèle est motivée par la mise en conformité des variables manipulées avec la réalité des données, on compare nécessairement les sorties du modèle sans le biais de mettre face à face des objets qui n'ont pas même signification. Le fait peut s'avérer un élément incontournable lors de phases d'inversion du modèle. Si introduire de multiples mécanismes et processus, donc alourdir la description physique, peut engendrer des modèles plus complexes, paradoxalement, une physique plus affûtée peut aussi aider à la réduction. Une approche simple « entassera » sur ses propres mécanismes les effets de mécanismes absents mais partiellement capturés par des données d'observation. Il peut alors ressortir d'un exercice de comparaison modèle – données, des paramétrisations lourdes, voire irréalistes tant sur le nombre de paramètres à identifier que sur les valeurs trouvées. Prenant l'exemple d'écoulements préférentiels dans un hydrosystème carbonaté karstifié, on montre qu'une physique qui intègre localement des écoulements bimodaux à la faveur d'un continuum fracturé et d'un continuum matriciel, engendre une paramétrisation simplifiée du modèle. Les exercices d'inversion deviennent abordables, en particulier dans une approche stochastique cherchant les valeurs des paramètres, leurs distributions statistiques et spatiales, ainsi que le caractère plausible ou non d'une prédiction. S'il est indéniable qu'il faille faire progresser les modèles à la mesure des avancées dans la compréhension du fonctionnement du milieu naturel, il faut également se poser la question de l'adéquation entre les sorties du modèle et ce que peut être sa documentation, en particulier sur des données d'observation. Un idéal (mais il peut y en avoir d'autres) qui ne construirait l'expérimentation et l'observation qu'en pleine connaissance des capacités du modèle et de ses besoins n'est en pratique jamais atteint. Par conséquent, le véritable progrès dans la modélisation est typiquement un exercice à double entrée, compliquant les approches pour améliorer la description des processus physiques en jeu d'une part, mais d'autre part et en connaissance de cause, sachant réduire cette complexité en fonction des besoins et/ou des informations disponibles.