Soutenance mémoire
Le forçage climatique : apport d’eau et d’énergie
Les précipitations sont le principal apport d’eau à un versant. Sous forme de pluie, elles participent directement au cycle de l’eau. Sous forme de neige, un retard lié à la fonte est à prendre en compte. La réponse d’un bassin dépend de la nature et du volume total des précipitations apportées par un événement pluvieux. L’intensité et la durée de l’événement ont également une grande influence sur la réponse des bassins, et contrôlent par exemple l’évolution des surfaces saturées (Fig 1.4). Les précipitations ont comme principale caractéristique d’être très variables dans le temps et l’espace. Au pas de temps journalier, elles dépendent des conditions météorologiques, au pas de temps annuel, des caractéristiques climatiques de la région considérée. La variabilité spatiale est très difficile à évaluer puisqu’elle dépend principalement des variations du climat local. Cette grande variabilité spatio-temporelle peut expliquer la diversité des réponses observées sur le terrain. Les apports d’énergie sont principalement de nature radiative ou advective. Ils contrôlent la demande évaporatoire (partie de la pluie qui repart directement dans l’atmosphère par évaporation) entre deux événements et donc l’évapotranspiration réelle. L’évaporation et l’évapotranspiration influent sur la redistribution des masses d’eau dans les versants en diminuant les volumes drainés par le sol et en modifiant les réserves et les conditions hydriques de surface. Ces apports d’énergie sont également très variables dans le temps et l’espace. Une de leurs caractéristiques propres est qu’elles sont décalées dans le temps par rapport aux apports d’eau et ont donc un effet important entre deux événements pluvieux. On peut notamment citer leur rôle dans les oscillations diurnes des débits observées dans le cas de nappes affleurantes [Ambroise, 1988].
Paramétrisation du sol
Un des problèmes majeurs dans l’utilisation de cette équation concerne la paramétrisation du sol, et en particulier de la zone non saturée. Dans cette zone, la teneur en eau θ ainsi que la conductivité hydraulique K sont des grandeurs qui dépendent nonlinéairement de la pression en eau h [Hillel, 2004] (voir Fig 1.1). La plupart des modèles distribués utilise des relations fonctionnelles entre teneur en eau, conductivité hydraulique et pression capillaire. Ces fonctions caractérisent le comportement des sols à partir de paramètres, qui sont nécessaires à la bonne définition du problème. Un grand nombre de ces modèles de comportement existent dans la littérature [Mualem, 1976] [Burdine, 1953]. Les plus utilisés ont été développés par Van Genuchten [Van Genuchten, 1980] et Brooks and Corey [Brooks et Corey, 1964]. Ils associent à une valeur de pression capillaire une valeur unique de teneur en eau ou de conductivité hydraulique. Néanmoins, ils ne sont qu’une approximation grossière du comportement complexe des sols. L’évolution de la conductivité hydraulique ou de la teneur en eau est très souvent différente si le sol est en phase de drainage ou d’imbibition. Ce phénomène, appelé hystérésis, est rarement pris en compte dans les modèles distribués. La définition de ces fonctions est, de plus, difficile car la mesure des courbes caractéristiques de sol représente un investissement expérimental non négligeable. Les études de laboratoires ont montré que ces courbes dépendent beaucoup de la texture des sols. La Fig 1.1 schématise de façon simple les lois de teneur et de conductivité hydraulique pour trois types de sol (une argile, un sable marneux et un sable). La grandeur notée Lc sur ces courbes est appelée longueur capillaire et caractérise la hauteur de la frange capillaire (zone de tension au dessus du niveau de la nappe). On peut voir sur cette illustration qu’un sable, qui a une structure grossière, a une conductivité hydraulique à saturation forte, une longueur capillaire faible et que l’évolution de la conductivité hydraulique et de la teneur en eau sont beaucoup plus non linéraires que pour les deux autres types de sol. En revanche, les sols à texture plus fine ont une conductivité hydraulique à saturation plus faible, une longueur capillaire plus grande et présentent des évolutions de teneur en eau et de conductivité hydraulique plus douces. Ceci est lié à la taille des particules solides et leur organisation qui déterminent la distribution des pores et donc le comportement du sol. Il faut noter qu’en terme de modélisation ces lois jouent un rôle primordial et ont une grande influence sur le comportement simulé des systèmes (e.g. [Vogel et al, 2001], [Van Genuchten et Nielsen, 1985]).
Historique des méthodes de couplage
Les premières approches pour prendre en compte les interactions surface/subsurface ont été implémentées dans des modèles ne prenant en compte que les processus de surface [Wooding, 1965] [Chen et Chow, 1968]. Les écoulements de subsurface sont pris en compte par l’intermédiaire de termes puits estimés à partir de théories simplifiées de l’infiltration. Ce genre de modèles est encore développé aujourd’hui. Le modèle de Esteves [Esteves et al, 2000], par exemple, résout les équations de Saint Venant complètes en évaluant le terme puits d’infiltration à partir de l’équation de Green et Ampt. Le but de ce modèle n’est pas de modéliser le cycle de l’eau à l’échelle d’un bassin versant mais d’étudier le ruissellement hortonien à l’échelle d’une parcelle. Ce modèle a également servi de base au développement d’un modèle distribué d’érosion [Nord et Esteves, 2004]. L’approche la plus courante pour traiter le couplage surface/subsurface s’inspire du formalisme de Freeze et Harlan [Freeze et Harlan, 1969] et s’appuie sur la méthode dite de changement de condition à la limite. Elle est notamment utilisée dans les travaux suivants [Smith et Woolhiser, 1971], [Freeze, 1972a, Freeze, 1972b], [Beven et Kirkby, 1979], [Smith et Herbert, 1983],[Abbott et al, 1986a, Abbott et al, 1986b], [Govindaraju et Kavvas, 1991], [Querner, 1997], [Branstert et Plate, 1997]. Dans ces travaux, la condition à la limite à la surface du sol est une condition de flux imposé, égal à la pluie, jusqu’à ce que l’accumulation d’eau en surface se produise. A ce moment de la simulation, cette condition à la limite est changée et une pression égale à la hauteur de lame d’eau est imposée à la surface du sol. Cette méthode a déjà été décrite et traduite en équations précédemment (voir Eq 1.13 et Fig 2.1) dans la partie traitant des écoulements de subsurface. Cette approche présente des difficultés tant numériques qu’algorithmiques [Vanderkwaak, 1999]. La réponse de la subsurface est en effet transitoire et très variable spatialement en fonction des propriétés du milieu. Déterminer la zone sur laquelle imposer une pression et la position du point de débordement représente donc un véritable défi numérique [Beven, 1985]. Des schémas de type prédicteur/correcteur sont souvent utilisés pour déterminer avec précision la position de ce point et les pas de temps employés sont souvent petits. De plus, avec cette méthode, l’évaluation de la partition infiltration/ruissellement au niveau de la surface est relativement complexe. Malgré ces difficultés, de nombreux travaux s’appuyant sur cette méthode ont été réalisés et vont maintenant être présentés plus en détails. Smith et Woolhiser [Smith et Woolhiser, 1971] sont les premiers à avoir tenté de coupler surface et subsurface sur la base des travaux de Freeze et Harlan. Dans leur modèle, le ruissellement est décrit par une équation d’onde cinématique et l’infiltration par une équation de Richards 1D verticale. L’algorithme de couplage s’appuie sur la méthode de changement de condition à la limite. Une fois que la surface est saturée et que le processus de ruissellement apparaît, le flux s’infiltrant i(t) est calculé en tout point à partir de l’équation de Richards avec une condition de pression imposée en surface. Le flux alimentant le processus de ruissellement q(t) est déduit de la pluie r(t) et du flux infiltré i(t) de façon simple puisque : q(t) = r(t) − i(t) (2.1) L’équation d’onde cinématique est alors résolue avec ce terme source. Un algorithme de choix de pas de temps a été développé pour déterminer avec précision le moment auquel il faut changer de condition à la limite en surface. Ce point représente une des difficultés majeures évoquées précédemment et l’algorithme mis en place dans cet article semble très coûteux. Ce modèle ne permet de décrire que le ruissellement de type hortonien et a été testé sur une expérience de laboratoire et sur un petit bassin versant. Freeze [Freeze, 1972a, Freeze, 1972b] étudie le rôle de la subsurface dans la genèse des écoulements dans les rivières à partir d’un modèle différences finies couplant une équation 3D de subsurface et une équation 1D décrivant les écoulements en rivières. Les deux équations sont reliées de façon externe : la solution de l’équation de subsurface donne les conditions aux limites pour l’équation d’écoulement dans la rivière et inversement. Le phénomène de ruissellement n’est pas modélisé puisque seul le domaine de subsurface est considéré. L’eau s’exfiltrant par les surfaces de suintement est directement routée vers la rivière en introduisant un retard lié au temps de transfert à la surface du sol [Freeze, 1971]. Ces simulations démontrent le rôle important de la subsurface dans la genèse des débits, l’influence de paramètres tels que l’intensité de la pluie, sa localisation ou les propriétés du sol et remettent en question le modèle conceptuel du ruissellement hortonien en montrant la rareté de son occurrence. Beven [Beven, 1977] a complété le travail de Freeze en implémentant un modèle similaire avec la méthode des éléments finis. Il a montré que la réponse d’un versant était fortement non-linéaire en fonction de la pluie et de différents paramètres et que la connaissance des conditions initiales, particulièrement dans la zone non saturée, était très importante pour déterminer la réponse d’un versant. Le Système Hydrologique Européen (SHE) [Beven, 1985] [Abbott et al, 1986a] [Abbott et al, 1986b] [Bathurst, 1986a] [Bathurst, 1986b] est considéré aujourd’hui comme le modèle hydrologique le plus accompli. SHE permet de décrire les principales composantes du cycle de l’eau à l’échelle d’un bassin versant : les écoulements en rivière (équation 1D), le ruissellement (équation d’onde diffusive 2D), les écoulements en zone non saturée (équation de Richards 1D) et les écoulements dans la zone saturée (équation de diffusivité 2D). Les phénomènes d’interception, d’évaporation et les effets liés à la prise en compte de la neige et de sa fonte peuvent être ajoutés. Le couplage entre les différents processus est réalisé à l’aide de la méthode de changement de condition à la limite. Ce modèle permet de simuler l’ensemble des processus de genèse de débits décrits précédemment (ruissellement hortonien, ruissellement sur surface saturée, exfiltration, …). Bien que ce modèle ait été testé et validé de façon intensive à différentes échelles, il semble que la convergence, notamment au niveau des conditions aux limites, soit difficile à atteindre [Beven, 1985]. Ce modèle est encore développé aujourd’hui pour introduire des processus de transport, d’érosion ou prendre en compte la géochimie. Les travaux de Govindaraju [Govindaraju et Kavvas, 1991] ont permis d’avancer considérablement dans la modélisation du couplage surface/subsurface. Dans son modèle, Govindaraju couple à l’aide du changement de condition à la limite des écoulements 2D en zone saturée/non saturée avec des écoulements de surface 1D (ruissellement et rivière). La nouveauté de son travail est qu’il introduit un modèle de transport à la fois dans son domaine de surface et de subsurface. Son étude de l’influence de la géométrie, du type de sol et d’autres paramètres est particulièrement complète et fait figure de référence dans le domaine de la modélisation des processus de surface et de subsurface. Il ressort de la plupart de ces travaux que la méthode de couplage par changement de condition à la limite est relativement compliquée à mettre en oeuvre. Les algorithmes nécessaires à la convergence et à la stabilité sont complexes comme peuvent l’illustrer les travaux de Beaugendre et al [Beaugendre et al, 2006]. De nouvelles méthodes de couplage ont donc été imaginées ces dernières années pour prendre en compte de façon intégrée l’ensemble des processus de surface et de subsurface.
Comparaison avec la solution analytique de l’approximation de l’onde cinématique
Deux approximations des équations de Saint Venant complètes sont couramment utilisées pour modéliser les écoulements de surface : les approximations de l’onde cinématique et de l’onde diffusive (cf partie 2, paragraphe 1.2.2). Le modèle implémenté dans Cast3M utilise l’approximation de l’onde diffusive. Malheureusement, cette équation ne présente pas de solution analytique. En revanche, l’équation de l’onde cinématique a une solution analytique dans un cas simple monodimensionnel [Eagleson, 1970]. Pour évaluer notre modèle de ruissellement, les résultats de notre modèle sont comparés à cette solution analytique. On simule la réponse d’une surface imperméable monodimensionelle (type surface goudronnée) à un événement pluvieux de durée finie. Le domaine considéré a une longueur de 183 mètres, une pente de 0.0016 et un coefficient de Manning de 0.025 sm−1/3. Une pluie d’intensité 1.4 × 105ms−1 et d’une durée de 30 minutes est imposée à la surface du domaine. L’influence des discrétisations spatiale et temporelle est étudiée. La figure 2.1 compare la solution analytique de l’onde cinématique et deux hydrogrammes simulés avec deux pas de temps différents (1 seconde et 10 secondes) et un nombre de mailles fixe (10 mailles). La figure 2.2 présente cette même solution analytique et deux hydrogrammes simulés avec le même pas de temps (1 seconde) mais deux discrétisations spatiales différentes (10 et 100 mailles). Ces figures montrent que l’accord entre les hydrogrammes simulés et la solution analytique de l’onde cinématique est bon. Les parties montantes et descendantes de l’hydrogramme sont très bien représentées par notre modèle de ruissellement. En revanche, le pic de ruissellement n’est pas très bien décrit. Comme prévu, notre approche de modélisation approche de façon plus douce le pic de ruissellement que l’approximation de l’onde cinématique. Ceci est lié au fait que l’on introduit un terme diffusif à l’équation de l’onde cinématique pour obtenir l’équation de l’onde diffusive. La Fig 2.2 montre que, pour une simulation de ruissellement pur, l’impact de la discrétisation spatiale est faible. En effet, les deux hydrogrammes obtenus avec deux discrétisations spatiales différentes sont très proches. En revanche, la Fig 2.1 montre qu’un choix judicieux de pas de temps est nécessaire à la bonne représentation du phénomène de ruissellement. En effet, quand on augmente le pas de temps, la précision de notre approche de modélisation diminue. Ceci est lié à la condition à la limite itérative employée en sortie de domaine. Comme prévu, on constate que le flux en sortie est sous-estimé dans la phase ascendante de l’hydrogramme et qu’il est surestimé dans la phase descendante. La condition à la limite itérative permet de suivre l’évolution de la lame d’eau en sortie mais elle introduit un retard qui peut expliquer les différences observées. Ces effets liés à la condition à la limite utilisée sont moins marqués lorsque le pas de temps diminue. Pour bien suivre l’évolution de la hauteur de lame d’eau à l’exutoire du système, un petit pas de temps est donc nécessaire. Si le pas de temps employé est trop grand, des gradients de hauteurs de lame d’eau peuvent se créer et affecter la qualité des hydrogrammes simulés.
|
Table des matières
Introduction générale
Partie I Les interactions surface/subsurface au coeur de la genèse des débits
Chapitre 1 Les processus de genèse des débits
1.1 Les processus souterrains
1.1.1 Écoulement hypodermique
1.1.2 L’effet piston
1.1.3 Intumescence de nappe
1.1.4 Écoulement en macropores
1.2 Les processus de surface
1.2.1 Ruissellement par dépassement d’infiltrabilité
1.2.2 Ruissellement sur surface saturée
Chapitre 2 Les facteurs de contrôle
2.1 Le forçage climatique : apport d’eau et d’énergie
2.2 Les conditions initiales
2.3 Les propriétés du milieu
2.4 La topographie
Partie II La modélisation distribuée des processus de surface et de subsurface
Chapitre 1 Modélisation des processus
1.1 Écoulements de subsurface
1.1.1 Écoulements dans la zone saturée
1.1.2 Écoulement en zone non saturée
1.1.3 Paramétrisation du sol
1.1.4 Conditions aux limites
1.1.5 Résolution numérique
1.2 Écoulements de surface
1.2.1 Les équations de Barré SaintVenant
1.2.2 Approximations dites de l’onde diffusive et de l’onde cinématique
1.2.3 Les conditions aux limites
1.2.4 Résolution numérique
Chapitre 2 Les différentes méthodes de couplage
2.1 Historique des méthodes de couplage
2.2 Le couplage du premier ordre
2.3 Le « nouveau » changement de condition à la limite
Partie III Le modèle développé
Introduction
Chapitre 1 Le modèle physique
1.1 Description unifiée des écoulements de subsurface et souterrains
1.2 Les écoulements de surface
1.3 L’approche darcéenne multidomaine
1.4 Le modèle de transport
Chapitre 2 L’environnement de modélisation
2.1 Cast3m
2.2 Les éléments finis mixtes hybrides
Chapitre 3 Le modèle numérique
3.1 La résolution numérique
3.2 Introduction d’une conductivité hydraulique résiduelle
3.3 Continuité surface/subsurface
3.4 Condition à la limite de charge imposée “itérative”
Partie IV Validation/Vérification
Introduction
Chapitre 1 Vérification du modèle d’écoulement/transport en milieu poreux
1.1 Description du cas test
1.2 Résultats
Chapitre 2 Validation du modèle de ruissellement
2.1 Comparaison avec la solution analytique de l’approximation de l’onde cinématique
2.2 Ruissellement bidimensionnel sur une surface de type « livre ouvert »’
2.2.1 Description du cas test
2.2.2 Résultats
Chapitre 3 Vérification de l’approche de couplage par tenseur
3.1 Régime d’infiltration pure
3.1.1 Description
3.1.2 Résultats
3.2 Régime de Horton
3.2.1 Description
3.2.2 Résultats
Chapitre 4 Vérification du modèle couplé
4.1 Le système d’Abdul et Gillham
4.2 Le système d’Ogden et Watts
Partie V Applications
Chapitre 1 Sensibilité du système d’Abdul et Gillham à différents paramètres physiques
1.1 Influence du type de sol
1.2 Influence de l’intensité de la pluie
1.3 Influence de la topographie
1.4 Influence de la condition initiale
Chapitre 2 Le régime de Horton
2.1 Influence de l’intensité de la pluie
2.2 Influence de la condition initiale
Chapitre 3 Application à la parcelle de Thies – Sénégal
3.1 La parcelle de Thies
3.2 L’expérience du 28 juin 2004
3.3 Travaux antérieurs
3.4 Application de notre modèle à la parcelle de Thies
3.4.1 Démarche de modélisation
3.4.2 Résultats
3.5 L’expérience du 24 juin
3.5.1 Démarche de modélisation
3.5.2 Résultats
Performances du modèle, convergence, stabilité
Conclusion générale et perspectives
Bibliographie
Télécharger le rapport complet
