These
These
These
par
Pierre SOCHALA
A mes parents,
Michel et Martine.
Remerciements
Je suis très reconnaissant à Peter Bastian et à Stéphane Cordier du temps qu’ils m’ont
consacré en tant que rapporteur sur cette thèse. Mes remerciements s’adressent également
à Robert Eymard, Cyril Kao et Frédéric Nataf pour leur présence dans le jury, témoignant
de l’intérêt qu’ils portent à ce travail.
Je remercie Karim Djadel et Daniele di Pietro avec lesquels j’ai travaillé lors de l’implémen-
tation des méthodes de Galerkine discontinues ainsi que Thomas Esclaffer et Bénédicte
Augeard qui m’ont permis de mieux comprendre l’hydrologie. Je salue aussi mes collègues
du CERMICS et de l’ENPC en particulier Jean–Philippe Chancelier et Roland Jarry
pour leur aide lors de mes soucis ponctuels en informatique. Pablo Tassi fut un excellent
camarade de bureau : sa bonne humeur et son café m’ont accompagnés quotidiennement
au cours de la troisième année.
Jean–Pierre Chaquin m’a toujours orienté avec attention lors des moments importants de
mon parcours. Mon intérêt pour la simulation numérique s’est révélé pendant son cours
de mathématiques générales à l’ESTP. Je lui suis sincèrement reconnaissant de son sens
pédagogique, de sa disponibilité et de sa gentillesse.
Cette thèse n’aurait pu être réalisée sans le soutien infaillible de mes parents qui m’ont
toujours encouragé dans mes choix et donné les moyens de faire ce que j’aime. Magalie m’a
considérablement épaulé par ses conseils, son écoute et sa patience de chaque jour. Je les
remercie profondément, ainsi que ma famille et mes amis, pour tout ce qu’ils m’apportent.
.
Table des matières
Résumé 1
Chapitre 1 Introduction 3
1.1 Contexte et motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Le ruissellement : processus et enjeux . . . . . . . . . . . . . . . . . . . . 3
1.1.2 Quels modèles en hydrologie ? . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Modélisation des écoulements en milieu poreux . . . . . . . . . . . . . . . . . . . 5
1.2.1 Milieu poreux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2.2 Ecoulements diphasiques eau-air . . . . . . . . . . . . . . . . . . . . . . . 7
1.2.3 Propriétés hydrodynamiques des sols . . . . . . . . . . . . . . . . . . . . . 9
1.3 Modélisation des écoulements à surface libre . . . . . . . . . . . . . . . . . . . . . 12
1.3.1 Equations de Saint–Venant . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3.2 Approximation de l’onde cinématique . . . . . . . . . . . . . . . . . . . . 12
1.4 Objectifs de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.5 Plan de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
Conclusion 95
Bibliographie 116
Méthodes numériques pour les écoulements souterrains
et couplage avec le ruissellement
Résumé. Des schémas numériques précis et robustes sont proposés pour modéliser les écoule-
ments souterrains et leur couplage avec le ruissellement surfacique. Les écoulements souterrains
sont décrits par l’équation de Richards (instationnaire) qui est discrétisée par une méthode BDF
en temps et une méthode de Galerkine discontinue à pénalisation intérieure symétrique en espace.
Des cas tests sur des colonnes d’infiltration confirment la robustesse des schémas choisis. Dans un
premier temps, nous considérons des conditions de Signorini pour l’équation de Richards afin de
modéliser la présence de drains en fond d’aquifère ou l’affleurement de la nappe en négligeant
le ruissellement, c’est-à-dire en supposant que l’eau exfiltrée est immédiatement évacuée du
système. Dans un second temps, nous prenons en compte le ruissellement par le biais de condi-
tions de couplage qui imposent l’égalité des flux d’eau échangés et la continuité de la pression
à l’interface. Les écoulements superficiels sont décrits par l’équation de l’onde cinématique qui
constitue une approximation des équations de Saint–Venant. L’équation de l’onde cinématique
est discrétisée par une méthode de Godunov. Les deux schémas, pour l’écoulement souterrain et
pour l’écoulement superficiel, sont conservatifs et peuvent être utilisés dans des algorithmes de
couplage faisant intervenir un ou plusieurs pas de temps. Pour assurer la conservation de la masse
d’eau totale du système couplé, les flux à l’interface doivent être convenablement choisis. Nous
donnons en particulier la construction de ces flux pour les schémas BDF1 et BDF2. La précision
et la robustesse de nos schémas sont évaluées sur plusieurs cas tests dont le drainage d’une lame
d’eau, deux cas d’exfiltration de nappe (l’un provoqué par la pluie et l’autre par une injection en
fond d’aquifère) et un ruissellement hortonien. Enfin, nous présentons une application concrète
portant sur le fonctionnement hydrologique d’un petit bassin versant drainé.
Mots-clefs : Equation de Richards, approximation de l’onde cinématique, éléments finis dis-
continus, BDF, schéma de Godunov, algorithmes de couplage, conservation de la masse.
Numerical methods for subsurface flows and coupling with surface runoff
Abstract. Accurate and robust numerical schemes are proposed to simulate subsurface flows
and their coupling with surface runoff. Subsurface flows are modelled by the (unsteady) Richards’
equation discretized by a BDF in time and a discontinuous Galerkin method with symmetric
interior penalty in space. Infiltration column test cases confirm the robustness of our schemes.
Firstly, we consider Signorini conditions for the Richards’ equation to model buried drains or the
water table reaching the ground surface ; thus we assume that exfiltrated water immediately exits
the system. Secondly, runoff flow is taken into account through coupling conditions enforcing
water flux equality and pressure continuity at the interface. Overland flows are modelled by
the kinematic wave approximation of the shallow water equations, which is discretized by a
Godunov method. Both schemes, that for the subsurface flow and that for the overland flow,
are mass conservative and can be coupled within a multiple time step procedure. To ensure
total water mass conservation for the whole system, interface fluxes must be carefully designed.
We specify the form of these fluxes for BDF1 and BDF2. Accuracy and robustness of our
schemes are assessed on drainage, exfiltration and hortonian runoff test cases. Finally, we apply
the methology to investigate the hydrological behavior of a small-scale drained watershed.
keywords : Richards’ equation, kinematic wave approximation, discontinuous finite elements,
coupling algorithms, mass conservation.
Chapitre 1
Introduction
Cette thèse porte sur la simulation numérique des écoulements souterrains en milieu poreux
ainsi que leur couplage avec les écoulements surfaciques. L’infiltration de l’eau dans le sol, son
transport dans la zone variablement saturée du sol et dans la nappe phréatique ainsi que son
ruissellement sur la surface du sol sont des processus complexes régis par des lois physiques
différentes. Ces phénomènes ont chacun un comportement propre et la compréhension de leurs
interactions est cruciale en hydrologie.
Dans ce chapitre introductif, nous commençons par indiquer le contexte et les motivations
de notre travail, puis nous présentons les différents modèles utilisés pour décrire les écoulements
en milieu poreux, en précisant leurs domaines de validité. Nous donnons ensuite les équations
régissant les écoulements à la surface du sol en indiquant les simplifications qui sont possibles.
Nous précisons également les objectifs de cette thèse en écrivant le système d’équations couplées
que nous proposons de résoudre numériquement. Nous concluons ce chapitre en détaillant le
plan de ce mémoire.
Suivant les conditions d’humidité et les propriétés des sols en aval, ce ruissellement peut se ré-
infiltrer et contribuer à une recharge de la nappe ou rejoindre une rivière. Une connaissance de
l’évolution du ruissellement nécessite donc de connaı̂tre la capacité d’infiltration des zones tra-
versées et plus généralement les interactions entre les écoulements de surface et les écoulements
souterrains.
Pendant les périodes de fortes pluies, le ruissellement alimente de façon plus importante les
cours d’eau et participe à la formation de crues. Le débat sur la provenance de l’eau des crues,
c’est-à-dire sur la répartition de l’eau de pluie ruisselée et de l’eau provenant de la nappe, a
toujours lieu comme le montre la synthèse proposée par Kienzler et Naef [KN08]. Les hydrologues
cherchent à mieux comprendre et anticiper l’importance du ruissellement qui peut engendrer des
dégâts dans les zones habitables ou cultivées et sur le milieu naturel. En effet, lors d’inondations,
les dommages sur les biens matériels sont souvent considérables et concernent parfois direc-
tement les populations. La pollution de l’eau est un autre risque accentué par le ruissellement.
En milieu urbain, lorsque les systèmes d’assainissement des eaux usées sont saturés par des
pluies d’intensité importante, le ruissellement des eaux pluviales non traitées détériore la qualité
de l’eau utilisée dans les villes et a un impact négatif sur le milieu naturel. De même, en milieu
rural, une partie des engrais et des pesticides répandus dans les champs peut se retrouver dans
les eaux de ruissellement et être nocive pour les écosystèmes en diminuant la qualité de l’eau.
Enfin, l’érosion des sols constitue un troisième risque rencontré en milieu agricole et amplifié
par le ruissellement. La conséquence principale de ce processus naturel est la diminution de
la couche superficielle du sol, ce qui peut engendrer la destruction des semis, la baisse de la
fertilité du sol et des coulées de boue en aval. L’étude du ruissellement nécessite par ailleurs
de prendre en compte les aménagements humains réalisés en milieux urbains et agricoles qui
modifient le chemin de l’eau (fossé, présence de haies, tuyau de drainage, imperméabilisation,
bassin de retenue, bandes enherbées ...).
L’intérêt de développer des modèles numériques en hydrologie est triple. Les modèles per-
mettent d’abord une meilleure compréhension des facteurs contrôlant les processus physiques en
s’affranchissant de certaines limites des mesures expérimentales notamment pour faire varier les
différents paramètres comme la géométrie, l’intensité de la pluie et les conditions initiales et aux
limites. Par ailleurs, ces modèles constituent dans certains cas des outils utilisables pour prévoir
1.2. Modélisation des écoulements en milieu poreux 5
le comportement de systèmes soumis à des sollicitations extrêmes qui favorisent l’apparition des
risques liés au ruissellement. Enfin, ils peuvent également permettre de tester l’effet de nouveaux
aménagements sur l’hydrologie. Ce développement de modèles numériques pour l’hydrologie est
d’actualité comme le montre le récent Projet ANR « Méthode » [Mét] dont le but est d’étudier
les effets des hétérogénéités de surface sur les écoulements. Plus précisément, il s’agit d’étudier
l’impact des sillons agricoles sur le ruissellement des eaux de pluie.
un fluide liquide ou gazeux. Ces deux éléments renvoient à deux caractéristiques essentielles dans
la description d’un milieu poreux :
– la porosité ϕ [m3 /m3 ] qui est définie comme le rapport entre le volume des pores et le vo-
lume total dans un volume élémentaire représentatif du milieu. Ce paramètre géométrique
est compris entre 0 et 1.
– la perméabilité intrinsèque ki [m2 ] qui indique l’aptitude du milieu à être traversé par
un écoulement. Cette grandeur ne dépend que de la structure et de la connectivité des
pores. Lorsque le milieu est isotrope, c’est-à-dire que ses propriétés sont indépendantes de
la direction, la perméabilité intrinsèque est un tenseur proportionnel au tenseur identité.
Lorsque le milieu est anisotrope, la perméabilité intrinsèque est un tenseur symétrique.
Plusieurs approximations analytiques de la perméabilité intrinsèque existent en fonction
du type de milieu poreux. Nous citons, entre autres, le modèle de Hazen [Haz11] pour les
milieux granulaires, celui de Saffmann [Saf59] pour les réseaux de canaux parallèles et celui
de Kozeny–Carmann [Koz27] pour les réseaux de canaux tortueux. Des calculs récents sur
des réseaux de sphères ont aussi été effectués par Tardif, Ern et Dormieux [Td06, TdED07].
+ + =
Deux autres éléments relatifs à chaque fluide présent dans les pores sont importants pour
comprendre les milieux poreux :
– la saturation s [m3 /m3 ] d’un fluide qui est le rapport entre le volume de ce fluide et
le volume des pores dans un volume élémentaire représentatif du milieu. Ce paramètre
représentant la fraction du fluide dans les pores est compris entre 0 et 1.
– la perméabilité relative kr [−] d’un fluide donné qui varie en fonction de la saturation dans
le milieu et qui traduit l’aptitude de ce fluide à s’écouler. La perméabilité relative est une
quantité scalaire. Plusieurs auteurs ont proposé des relations empiriques pour déterminer
la perméabilité relative pour l’eau et l’air. Nous citons entre autres Corey [Cor54], Fatt et
Klikoff [FK59] ainsi que Bazant et Thonguthai [BT78].
Les caractéristiques intrinsèques des fluides s’écoulant dans les pores sont la viscosité dynamique
µ [P a.s] et la masse volumique ρ [kg.m−3 ]. Les grandeurs inconnues sont la saturation s et la
pression p [P a]. Dans les milieux poreux, l’écoulement est multiphasique quand il fait intervenir
plusieurs phases et monophasique quand il fait intervenir une seule phase. Par ailleurs, le milieu
est dit saturé lorsque les pores sont complètement occupés par l’eau en phase liquide et insaturé
dans le cas contraire.
1.2. Modélisation des écoulements en milieu poreux 7
Equations diphasiques
Les variables concernant l’air portent l’indice a et celles concernant l’eau portent l’indice w.
Dans notre cas, nous supposons que l’espace libre du sol est occupé par de l’air et de l’eau si
bien que la saturation en eau sw et la saturation en air sa vérifient
sw + sa = 1.
Les écoulements diphasiques sont gouvernés par les équations de continuité des deux phases qui
traduisent la conservation de la masse de l’eau et de l’air
Les vitesses de chaque phase sont données par la loi de Darcy. Cette loi expérimentale, fonda-
mentale en hydraulique des sols, exprime une relation affine entre la vitesse de chaque phase
et le gradient de pression pj de cette phase. Elle fut proposée par Darcy [Dar56] et la version
tridimensionnelle est
krj (x, pj )
vj (x, t) = −ki(x) ∇pj + ρj ~g pour j ∈ {a, w}.
µj
Equation de Richards
Le système d’équations précédent peut être simplifié sous les trois hypothèses suivantes :
1. le squelette est indéformable,
2. l’eau est incompressible,
3. l’air circule librement dans le milieu poreux.
8 Chapitre 1. Introduction
La première hypothèse implique que la porosité est constante. La deuxième hypothèse signifie que
la masse volumique de l’eau est constante. La troisième hypothèse se traduit par une viscosité
de l’air nulle ce qui permet de déduire d’après la loi de Darcy que la pression de l’air est
hydrostatique :
µa = 0 ⇒ ∇(pa + ρa gz) = 0, (1.1)
où z désigne la coordonnée verticale orientée vers le haut. Comme la quantité ρa gz est négligeable
devant la pression atmosphérique (pour z de l’ordre de quelques dizaines mètres), l’équation (1.1)
implique en général que la pression de l’air est constante et égale à la pression atmosphérique.
La masse volumique de l’eau et la porosité étant constantes, l’équation de continuité de l’eau se
simplifie
ϕ∂t sw + ∇ · vw = 0. (1.2)
La loi de Darcy portant sur l’eau peut se réécrire, si l’on tient compte du fait que le terme de
gravité dérive du potentiel gravitaire, de la façon suivante
krw (pw )
vw = −ki ∇(pw + ρw gz). (1.3)
µw
Pour alléger les notations, nous omettons la dépendance en la variable d’espace.
En hydrologie, il est courant d’introduire les trois grandeurs suivantes :
def
– la teneur en eau θ [m3 /m3 ] regroupant la porosité et la saturation en eau : θ = ϕsw ,
def
– la conductivité hydraulique K [m.s−1 ] définie par K = ki.krw (pw )ρw g/µw ,
def
– la charge hydraulique ψ [m] qui est reliée à la pression par ψ = pw /(ρw g).
La teneur en eau étant reliée à la charge hydraulique par la courbe de pression capillaire, nous
écrivons θ(ψ) pour indiquer cette dépendance (voir la section 1.2.3 pour plus de détails). Avec
ces notations et en désignant par v(ψ) la vitesse de l’eau, les équations (1.2) et (1.3) s’écrivent
∂t [θ(ψ)] + ∇ · v(ψ) = 0,
v(ψ) = −K(ψ)∇(ψ + z).
L’élimination de la vitesse aboutit à l’équation de Richards [Ric31] où l’inconnue est la charge
hydraulique
∂t [θ(ψ)] − ∇ · (K(ψ)∇(ψ + z)) = 0. (1.4)
Si l’on suppose que la teneur en eau est une fonction dérivable de la charge hydraulique, il est
possible de réécrire la dérivée temporelle de l’équation précédente en introduisant la capacité
capillaire C [m−1 ] qui est la dérivée de la teneur en eau par rapport à la charge. L’équation de
Richards s’écrit alors
C(ψ)∂t ψ − ∇ · (K(ψ)∇(ψ + z)) = 0. (1.5)
L’équation de Richards peut également avoir comme inconnue la teneur en eau. Cette dernière
formulation fait intervenir la diffusivité hydraulique D [m2 .s−1 ] qui est le rapport entre la conduc-
tivité hydraulique et la capacité capillaire
la conductivité hydraulique étant alors une fonction de la teneur en eau. Ginting [Gin04] présente
les avantages et les inconvénients des formulations (1.5) et (1.6). La formulation en θ correspond
1.2. Modélisation des écoulements en milieu poreux 9
à une loi de conservation avec terme source mais cette formulation n’est plus applicable si le
milieu devient saturé puisque la teneur en eau est constante et la diffusivité tend vers l’infini.
Par ailleurs, la teneur en eau n’étant pas nécessairement continue entre deux sols de nature
différente, cette forme n’est valable que pour des sols homogènes. Par contre, la charge hydrau-
lique étant continue dans le sol, la formulation (1.5) peut s’employer dans des sols hétérogènes.
Cette version s’étend également aux milieux contenant des parties insaturées et saturées, mais
elle n’est pas conservative. La forme (1.4) présente les trois avantages précédents puisqu’elle est
conservative et s’applique dans des milieux hétérogènes pouvant être saturés.
Modèle de Green–Ampt
Un modèle proposé par Green et Ampt [GA11] est parfois utilisé en hydrologie, comme
dans les travaux de Govindaraju et al. [GKJR96] et de Esteves et al. [EFGV00], lorsque le sol
est supposé homogène, de perméabilité constante et de texture grossière. Les deux hypothèses
fondamentales de ce modèle sont d’une part que les transferts horizontaux sont négligeables si
bien que le milieu est considéré comme une juxtaposition de colonnes d’infiltrations verticales et
d’autre part que chaque colonne peut être séparée en deux zones distinctes par un front raide.
En faisant tendre l’épaisseur du front vers zéro, on obtient une discontinuité dans la teneur en
eau et dans la perméabilité. La vitesse d’infiltration v de l’eau dans une colonne de sol est alors
simplement donnée par la loi de Darcy,
ψ 0 − ψ f + zf
v = Kf ,
zf
Le paramètre ψe est positif ou nul et désigne la pression d’entrée d’air. Une conséquence des
relations (1.7) et (1.8) est un changement de la nature mathématique de l’équation de Richards
en fonction de l’état de saturation du sol. Quand le milieu est insaturé (ψ < −ψe ), les propriétés
hydrodynamiques dépendent de la charge et l’équation (1.4) est de type parabolique non linéaire.
En revanche, lorsque le milieu est saturé (ψ ≥ −ψe ), les propriétés hydrodynamiques sont
constantes et l’équation (1.4) se simplifie en une équation elliptique linéaire.
Les deux modèles les plus utilisés en hydrologie pour représenter les propriétés hydrodyna-
miques des sols sont celui de Brooks–Corey [BC64] et celui de Van Genuchten [Gen80] dans
lesquels l’expression de la perméabilité relative de l’eau (l’indice w est dorénavant omis) est
obtenue à partir de la relation de Mualem [Mua76] qui exprime cette quantité en fonction de la
teneur en eau réduite θ̃ ∈ [0, 1] sous la forme
R θ̃ !2
l dτ /ψ(τ )
kr(θ̃) = θ̃ R01 ,
0 dτ /ψ(τ )
où θ̃ 7→ ψ(θ̃) est la fonction réciproque de ψ 7→ θ̃(ψ), l est un paramètre traduisant la connectivité
des pores souvent pris égal à 1/2 et θ̃ représente la teneur en eau réduite définie par
θ(ψ) − θr
θ̃(ψ) = ,
θs − θr
avec θr et θs désignant respectivement la teneur en eau résiduelle et la teneur en eau à saturation.
La quantité θr est la valeur limite de la teneur en eau lorsque ψ 7→ −∞. Cette valeur est non-nulle
pour des raisons analogues à celles impliquant que la teneur en eau à saturation est strictement
inférieure à la porosité (θs < ϕ).
Modèle de Brooks–Corey
−λ −n
|ψ| |ψ|
θ̃(ψ) = et K(ψ) = Ks ,
ψe ψe
Etant donné les expressions des propriétés hydrodynamiques des sols proposées par les
modèles, la non-linéarité de l’équation de Richards est dûe intrinsèquement aux variations de
la teneur en eau et de la conductivité hydraulique en fonction de la charge hydraulique (voir
la figure 1.3 et les courbes représentant les propriétés hydrodynamiques des sols utilisés dans
les cas tests). Des estimations de la conductivité hydraulique à saturation en fonction du type
de sol sont indiquées dans le tableau 1.1. Les valeurs sont reprises des travaux de Carsel et
Parrish [CP88].
1.2. Modélisation des écoulements en milieu poreux 11
θ θ
θs θs
drainage drainage
rétention rétention
θr ψ θr ψ
∂t h + ∂x q = 0, (1.9)
q 2 gh2
∂t q + ∂ x + = gh(S − J), (1.10)
h 2
où g [m.s−2 ] représente la gravité, S [−] est la pente du fond et J [−] prend en compte les
frottements sur le fond et les berges. Le principe pour obtenir les équations de Saint–Venant est
d’effectuer une intégration suivant la verticale (pour la version bidimensionnelle) ou suivant la
section (pour la version monodimensionnelle) des équations de Navier–Stokes en supposant que
l’écoulement est horizontal (typiquement une pente inférieure à 10%), la pression est hydrosta-
tique, les variations de la surface libre sont faibles et la turbulence est négligeable. On suppose
également qu’il n’y a pas de transfert de masse à travers le fond et la surface, que la vitesse
verticale est nulle au fond et que la cote du fond est indépendante du temps. Une démonstration
complète est par exemple donnée par Gerbeau et Perthame [GP01], Hervouet [Her03, p. 35] et
Viollet et al. [VCEL03, p. 215].
Plusieurs formules empiriques existent pour estimer le terme de frottement J. Les plus
connues sont celles de Chézy, de Strickler et de Darcy–Weisbash. Nous utilisons ici la relation
de Manning–Strickler qui exprime le débit en fonction du coefficient de Strickler K [m1/3 .s−1 ],
de la hauteur d’eau h et du frottement J sous la forme
L’écoulement est ici supposé unidirectionnel de sorte que le débit est positif ou nul (q ≥ 0).
S = J.
1.4. Objectifs de la thèse 13
En introduisant la relation (1.12) dans l’équation (1.9), nous obtenons l’équation de l’onde
cinématique sous la forme d’une loi de conservation scalaire,
∂t h + ∂x ϕ(h, S) = 0, (1.13)
où ϕ représente une fonction de flux pour la hauteur d’eau h. Cette équation sera complétée
par des termes sources pour tenir compte des échanges de masse avec le sol et de la présence
éventuelle d’un apport de masse par la pluie.
munie de conditions aux limites et d’une condition initiale. Ω représente un domaine en espace
et T désigne un temps fini de simulation. De nombreuses méthodes peuvent être employées pour
la discrétisation en espace de l’équation de Richards. Celia et Bouloutas [CB90], Woodward
et Dawson [WD00], Jones et Woodward [JW00] utilisent les différences finies (DF). Manzini
et Ferraris [MF04] choisissent les volumes finis (VF). Celia et Bouloutas [CB90], Lehmann et
Ackerer [LA98], Kavanagh et al. [KKB+ 02] travaillent avec les éléments finis (EF). Diersch et
Perrochet [DP99], Knabner et Schneid [KS02], Bause et Knabner [BK04] utilisent les éléments
finis mixtes (EFM). D’autres techniques ont été testées plus récemment, les éléments finis multi-
échelles pour des applications en milieu hétérogène par He et Ren [HR06] et une méthode
de Lattice Boltzmann par Ginzburg, Carlier et Kao [GCK04]. Pour étudier la concentration
d’une espèce chimique dans le sol, un couplage de l’équation de Richards avec une équation
de transport a aussi été effectué par Gabbouhy et al. [GMMS02] et par Diaw, Lehmann et
Ackerer [DLA01]. Les méthodes de discrétisation en espace sont souvent combinées avec un
schéma d’Euler implicite pour traiter le terme instationnaire.
14 Chapitre 1. Introduction
Les méthodes de Galerkin discontinues (DG) constituent une approche différente pour la
discrétisation en espace et offrent de nombreux avantages puisqu’elles sont localement conserva-
tives (comme les VF et les EFM), permettent l’utilisation de polynômes d’interpolation d’ordre
élevé (comme les EF et les EFM) éventuellement variable d’un élément à un autre (p-raffinement)
et offrent une grande flexibilité dans l’usage de maillages non-coı̈ncidants (comme les VF),
facilitant le h-raffinement. Plusieurs méthodes DG peuvent être utilisées pour l’équation de
Richards ainsi que pour les écoulements diphasiques en milieux poreux. La méthode « Local
Discontinuous Galerkin » (plus connue sous l’abréviation LDG en anglais) est considérée par
Fagherazzi et al. [FFRH04] et par Bastian et al. [BIR+ 07]. Les méthodes « non-symmetric
and symmetric interior penalty » (désignées respectivement par les acronymes anglais NIPG et
SIPG), développées par Baker [Bak77], Wheeler [Whe78], Arnold [Arn82], Oden et al. [OBB98],
Rivière, Wheeler et Girault [RWG99], Arnold et al. [ABCM01] et Ern et Guermond [EG06], sont
considérées par Klieber et Rivière [KR06] et par Bastian et Rivière [BR04] pour les écoulements
en milieu poreux. Dans cette thèse, nous choisissons la méthode SIPG puisqu’elle préserve la
symétrie dans l’opérateur de diffusion au niveau discret. Concernant la discrétisation temporelle
associée aux méthodes DG, les schémas les plus souvent rencontrés sont ceux de Runge–Kutta
(RK) explicites, développés par Cockburn et Shu [CS98a] pour les systèmes hyperboliques,
et parfois diagonalement implicites [Bas03]. Nous proposons d’utiliser plutôt les formules de
différentiation rétrograde (BDF pour « backward differentiation formulae ») dans lesquelles un
ordre élevé peut être obtenu sans condition de stabilité dans le cas linéaire lorsque les valeurs
propres de l’opérateur discret sont réelles négatives. Cela est notre cas si l’on suppose que la
conductivité hydraulique ne varie pas trop fortement dans le sol. De plus, les BDF permettent de
s’affranchir de la condition CFL classiquement recontrée avec les schémas explicites (cette condi-
tion pouvant être très restrictive quand la diffusion est dominante) et sont moins coûteuses que
les schémas RK diagonalement implicites notamment pour les problèmes non-linéaires. Lorsque
des polynômes de degré p sont utilisés dans les méthodes DG, une BDF d’ordre p + 1 est utilisée
en temps.
Nous effectuons ensuite un couplage entre l’équation de Richards (1.14) et l’équation de l’onde
cinématique (1.13) dans le cadre d’un ruissellement surfacique sur un sol insaturé. Le couplage
se produit sur une interface I qui désigne la partie de la frontière du domaine Ω représentant
la surface supérieure du sol. L’approche que nous choisissons repose sur l’égalité des flux d’eau
échangés et la continuité de la pression à l’interface. Cela signifie d’une part que la charge hy-
draulique du sol est égale à la hauteur d’eau ruisselante sur le sol et d’autre part que la vitesse
normale de l’écoulement souterrain est introduite comme terme source dans l’équation de conser-
vation de la masse des écoulements superficiels. Ces conditions de couplage sont généralement
valables quand l’écoulement superficiel est produit par une exfiltration d’eau contenue dans le
sol, un équilibre des flux d’eau et des pressions étant alors attendu. Parmi les travaux basés sur
cette méthode, nous pouvons citer Esteves et al. [EFGV00] couplant un écoulement surfacique
bidimensionnel avec des colonnes de sols verticales en utilisant le modèle de Green–Ampt, Singh
et Bhallamudi [SB98] couplant un écoulement surfacique unidimensionnel avec des colonnes
de sols verticales, Kollet et Maxwell [KM06] ainsi que Beaugendre et al. [BEE+ 06] couplant
l’équation de Richards en deux dimensions d’espace avec l’approximation de l’onde diffusive ou
cinématique en une dimension, et Dawson [Daw06] couplant l’équation de Darcy en deux dimen-
sions avec les équations de Saint–Venant en une dimension. D’autres approches existent pour
coupler les écoulements superficiels avec les écoulements souterrains. Une approche utilisée pour
coupler notamment les écoulements de Darcy et de Stokes, suppose la continuité des vitesses nor-
males et une discontinuité des pressions résultant de la condition de Beavers–Joseph–Saffman
qui fut introduite par Beavers et Joseph [BJ67] ainsi que Saffman [Saf71]. Cette condition à
l’interface modélise l’adhérence des composantes tangentielles du fluide s’écoulant en surface.
1.5. Plan de la thèse 15
Jäger et Mikelic [JM00] l’utilisent par exemple dans le cadre d’écoulements visqueux laminaires.
Cette condition a été utilisée par Discacciati et al. [DMQ02, MQS03] dans le cadre d’une étude
mathématique et numérique du couplage de l’écoulement de Darcy avec un modèle de Saint–
Venant non hydrostatique. Une autre approche rencontrée en hydrologie numérique et utilisée
entre autres par VanderKwaak et Loague [VL01] considère que la pression est discontinue et
évalue le flux à l’interface comme le produit de la différence de pression et d’un coefficient
empirique d’échange qui dépend de la perméabilité relative du sol.
Dans le cadre de notre travail, deux termes sources supplémentaires sont présents dans
l’équation (1.13) : l’un provient de la pluie et l’autre provient des interactions avec les
écoulements souterrains. Nous résolvons ainsi l’équation
munie de conditions aux limites et d’une condition initiale et dans laquelle v(ψ) représente le
terme de couplage, vr est l’intensité de la pluie et nΩ désigne la normale unitaire sortante à Ω.
Notons que la quantité v(ψ) · nΩ est positive quand de l’eau s’exfiltre du sol et négative quand
de l’eau s’y infiltre. L’équation de l’onde cinématique est discrétisée par un schéma de Godunov
avec un pas de temps choisi éventuellement plus petit que celui utilisé dans la discrétisation de
l’équation de Richards. Ce choix permet de respecter la condition CFL du schéma explicite sans
diminuer l’efficacité de la simulation numérique dans le sol. Nous avons choisi les volumes finis
pour faciliter la mise en œuvre du couplage mais une méthode DG (d’ordre supérieur ou égal
à un) avec des flux de Godunov est également possible. Cela nécessite néanmoins l’utilisation
de limiteurs de pente qui diminuent les oscillations et assurent la positivité de la hauteur d’eau
en particulier en présence de lits secs.
Les algorithmes de couplage sont conçus pour vérifier deux critères. Le premier est de sa-
tisfaire aussi précisément que possible les conditions de couplage qui imposent des contraintes
spécifiques de type égalité et inégalité sur la pression et la vitesse normale, comme pour les
conditions aux limites rencontrées dans les problèmes de Signorini. Le second est d’assurer la
conservation de la masse totale du système couplé. En effet, bien que des schémas conserva-
tifs soient utilisés pour les écoulements souterrains et pour les écoulements superficiels, le flux
à l’interface doit être convenablement adapté lorsque des méthodes à plusieurs pas de temps
(comme les BDF) sont choisies. Nous détaillons le choix du flux pour la BDF1 et la BDF2 mais
l’extension aux schémas d’ordre supérieur s’effectue sans difficulté.
doivent être déterminées pendant la simulation. Nous commençons par formaliser le problème
continu en indiquant les inégalités sur la charge hydraulique et la vitesse normale à vérifier sur
la frontière libre. Nous présentons ensuite l’algorithme implémenté pour satisfaire les contraintes
ainsi que les propriétés sur la conservation de la masse du système étudié. Enfin, un cas test avec
affleurement de nappe et un cas test avec un drain dans le sol permettent d’évaluer l’algorithme
proposé pour traiter les frontières libres.
Le chapitre 4 est consacré au couplage entre les écoulements souterrains et surfaciques.
Le système couplé consiste à résoudre l’équation de Richards dans un domaine Ω et l’équation
de l’onde cinématique sur une partie I de la frontière de Ω en respectant certaines contraintes
physiques. Contrairement au chapitre précédent, où les valeurs à imposer sur le partitionnement
de la frontière libre sont connues, les données dans les conditions aux limites de type Dirichlet
et de type Neumann sur la zone de couplage sont liées à la résolution de l’équation de l’onde
cinématique. La discrétisation par un schéma de Godunov de l’équation de l’onde cinématique est
décrite ainsi que les algorithmes permettant de satisfaire les contraintes physiques de couplage et
de conserver la masse totale du système. Notre algorithme conçu avec la BDF2 est évalué dans
quatre cas tests aux configurations variées. Nous présentons également une application concrète
en hydrologie dans le cadre d’une collaboration avec l’Institut de recherche pour l’ingénierie de
l’agriculture et de l’environnement (désigné par l’acronyme CEMAGREF). Le but de l’étude est
de prévoir le comportement hydrologique d’un système correspondant à une partie d’un bassin
versant expérimental situé au nord de Coulommier.
L’annexe A présente la discrétisation des équations de Saint–Venant par une méthode de
volumes finis. La modélisation des écoulements superficiels par ces équations permet de traiter
un plus grand nombre de cas tests qu’avec l’approximation de l’onde cinématique. Un cas test
de transfert entre sillons est présenté.
L’annexe B expose le logiciel conçu, developpé, testé et mis en œuvre au cours de cette thèse.
Nous l’avons appelé R+R DG en référence à Richards et Ruissellement par une méthode DG.
Nous en indiquons l’organisation générale et précisons l’assemblage des matrices élémentaires et
du second membre. Les éléments spécifiques à l’implémentation des conditions aux limites pour
le couplage sont aussi détaillés.
Chapitre 2
Nous présentons dans ce chapitre la méthode numérique utilisée pour résoudre l’équation
de Richards lorsque le milieu est variablement saturé (c’est-à-dire pouvant être aussi bien in-
saturé que saturé). Notre méthode est basée sur la discrétisation de l’équation de Richards
instationnaire et reste valable lorsque certaines parties du sol sont saturées puisqu’elle conduit
à la discrétisation de l’équation de Darcy.
Ce chapitre comporte deux parties. Tout d’abord, nous détaillons la discrétisation de
l’équation de Richards en exposant la méthode de Galerkine discontinue SIPG en espace,
l’intégration temporelle par les BDF et la linéarisation utilisée pour la résolution du problème
non-linéaire. Nous précisons aussi l’initialisation du schéma en temps par un schéma de Runge–
Kutta diagonalement implicite ou par le schéma de Crank–Nicolson. Dans la deuxième partie,
nous validons la méthode SIPG et son implémentation dans le logiciel R+R DG en réalisant plu-
sieurs cas tests. Deux cas tests analytiques permettent de vérifier les ordres de convergence en
milieu insaturé puis en milieu partiellement saturé. Plusieurs cas tests de colonne d’infiltration
permettent ensuite de tester et de confirmer la robustesse de notre méthode.
Les fonctions ψD et vN sont des données pour les conditions aux limites, imposées respectivement
sur ∂ΩD et ∂ΩN , la fonction ψ 0 est la condition initiale et T est un temps fini de simulation.
Dans le système (2.1), la charge hydraulique ψ est une quantité scalaire et la vitesse v(ψ) est
une quantité vectorielle. En posant u(ψ) = v(ψ) + K(ψ)ez , le système (2.1) est équivalent à
∂t [θ(ψ)] + ∇ · u(ψ) = ∇ · (K(ψ)ez ) dans Ω × [0, T ],
u(ψ) = −K(ψ)∇ψ dans Ω × [0, T ],
ψ = ψD sur ∂ΩD × [0, T ], (2.2)
N × [0, T ],
u(ψ) · n Ω = v N + K(ψ)ez · n Ω sur ∂Ω
ψ(., 0) = ψ 0 dans ∂Ω,
où ez = (0, 0, 1) désigne le vecteur vertical unitaire orienté vers le haut du repère cartésien.
dans lequel Pp (τ ) est l’ensemble des polynômes de degré inférieur ou égal à p sur un élément
quelconque τ . Les fonctions de l’espace Vh n’étant pas continues, cela permet de choisir des
fonctions de base dont le support correspond à un élément du maillage.
L’ensemble Fh des faces du maillage est partitionné en Fhi ∪ FhD ∪ FhN où Fhi est l’ensemble
des faces internes, FhD regroupe l’ensemble des faces sur lesquelles est imposée une condition
de Dirichlet et FhN désigne l’ensemble des faces sur lesquelles est imposée une condition de
Neumann. Nous supposons qu’un seul type de condition aux limites est imposée sur une face,
c’est-à-dire que ∀F ⊂ ∂Ω, F ⊂ ∂ΩD ou F ⊂ ∂ΩN . Pour une face interne F ∈ Fhi , il existe τ + et
τ − dans Th tels que F = ∂τ + ∩∂τ − et nous définissons l’opérateur de moyenne {}F et l’opérateur
de saut [[]]F de la façon suivante : pour une fonction ξ pouvant prendre deux valeurs sur F ,
def 1 + def
{ξ}F = (ξ + ξ − ) et [[ξ]]F = ξ − − ξ + ,
2
avec ξ ± = ξ|τ ± . Pour une fonction vectorielle, les opérateurs de moyenne et de saut sont définis
composante par composante. Nous désignons par nF la normale unitaire à F dirigée de τ −
vers τ + . Notons que le signe du saut est arbitraire mais sans conséquence pour la suite puisque
les sauts seront toujours multipliés par la normale nF .
Les notations ψh et uh désignent respectivement les approximations discrètes en espace et
continues en temps de ψ et u. Pour obtenir la formulation discrète du système (2.2), nous
constatons que la solution exacte (u, ψ) vérifie
Z Z
d
∀τ ∈ Th, ∀φ̃ ∈ [Pp (τ )] , u · φ̃ = − K(ψ)∇ψ · φ̃,
τ
Z Zτ Z (2.3)
∀τ ∈ Th, ∀φ ∈ Pp (τ ), ∂t [θ(ψ)]φ + ∇·u φ = ∇ · (K(ψ)ez )φ,
τ τ τ
Une intégration par parties du second membre de la première équation du système (2.3) donne
Z Z Z
∀τ ∈ Th, ∀φ̃ ∈ [Pp (τ )]d , u · φ̃ = ψ∇ · (K(ψ)φ̃) − K(ψ) φ̃ · nτ ψ.
τ τ ∂τ
2.1. Discrétisation de l’équation de Richards 19
Il n’est pas possible de remplacer directement (u, ψ) par (uh, ψh) dans cette équation car les
fonctions discrètes ne sont pas univaluées sur les interfaces. Nous introduisons dans l’intégrale
de bord un flux numérique scalaire Fψ (ψh) et écrivons au niveau discret
Z Z Z
∀τ ∈ Th, ∀φ̃ ∈ [Pp (τ )]d , uh · φ̃ = ψh∇ · (K(ψh)φ̃) − K(ψh|τ ) φ̃ · nτ Fψ (ψh). (2.4)
τ τ ∂τ
Le flux numérique Fψ (ψh) défini sur les faces du maillage est associé à la charge hydraulique ψ
et varie suivant que la face est interne, de bord avec une condition de Dirichlet ou de bord avec
une condition de Neumann,
i
{ψh}F si F ∈ Fh,
def
∀F ∈ Fh, ∀ψh ∈ Vh, Fψ (ψh)|F = ψD si F ∈ FhD ,
ψh si F ∈ FhN .
Notons que nous utilisons la donnée de la condition aux limites de Dirichlet pour la définition
de ce flux sur FhD . Après une contre-intégration par parties, l’équation (2.4) devient
Z Z Z
d
∀τ ∈ Th, ∀φ̃ ∈ [Pp (τ )] , uh · φ̃ = − K(ψh)∇ψh · φ̃ − K(ψh|τ ) φ̃ · nτ (Fψ (ψh) − ψh|τ ).
τ τ ∂τ
(2.5)
Nous effectuons maintenant une intégration par parties dans le deuxième terme de la seconde
équation du système (2.3),
Z Z Z Z
∀τ ∈ Th, ∀φ ∈ Pp (τ ), ∂t [θ(ψ)]φ − u · ∇φ + u · nτ φ = ∇ · (K(ψ)ez )φ.
τ τ ∂τ τ
Pour les mêmes raisons que précédemment, nous introduisons dans l’intégrale de bord un flux
numérique vectoriel Fu (ψh) associé à la grandeur u et écrivons au niveau discret
Z Z Z Z
∀τ ∈ Th, ∀φ ∈ Pp (τ ), ∂t [θ(ψh)]φ− uh ·∇φ+ Fu (ψh)·nτ φ = ∇·(K(ψh)ez )φ. (2.6)
τ τ ∂τ τ
Le flux numérique Fu (ψh) défini sur les faces du maillage varie suivant que la face est interne,
de bord avec une condition de Dirichlet ou de bord avec une condition de Neumann,
−1 i
−{K(ψh)∇ψh}F + ηKs dF [[ψh]]F nF si F ∈ Fh,
def
∀F ∈ Fh, ∀ψh ∈ Vh, Fu (ψh)|F = − K(ψh)∇ψh + ηKs d−1 F (ψh − ψD )nΩ si F ∈ Fh ,
D
vN nτ + K(ψh)ez si F ∈ FhN ,
où η est un paramètre positif qui sert à la fois à pénaliser les sauts aux interfaces par moindres
carrés et à imposer de manière faible les conditions aux limites de Dirichlet sur les faces dans FhD .
Ce paramètre η doit être plus grand qu’une valeur minimale dépendant de la régularité des
maillages {Th}h>0 . Notons également que nous utilisons la donnée de la condition aux limites de
Neumann sur les faces de FhN . Par ailleurs, Ks est la conductivité à saturation et dF une échelle
de longueur de la face F qui sera ici évaluée en utilisant les diamètres des éléments pour lesquels
F est une face,
(
def max(dτ + , dτ − ), F = ∂τ + ∩ ∂τ − , si F ∈ Fhi ,
∀F ∈ Fh, dF =
dτ , F ⊂ ∂τ, si F ∈ FhD ,
où dτ désigne la plus grande arête de la maille τ. Pour un milieu poreux ayant une conductivité
variable (comme c’est le cas dans les écoulements variablement saturés puisque la conductivité
20 Chapitre 2. Ecoulements en milieu poreux variablement saturé
Nous avons intégré par parties le membre de droite de l’équation (2.7). Cette manipulation
s’avère importante pour améliorer la robustesse de la méthode dans certains cas ; nous y revien-
b ψ (ψh)
drons à la fin de ce chapitre (section 2.2.5). Dans l’expression de aτ (ψh, φ), la quantité F
est définie par
F (ψ ) si F ∈ Fhi ,
ψ h
∀F ∈ Fh, ∀ψh ∈ Vh, b ψ (ψh)|F def
F = 0 si F ∈ FhD ,
F (ψ ) si F ∈ F N ,
ψ h h
Dans l’expression de bτ (t, ψh, φ), la quantité b̃τ (t, ψh, φ) est définie par
Z Z
def
b̃τ (t, ψh, φ) = − K(ψh) ∇φ · nΩ + ηKs d−1
F φ ψD − (vN + K(ψh)ez · nΩ )φ.
∂τ ∩FhD ∂τ ∩FhN
La dépendance par rapport au temps de b̃τ (t, ψh, φ) vient de la dépendance par rapport au temps
des données des conditions aux limites ψD et vN . Cette dépendance explicite nous sera utile dans
la description des schémas DIRK3 et de Crank–Nicolson aux sections 2.1.4 et 2.1.5. Enfin, en
sommant la forme locale aτ (ψ, φ) sur l’ensemble des éléments du maillage, nous obtenons la
2.1. Discrétisation de l’équation de Richards 21
def
X
ah(ψh, φ) = aτ (ψh, φ)
τ ∈Th
XZ
= K(ψh)∇ψh · ∇φ
τ ∈Th τ
X Z
− {K(ψh)∇φ}F [[ψh]]F · nF + {K(ψh)∇ψh}F [[φ]]F · nF − ηKs d−1
F [[ψh]]F [[φ]]F
F
F ∈Fhi
X Z
− K(ψh)∇φ ψh · nΩ + K(ψh)∇ψh φ · nΩ − ηKs d−1
F ψhφ . (2.11)
F
F ∈FhD
De manière analogue à l’approximation du Laplacien par la méthode SIPG (voir par exemple
Arnold et al. [ABCM01]), le paramètre η doit être choisi suffisamment grand pour assurer que
la forme globale ah soit coercive, c’est-à-dire qu’il existe α > 0 tel que
XZ X Z X Z
∀φh ∈ Vh, ah(φh, φh) ≥ α K(φh)|∇φh|2 + Ks d−1
F [[φh]]2F + Ks d−1
F φ2h .
τ ∈Th τ F F
F ∈Fhi F ∈FhD
(2.12)
Le seuil minimal de η dépend de la régularité des maillages {Th}h>0 .
n 1 3 n 1
∂t χ = χ − 2χn−1 + χn−2 + O(δt2 ), (2.15)
δt 2 2
ainsi que la BDF d’ordre 3 pour laquelle α03 = 11/6, α13 = −3, α23 = 3/2 et α33 = −1/3,
n 1 11 n 3 1
∂t χ = χ − 3χn−1 + χn−2 − χn−3 + O(δt3 ). (2.16)
δt 6 2 3
22 Chapitre 2. Ecoulements en milieu poreux variablement saturé
Pour les premiers pas de temps, des BDF d’ordres inférieurs ou des schémas implicites à pas
unique sont utilisés (voir section 2.1.4).
2.1.3 Linéarisation
L’équation non-linéaire (2.17) est résolue itérativement par l’algorithme 1 présenté ci-après.
Nous ajoutons un argument à la forme aτ ainsi qu’au flux F b u pour indiquer le traitement des
termes non-linéaires. Ainsi, pour tout (ζ, ψ, φ) ∈ Vh × Vh × Pp (τ ), nous posons
Z Z Z
def b ψ (ψ) − ψ) + b u (ζ, ψ) · nτ φ,
aτ (ζ, ψ, φ) =
b K(ζ)∇ψ · ∇φ + K(ζ|τ ) ∇φ · nτ (F F
τ ∂τ ∂τ
Nous rappelons que la quantité C(ψ) est la dérivée de la teneur en eau par rapport à la charge
hydraulique et s’appelle la capacité capillaire. En posant δψhn,m = ψhn,m+1 −ψhn,m et en définissant
cτ par Z
def
cτ (ζ, ψ, φ) = C(ζ)ψφ,
τ
α0q
∀τ ∈ Th, ∀φ ∈ Pp (τ ), cτ (ψhn,m , δψhn,m , φ) + b
aτ (ψhn,m , δψhn,m , φ) = bτ (nδt, ψhn,m , φ)
δt
q Z Z
X αrq αq
− θ(ψhn−r )φ − 0 θ(ψhn,m )φ − b aτ (ψhn,m , ψhn,m , φ). (2.19)
δt τ δt τ
r=1
Il est important de préciser que notre méthode permet de traiter automatiquement la saturation
complète du milieu poreux. En effet, dans ce cas, la teneur en eau est constante et sa dérivée
2.1. Discrétisation de l’équation de Richards 23
par rapport à la charge hydraulique est nulle, ce qui modifie l’équation (2.19) puisque le terme
cτ est nul,
L’équation (2.20) est vérifiée quand le milieu est saturé à l’instant nδt mais insaturé aux instants
antérieurs. Dans le cas où la saturation est obtenue pendant suffisamment de pas de temps,
les termes du second membre provenant de la BDF s’annulent (car la somme des coefficients
{αrq }0≤r≤q est nulle) et l’équation (2.20) se simplifie,
Algorithme 1 Algorithme itératif à chaque pas de temps pour résoudre l’équation de Richards
les βjr et les cr étant des coefficients donnés dans le tableau de Butcher ci-après. La nature impli-
(j)
cite du schéma provient de la présence de ψh dans le second membre de l’étape j. La résolution
(j)
de l’équation (2.23) s’effectue en décomposant le terme fτ (n − 1 + cj )δt, ψh , φ ,
∀ 1 ≤ j ≤ 3, ∀τ ∈ Th, ∀φ ∈ Pp (τ ),
Z Z j−1
X
(j) (j) (r)
θ(ψh )φ + δtβjj aτ (ψh , φ) = θ(ψhn−1 )φ + δt βjr fτ (n − 1 + cr )δt, ψh , φ
τ τ r=1
(j)
+ δtβjj bτ (n − 1 + cj )δt, ψh , φ .
En divisant cette équation par δtβjj , nous obtenons une équation similaire à l’équation (2.17)
∀ 1 ≤ j ≤ 3, ∀τ ∈ Th, ∀φ ∈ Pp (τ ),
−1 Z −1 Z
βjj (j) (j) (j) βjj
θ(ψh )φ + aτ (ψh , φ) = bτ (n − 1 + cj )δt, ψh , φ + θ(ψhn−1 )φ
δt τ δt τ
j−1
X βjr (r)
+ fτ (n − 1 + cr )δt, ψh , φ , (2.24)
βjj
r=1
où les coefficients réels γ1 , γ2 et γ3 ont une somme égale à 1 et sont donnés dans le tableau 2.1
ci-dessous. Comme précédemment cette équation non-linéaire est résolue de façon itérative en
calculant des approximations successives ψhn,m de ψhn en utilisant la linéarisation suivante,
∀τ ∈ Th, ∀φ ∈ Pp (τ ), (2.25)
Z Z 3
X
(r)
θ(ψhn,m ) + C(ψhn,m )δψhn,m φ = θ(ψhn−1 )φ + δt γr fτ (n − 1 + cr )δt, ψh , φ ,
τ τ r=1
ou encore
∀τ ∈ Th, ∀φ ∈ Pp (τ ),
Z 3
X
(r)
cτ (ψhn,m , δψhn,m , φ) = θ(ψhn−1 ) − θ(ψhn,m ) φ + δt γr fτ (n − 1 + cr )δt, ψh , φ .
τ r=1
Les différents coefficients du schéma DIRK3 sont indiqués dans le tableau 2.1 en fonction
d’un paramètre
√ µ vérifiant l’équation
√ −3µ3 + 3µ + 1 = 0 √ et dont les trois racines sont
µ1 = 2/ 3 cos(π/18), µ2 = −2/ 3 cos(5π/18) et µ3 = −2/ 3 cos(7π/18). Le domaine de
stabilité de la méthode est infini si µ = µ1 alors qu’il est borné si µ = µ2 et µ = µ3 . Bien qu’avec
µ = µ3 , la méthode soit d’ordre 4, l’étude du domaine de stabilité en fonction de la valeur du
paramètre µ nous a amené à choisir la valeur µ1 .
Notons que la première étape du schéma DIRK3 conduit à résoudre trois systèmes linéaires
du même type que ceux rencontrés avec les méthodes BDF puisque l’équation (2.24) a la même
forme que l’équation (2.17). En revanche la seconde étape est différente puisqu’il n’y a pas
de terme implicite provenant de la discrétisation en espace dans l’équation (2.25). Cela rend
inutilisable le schéma DIRK3 quand une partie du domaine est saturée au début de la simulation.
C’est la raison pour laquelle nous avons également implémenté le schéma de Crank–Nicolson
qui est présenté à la section 2.1.5. Lors de l’initialisation, lorsque les conditions initiales et les
conditions aux limites sont incompatibles (comme c’est le cas notamment pour les cas tests
proposés par Haverkamp et Polmann et présentés dans les sections 2.2.3 et 2.2.4), nous avons
remarqué que les pas de temps du schéma DIRK3 doivent être plus petits que les pas de temps
des schémas BDF pour que l’algorithme 1 converge. Nous avons ainsi implémenté un algorithme
avec des sous-pas de temps adaptatifs pour le premier pas de temps. Nous notons δtp le p-ième
sous-pas de temps pour le schéma DIRK3. L’idée de base de cet algorithme est d’augmenter le
pas de temps δtp si le critère de convergence est vérifié et de le diminuer dans le cas contraire.
Cet algorithme est décrit ci-après dans le cas n = 1, ce qui correspond au premier pas de temps.
26 Chapitre 2. Ecoulements en milieu poreux variablement saturé
Donnons quelques précisions sur cet algorithme (voir également la figure 2.1 ci-après). Lors de
l’initialisation il est possible que le critère de convergence E ≤ ǫalg1 de l’algorithme 1 ne soit pas
vérifié quand le pas de temps δt est trop grand. Pour résoudre ce problème, nous introduisons un
nombre d’itérations maximales qui permet d’arrêter la boucle de l’algorithme 1 le cas échéant.
Le calcul est alors repris avec un pas de temps dix fois plus petit. La valeur du premier pas
de temps δt1 résulte des observations faites au cours des différentes simulations. Lorsqu’il y
a compatibilité entre la condition initiale et les conditions aux limites, le premier sous-pas de
temps δt1 peut être pris égal à δt. En revanche lorsqu’il y a discontinuité entre les différentes
conditions, il est nécessaire de prendre le premier sous-pas de temps δt1 plus petit que δt (δt1
pouvant être de l’ordre de 10−5 δt dans les cas les plus défavorables). Le calcul de la solution
à un sous-pas de temps δtp avec le schéma DIRK3 nécessite la résolution d’au moins quatre
systèmes non-linéaires. Pour avancer d’un pas de temps, le nombre de systèmes à résoudre est
0,(p)
ainsi égal à quatre fois le nombre de « solutions provisoires » ψh calculées. Par ailleurs dans
les arguments de la fonction « min » utilisés pour le calcul de δtp+1 , le coefficient 2 permet une
augmentation progressive des sous-pas de temps, le deuxième argument δt/4 impose un minimum
de quatre sous-pas de temps par pas de temps δt et le dernier argument trestant permet de finir
l’algorithme sur le calcul de la solution au premier temps discret t1 = δt. La figure 2.1 illustre
l’initialisation du schéma d’intégration temporelle avec la BDF2 : ψh0 est donné, ψh1 est déterminé
par l’algorithme 2 à partir de ψh0 en utilisant plusieurs fois le schéma DIRK3 et ψh2 est calculé à
partir de ψh0 et ψh1 en utilisant la BDF2.
2.1. Discrétisation de l’équation de Richards 27
0 δt 2δt
δt1 δt2 δt3 δt4 δt5
| × × × × × | |
ψh0,1 ψh0,2 ψh0,3 ψh0,4 ψh0,5
ψh0 ψh0,6 = ψh1 ψh2
ou encore
Z
δt
∀τ ∈ Th, ∀φ ∈ Pp (τ ), cτ (ψhn,m , δψhn,m , φ) + aτ (ψhn,m , δψhn,m , φ) =
b θ(ψhn−1 ) − θ(ψhn,m ) φ
2 τ
δt
+ bτ (nδt, ψhn,m , φ) − b
aτ (ψhn,m , ψhn,m , φ) + fτ ((n − 1)δt, ψhn−1 , φ) ,
2
que nous résolvons par l’algorithme itératif 1. Plus généralement, pour initialiser la BDFq si
une zone du sol est saturée, il faut implémenter un schéma d’ordre q dans lequel b aτ dépend
de la quantité δψhn,m . Nous n’avons pas utilisé la BDF3 excepté pour les tests de convergence
(section 2.2.2), mais nous avons remarqué que la combinaison du schéma de Milne–Simpson
d’ordre 4 (MS4) et du schéma de Crank–Nicolson est d’ordre 3. Le schéma MS4 consiste à
résoudre ∀τ ∈ Th, ∀φ ∈ Pp (τ ),
Z Z
δt
θ(ψhn )φ = θ(ψhn−2 )φ + fτ (nδt, ψhn , φ) + 4fτ ((n − 1)δt, ψhn−1 , φ) + fτ ((n − 2)δt, ψhn−2 , φ) .
τ τ 3
Ainsi une façon d’initialiser la BDF3 lorsque le milieu est partiellement saturé est de déterminer
1/2
ψh1 en combinant ces deux schémas (ψh étant calculé par CN2 et ψh1 par MS4) puis de déterminer
ψh2 par MS4.
28 Chapitre 2. Ecoulements en milieu poreux variablement saturé
2.2 Résultats
Avant de présenter les résultats, nous commençons par décrire les maillages utilisés.
Nous avons construit un indicateur qui mesure la qualité des triangles de chaque maillage.
De plus nous avons utilisé une renumérotation des triangles pour réduire la largeur de bande
des matrices. Ensuite deux cas tests analytiques valident les ordres de convergence de notre
méthode. Trois cas tests de colonne d’infiltration unidimensionnelle (traités en dimension deux)
sont aussi réalisés pour évaluer la robustesse de notre approche par rapport à l’état hydrologique
et aux propriétés hydrodynamiques du sol.
2.2.1 Maillages
Caractéristiques des maillages
Dans les cas tests réalisés nous utilisons des maillages non structurés générés par le logiciel
c
FreeFEM
développé par Danaila, Hecht et Pironneau [Fre]. Pour déterminer la qualité d’un
maillage Th, nous calculons pour chaque triangle τ ∈ Th l’indicateur Qτ suivant
def √ ρτ
Qτ = 2 3 ,
dτ
dans lequel ρτ désigne le rayon du cercle inscrit et dτ est la plus grande arête du triangle τ .
Cet indicateur est égal à 1 pour un triangle équilatéral et tend vers zéro lorsque le triangle
est de plus en plus aplati. Comme nous effectuons en dimension deux des simulations sur des
colonnes d’infiltration (cas réalisables en dimension un), le domaine utilisé est un rectangle.
Sa largeur est égale au cinquième de la longueur, l’écoulement se faisant suivant la longueur.
La taille du rectangle initial est [0, 0.2] × [0, 1] et un coefficient de dilatation permet de trans-
former ces dimensions de façon isotrope pour obtenir la longueur souhaitée suivant les cas tests.
Pour les six maillages nous précisons dans le tableau 2.2 le nombre de nœuds Nn , le nombre de
triangles Nt , le diamètre minimal hmin , moyen hmoy et maximal hmax des mailles ainsi que la
qualité minimale Qmin , moyenne Qmoy et maximale Qmax .
Tab. 2.2: Caractéristiques des maillages du rectangle [0, 0.2] × [0, 1].
c
Fig. 2.2: Maillages triangulaires non-structurés générés par le logiciel FreeFEM
.
bande initiale lbinit et la largeur de bande après renumérotation lbrenu . La figure 2.3 montre, pour
le maillage M4 , la position des blocs non nuls dans la matrice avec et sans la renumérotation.
Comme attendu nous pouvons observer sur cette figure la structure symétrique des matrices.
M1 M2 M3 M4 M5 M6
lbinit 8 31 195 536 3121 12695
lbrenu 4 5 8 19 32 66
Un point particulièrement agréable des méthodes d’éléments finis discontinus est qu’elles ne
nécessitent pas de renumérotation quand on change le degré des polynômes utilisés. En effet, pour
les éléments finis discontinus, la renumérotation est valable quel que soit l’ordre d’interpolation
puisqu’elle porte sur la position des éléments les uns par rapport aux autres alors que, pour les
éléments finis continus, la renumérotation porte sur les degrés de liberté et change donc suivant
l’ordre utilisé.
Fig. 2.3: Influence de la renumérotation sur la position des blocs non nuls (maillage M4 ).
Il en est de même pour la BDF3 dans laquelle les deux premiers pas de temps sont calculés par
le schéma DIRK3 et une initialisation d’ordre 3 est réalisée lorsque n ≥ 4,
Nous utilisons une formule d’intégration numérique exacte pour les polynômes P5 . Quatre erreurs
sont calculées. Pour toute fonction χ ∈ L2 (Ω), nous désignons par ||χ||L2 (Ω) la norme L2 de χ
sur Ω. L’erreur relative L2 en espace et L∞ en temps pour la charge hydraulique ψ est définie
par
max ||ψ n − ψhn ||L2 (Ω)
1≤n≤NT
eψ,Ω = ,
max ||ψ n ||L2 (Ω)
1≤n≤NT
et son ordre de convergence théorique est p+1. Cela signifie qu’en décroissant la taille des mailles
et le pas de temps par 2, l’erreur est divisée par 2p+1 . L’erreur relative L2 en espace et L2 en
temps pour la vitesse v est définie par
||v − vh||L2 (Ω×[0,T ])
ev,Ω = ,
||v||L2 (Ω×[0,T ])
et son ordre de convergence théorique est p. Enfin, les deux dernières erreurs concernent res-
pectivement les sauts de ψh sur les faces internes et les faces de Dirichlet. Elles sont définies
par
Z T X 1
ηKs 2
2
eψ,F i = ||[[ψ − ψh ]]||L2 (F ) ,
h
0 i
dF
F ∈Fh
Z T X ηKs 1
2
2
eψ,F D = ||ψ − ψh ||L2 (F ) ,
h
0 D
dF
F ∈Fh
et leur ordre de convergence théorique est p. L’intégrale en temps est calculée par la méthode
des rectangles.
L’équation de Richards est résolue avec les conditions aux limites de Dirichlet et le terme source
nécessaire pour que ψ donnée par (2.26) soit solution. La teneur en eau et la conductivité
hydraulique sont données par
θs − θr Ks
θ(ψ) = + θr et K(ψ) = , (2.27)
1 + |α̃ψ|β 1 + |Ãψ|γ
Les tableaux 2.4 et 2.5 présentent les résultats obtenus avec le logiciel R+R DG sous raffinement
spatial et temporel. Les taux de convergence concordent avec ceux prédits par l’analyse d’erreur
a priori et qui ont été rappelés précédemment. Nous pouvons observer la superconvergence
de l’erreur eψ,F D . Les calculs sont effectués avec le paramètre de pénalisation η = 4 pour les
h
polynômes P1 et η = 12 pour les polynômes P2 . Pour les tests réalisés sur l’ensemble des maillages
(Mi )2≤i≤6 , nous donnons également la quantité adimensionnelle t∗CPU définie comme le rapport
entre le temps obtenu pour le calcul et le temps obtenu pour le calcul précédent (maillage plus
grossier et pas de temps deux fois plus important). Une valeur de 8 pour t∗CPU est raisonnable en
dimension deux d’espace puisqu’elle indique un coût linéaire en le nombre de degrés de liberté.
Les cas tests sont désormais réalisés uniquement avec la méthode P1 -BDF2. La frontière du
domaine est divisée en trois parties : I est la partie supérieure, W sont les parois latérales et
B désigne le fond du domaine (voir la figure 2.4). Ces trois lettres sont les premières lettres des
mots anglais interface, walls et bottom. Une pression de −61.5cm est imposée à l’instant initial
2.2. Résultats 33
Tab. 2.6: Résultats de convergence du cas test variablement saturé - méthode P1 –BDF2.
Tab. 2.7: Résultats de convergence du cas test variablement saturé - méthode P2 –BDF3.
34 Chapitre 2. Ecoulements en milieu poreux variablement saturé
ainsi que sur le fond pendant toute la durée de la simulation. De l’eau est injectée par le haut
du domaine en imposant une pression de −20.7cm et les parois latérales sont imperméables,
ψ 0 = −61.5cm dans Ω,
ψD = −61.5cm sur B × [0, T ],
ψD = −20.7cm sur I × [0, T ],
vN = 0 sur W × [0, T ].
ψD = −20.7cm
I
W
vN = 0 0.4m
W
(0, 0) • B
ψD = −61.5cm
Le sol est un sable dont les propriétés hydrodynamiques sont définies par les relations (2.27) et
(2.28). Comme le montrent les deux courbes de la figure 2.5, la teneur en eau et la conductivité
ne présentent pas de forts gradients.
−2
−4
x 10
0.3 1
Conductivité hydraulique K [cm.s−1 ]
0.25 0.8
Teneur en eau θ [−]
0.2 0.6
0.15 0.4
0.1 0.2
0.05 0
−80 −70 −60 −50 −40 −30 −20 −10 0 −80 −70 −60 −50 −40 −30 −20 −10 0
Charge hydraulique ψ [cm] Charge hydraulique ψ [cm]
Toutes les 120s, nous représentons sur la figure 2.6 la charge hydraulique moyenne suivant
la direction x définie par
Z
1 l
ψ¯h(z) = ψh(x, z)dx,
l 0
où l est la largeur du domaine.
2.2. Résultats 35
−40 −40
600s 600s
−50 480s −50 480s
360s 360s
240s 240s
−60 120s −60 120s
0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 30 35 40
z [cm] z [cm]
Nous rappelons que I, W et B désignent respectivement la partie supérieure, les parois latérales
et le fond du domaine. Une pression de -10m est imposée à l’instant initial ainsi que sur le fond
pendant toute la durée de la simulation. Une pression de -75cm (ce qui équivaut à une surpression
de 9.25m) est imposée en haut de la colonne et les parois latérales sont imperméables,
ψ 0 = −10m dans Ω,
ψD = −10m sur B × [0, T ],
ψD = −75cm sur I × [0, T ],
vN = 0 sur W × [0, T ].
Les propriétés hydrodynamiques du sol sont données par les relations de Van Genuchten
2
(θs − θr ) 1 − (ǫ|ψ|)n−1 (1 + (ǫ|ψ|)n )−m
θ(ψ) = m + θr et K(ψ) = Ks m , (2.30)
1 + (ǫ|ψ|)n 1 + (ǫ|ψ|)n 2
36 Chapitre 2. Ecoulements en milieu poreux variablement saturé
ψD = −75cm
I
W
vN = 0 1m
W
(0, 0) • B
ψD = −10m
−2
−5
x 10
0.4 10
Conductivité hydraulique K [cm.s−1 ]
0.35 8
Teneur en eau θ [−]
0.3
6
0.25
4
0.2
2
0.15
0
0.1
−1000 −800 −600 −400 −200 0 −1000 −800 −600 −400 −200 0
Charge hydraulique ψ [cm] Charge hydraulique ψ [cm]
Les résultats sont présentés avec le maillage M4 et un pas de temps δt = 1min. La valeur
moyenne de la charge ψ¯h est tracée sur la figure 2.9 à 12h, 24h, 36h et 48h. La partie gauche
de la figure 2.9 présente les résultats obtenus. Comme nous le voyons sur cette figure, la rai-
deur du front d’infiltration provoque l’apparition d’oscillations parasites. C’est la raison pour
laquelle nous avons implémenté le limiteur de pente proposé par Cockburn et Shu [CS98b] pour
les systèmes hyperboliques. La partie droite de la figure 2.9 montre que l’utilisation du limiteur
diminue considérablement les oscillations (courbes continues) même si cela engendre un retard
par rapport à la solution non limitée (courbes en traits pointillés, reprises de la figure gauche).
Notons que ce retard n’est pas incompatible avec la conservation de la masse lors de l’utilisa-
tion du limiteur puisque l’eau dans le sol est injectée en imposant une condition sur la charge
hydraulique.
2.2. Résultats 37
−4 −4
−6 −6
−8 −8
−10 −10
vN = 0 ψD = 0 vN = −0.5Ks
I I I
W W W
vN = 0 1m vN = 0 1m vN = 0 1m
W W W
CL 1 CL 2 CL 3
Les propriétés hydrodynamiques du sol sont données par les relations de Van Genuchten mo-
difiées
1/m )m 2
1 1 − (1 − (θ̃/β)
θ(ψ) = θ̃(θs − θr ) + θr et K(ψ) = Ks θ̃ 2 2 , (2.31)
1 − (1 − (1/β)1/m )m
où les deux quantités θ̃ et β sont données pour ψ ≤ −hs par
(1 + (ǫhs )n )m
θ̃ = et β = (1 + (ǫhs )n )m ,
(1 + (ǫ|ψ|)n )m
La figure 2.11, qui représente les propriétés hydrodynamiques pour hs = 10−3 cm en trait dis-
continu et pour hs = 2cm en trait continu, montre que la conductivité hydraulique est beaucoup
plus sensible que la teneur en eau à la variation de la hauteur minimum capillaire.
−5
x 10
0.39 6
Conductivité hydraulique K [cm.s−1 ]
5
Teneur en eau θ [−]
0.38 4
hs = 2cm
hs = 10−3 cm
3
0.37 2
hs = 2cm
1
hs = 10−3 cm
0.36 0
−100 −80 −60 −40 −20 0 20 −100 −80 −60 −40 −20 0 20
Charge hydraulique ψ [cm] Charge hydraulique ψ [cm]
hs = 10−3 cm hs = 2cm
2 2
0 0
Charge hydraulique ψ¯h [cm]
−4 −4
−6 −6
−8 −8
−10 −10
0.2j 0.6j 1j 0.2j 0.4j 0.6j 0.8j 1j
−12 −12
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100
z [cm] z [cm]
2 2
0 0
Charge hydraulique ψ¯h [cm]
−2 −2
−4 −4
−6 −6
−8 −8
−10 −10
0.5j 0.4j 0.3j 0.2j
0.5j 0.3j 0.1j 0.1j
−12 −12
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100
z [cm] z [cm]
2 2
0 0
Charge hydraulique ψ¯h [cm]
−2 −2
−4 −4 2j
−6 −6
−8 −8
−10 −10
2j 1.5j 1.5j 1j
1j 0.5j 0.1j 0.5j 0.1j
−12 −12
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100
z [cm] z [cm]
Les figures 2.12, 2.13 et 2.14 représentent les résultats obtenus avec le maillage M5 et un
pas de temps δt = 2min. Les valeurs moyennes de la charge hydraulique ψ¯h sont tracées en
fonction de la profondeur sur les figures de droite pour hs = 10−3 cm et sur les figures de gauche
pour hs = 2cm. Les courbes sont tracées toutes les 2h48min (=0.2j) pour la CL 1, toutes les
1h24min (=0.1j) pour la CL 2 et toutes les 6h (=0.5j) pour la CL 3.
Les résultats obtenus concordent avec ceux de Vogel, Van Genuchten et Cislerova [VGC01]
ainsi que ceux de Ginzburg, Carlier et Kao [GCK04] lorsque hs = 2cm. En revanche nous
avons observé que les évolutions du front d’infiltration sont fortement sensibles au paramètre hs
lorsque celui-ci est proche de zéro. Cela s’explique par les très fortes variations de la conductivité
hydraulique quand le milieu devient saturé si la hauteur minimum capillaire est nulle. De plus,
l’algorithme 1 ne converge pas (ou converge très lentement) si la divergence de la conductivité
apparaı̂t dans le second membre. Nous avons remarqué que l’intégration par parties de l’intégrale
surfacique contenant la divergence de la conductivité permet d’améliorer considérablement la
convergence, Z Z Z
∇ · (K(ψ)ez )φ = − K(ψ)ez ∇φ + K(ψ|τ )ez φ nτ .
τ τ ∂τ
Les résultats obtenus avec hs = 10−3 cm (valeur minimale de hs pour laquelle nos temps de calcul
sont raisonnables) s’éloignent un peu de ceux de Vogel et al. obtenus avec hs = 0 confirmant la
très forte sensibilité du système à la hauteur capillaire minimum quand celle-ci est faible.
Comme attendu pour les trois cas, l’infiltration est plus rapide avec hs = 2cm puisque la
conductivité hydraulique est supérieure au cas où hs = 10−3 (voir la figure 2.11). L’influence de la
gravité sur l’écoulement est également visible car l’infiltration par le haut (CL 2) est plus rapide
que l’infiltration par le bas (CL 1). Ces trois cas sont particulièrement intéressants car ils nous ont
permis de tester notre méthode, avec une zone dans le sol où la variation de la conductivité est
très grande, et de dégager une information importante concernant l’implémentation du second
membre que nous avons pris en compte dans le logiciel R+R DG.
Chapitre 3
Affleurement et drains
Nous avons présenté dans le chapitre précédent les écoulements en milieu poreux variablement
saturé, c’est-à-dire lorsque l’eau remplit partiellement ou complètement les pores du sous-sol,
avec des conditions aux limites déterminées de façon explicite sur la frontière du domaine.
Il existe cependant des situations dans lesquelles les conditions aux limites sont définies par des
inégalités et les morceaux de frontière, sur lesquels celles-ci s’imposent, ne sont pas explicitement
déterminés : on parle de problèmes à frontière libre. En hydrologie ce type de problème apparaı̂t
lorsqu’une nappe affleure par exemple. La détermination de la position de la nappe est un
problème de type obstacle régi par des conditions aux limites de Signorini. La modélisation des
drains fait aussi appel à ce type de conditions aux limites.
Ce chapitre est organisé en deux parties. La première partie expose les conditions de Signorini
et l’algorithme proposé pour résoudre les problèmes avec ce type de conditions aux limites.
La seconde partie présente les résultats obtenus dans les deux exemples évoqués précédemment :
d’abord la détermination de la position d’une nappe affleurante, puis la modélisation de drains.
Les exposants w et d ainsi que l’indice r font référence aux premières lettres des mots anglais
« wet », « dry » et « rain ». Par abus de langage nous appelons « partie sèche » la partie non
saturée et « partie mouillée » la partie saturée. Les partitions des ensembles I et D ne sont pas
connues a priori et dépendent du temps. Ainsi les deux doublets {I w,t , I d,t } et {Dw,t , Dd,t } sont
à déterminer en même temps que nous résolvons l’équation de Richards.
42 Chapitre 3. Affleurement et drains
vr
I
▽
ψ>0 ψ<0 ψ>0
W W nΩ
Ω
B D
L’affleurement d’une nappe dans une parcelle soumise à une pluie apparaı̂t lorsque le milieu
poreux n’est plus capable d’absorber l’eau provenant de cette pluie. Pour modéliser ce type de
phénomène, les conditions aux limites sont des conditions de Signorini (aussi appelées condition
de complémentarité) portant sur la charge hydraulique et la vitesse normale à la surface :
ψ ≤ 0,
v(ψ) · nΩ ≥ vr · nΩ , sur I × [0, T ]. (3.1)
ψ(v(ψ) − vr ) · nΩ = 0,
La première inégalité signifie que la pression au niveau du sol est inférieure ou égale à la pres-
sion atmosphérique. En particulier, nous supposons qu’en cas d’affleurement, l’eau ruisselante
est évacuée instantanément ou qu’elle n’a pas d’influence sur l’écoulement souterrain. Cette
hypothèse simplificatrice sera levée au chapitre 4. Toutefois nous pouvons d’ores et déjà obser-
ver que dans de nombreuses configurations où la lame d’eau ruisselante a une épaisseur très
faible, l’hypothèse retenue dans ce chapitre s’avère suffisamment précise ; des comparaisons sont
effectuées par exemple par Beaugendre et al. [BEE+ 06]. La seconde inégalité signifie que la vi-
tesse normale est supérieure à l’intensité de la pluie (noter que vr · nΩ ≤ 0) ce qui traduit que
la quantité d’eau maximale que le sol peut infiltrer est la quantité d’eau provenant de la pluie.
La troisième équation précise que la charge hydraulique est nulle lorsque le milieu est saturé et
que la vitesse normale est égale à l’intensité de la pluie lorsque le milieu est insaturé. Il est clair
que cette dernière équation indique le type de condition aux limites à imposer sur le doublet
{I w,t , I d,t } : condition de Dirichlet homogène sur I w,t et condition de Neumann non homogène
sur I d,t . Les deux inégalités sont des conditions complémentaires permettant de déterminer la
décomposition de l’intervalle I.
De nombreuses méthodes existent pour modéliser les drains, parmi lesquelles une méthode
« globale » se distingue des méthodes « locales ». La méthode globale, utilisée par exemple par
Carlier, Kao et Ginzburg [CKG07], consiste à représenter un ensemble de drains par un milieu po-
reux homogène équivalent. Une direction préférentielle est donnée à l’écoulement dans la couche
qui représente le réseau de drains par un tenseur de conductivité anisotrope. Les propriétés hy-
drodynamiques de ce milieu équivalent dépendent des propriétés du sol et des caractéristiques
géométriques du réseau de drains (espacement entre les drains et diamètre des drains). D’autres
méthodes représentent les drains de façon locale en les modélisant soit par des nœuds du maillage,
soit géométriquement (en faisant des trous dans le maillage). Dans les approches nodales, le drain
est assimilé à un terme puits en imposant une charge comme le font Barcelo et Nieber [BN82] ou
une vitesse comme le font Gaharaaty-Sani et King [GSK78]. Dans les sols entièrement saturés,
3.1. Problème de frontière libre 43
une analogie avec un réseau de résistances électriques est utilisée par Vimoke et Taylor [VT62].
Dans la modélisation géométrique, Gureghian et Youngs [GY75] déterminent la vitesse à imposer
sur le bord du drain en supposant que celle-ci varie logarithmiquement avec la distance radiale
au centre du drain. Pour une comparaison de ces quatre méthodes locales, voir les travaux de
Fipps et Skaggs [FS86]. Dans notre cas, nous choisissons une approche plus détaillée donc plus
coûteuse numériquement mais beaucoup plus précise, en modélisant géométriquement les drains
et en utilisant des conditions aux limites du même type que celles utilisées pour l’affleurement
et qui s’écrivent :
ψ ≤ 0,
v(ψ) · nΩ ≥ 0, sur D × [0, T ]. (3.2)
ψv(ψ) · nΩ = 0,
Ces équations sont les équations (3.1) obtenues avec vr = 0. La première inégalité reste in-
changée : la charge hydraulique est négative ou nulle. En revanche la seconde inégalité est
modifiée : la vitesse doit rester positive ou nulle. Nous rappelons qu’avec les conventions choi-
sies, une vitesse normale négative traduit une infiltration dans le milieu poreux, ce que nous
excluons pour un drain. De même la relation de complémentarité est différente : soit la charge
hydraulique est nulle si le milieu est saturé, soit la vitesse normale est nulle si le milieu est insa-
turé. Ainsi le type de condition aux limites à imposer sur le doublet {Dw,t , Dd,t } est le suivant :
condition de Dirichlet homogène sur Dw,t et condition de Neumann homogène sur Dd,t . Notons
qu’ici aussi, nous supposons que les drains évacuent instantanément l’eau.
Nous définissons l’ensemble E comme la réunion de l’interface et des drains
def def def
E = I ∪D ⇒ E w,t = I w,t ∪ Dw,t et E d,t = I d,t ∪ Dd,t ,
Avec les conditions aux limites (3.1) et (3.2) ainsi que les notations précédentes, nous obtenons
le système régissant les écoulements en milieux poreux décrits par l’équation de Richards avec
la prise en compte de l’affleurement et de drains,
∂t [θ(ψ)] + ∇ · v(ψ) = 0 dans Ω × [0, T ],
v(ψ) = −K(ψ)∇(ψ + z) dans Ω × [0, T ],
v(ψ) · nΩ = vN sur (W ∪ B) × [0, T ],
v(ψ) · nΩ = ωv sur {(x, t), x ∈ E d,t }, (3.3)
w,t },
ψ = 0 sur {(x, t), x ∈ E
ψ(., 0) = ψ 0 dans Ω,
(ψ, v(ψ) · nΩ ) ∈ A1 dans E × [0, T ],
où les données sont la condition initiale ψ 0 et la vitesse normale vN . Il s’agit d’un problème
à frontière libre dans lequel les couples (ψ, v(ψ) · nΩ ) pour (x, t) ∈ E × [0, T ] doivent rester
dans un ensemble d’états admissibles du point de vue physique. Cet ensemble admissible A1 ,
représenté à la figure 3.2, est composé de deux branches dans le plan (ψ, vn ) où vn = v(ψ) · nΩ .
La branche horizontale {vn = ωv } est associée à un sol insaturé sur lequel la charge hydraulique
est inférieure ou égale à zéro. La branche verticale {ψ = 0} est associée à un sol saturé où la
vitesse normale est supérieure ou égale à ωv . L’ensemble A1 est donc défini par
def
A1 = {(ψ, vn ) ∈ R2 , ψ ≤ 0, vn ≥ ωv , ψ(vn − ωv ) = 0}. (3.4)
44 Chapitre 3. Affleurement et drains
vn
sol mouillé
sol sec
× ψ
(0, ωv )
3.1.2 Algorithme
Nous donnons dans cette section l’algorithme proposé pour traiter les conditions aux limites
def
de Signorini. L’ensemble Eh est l’ensemble des faces appartenant à E (Eh = {F ∈ Fh, F ⊂ E})
et les deux ensembles Ehd,n et Ehw,n sont définis par
def def
Ehd,n = {F ∈ Fh, F ⊂ E d,nδt } et Ehw,n = {F ∈ Fh, F ⊂ E w,nδt }.
Nous pouvons noter que Ehd,n ∩ Ehw,n = ∅ par construction. Les sous-ensembles Ehd,n et Ehw,n
sont déterminés au cours du calcul et notre algorithme impose que toute arête de l’ensemble Eh
appartient à un de ces sous-ensembles de sorte que
Autrement dit, l’interface entre une zone sèche et une zone mouillée correspond à une interface
surfacique du maillage (un nœud en dimension deux). Pour plus de concision dans la présentation
de l’algorithme, nous introduisons les deux notations suivantes :
– ψhn ← Richards(Ehd,n , Ehw,n , ψhn−1 , ψhn−2 ) désigne la résolution du problème (3.3) sur un pas
de temps en utilisant la méthode SIPG en espace, le schéma BDF1 ou BDF2 en temps et
un partitionnement {Ehd,n , Ehw,n } donné (ψhn−2 n’est pas nécessaire en donnée d’entrée pour
la BDF1),
– vh⋆,n ← Vitesse Normale(Ehd,n , Ehw,n , ψhn ) désigne la détermination de la vitesse normale
vh⋆,n sur l’ensemble E ; celle-ci est définie par
(
def ωvn |F si F ⊂ Ehd,n ,
vh⋆,n |F =
v(ψhn )|F · nΩ + ηKs d−1 n
F ψ h |F si F ⊂ Ehw,n .
Dans l’algorithme 3 nous commençons par supposer que la frontière E est complètement saturée
puisque nous prenons comme premier partitionnement Ehw,n,1 = E et Ehd,n,1 = ∅. Les égalités
3.1. Problème de frontière libre 45
ψhn,1 = 0 et vh⋆,n,1 = ωvn correspondent à ce partionnement initial. Cela implique qu’une condition
de Dirichlet est d’abord imposée sur l’ensemble de la frontière E. L’équation de Richards est alors
résolue et la nouvelle vitesse normale vh⋆,n,p est déterminée. Nous vérifions ensuite que la condition
sur la vitesse vh⋆,n,p ≥ ωvn |F est satisfaite sur toutes les faces F ∈ Eh. Si c’est le cas, la solution et
le partionnement sont acceptés. Dans le cas contraire, une condition de Neumann est imposée
sur les arêtes ne vérifiant pas cette condition et l’équation de Richards est à nouveau résolue.
L’algorithme converge puisque l’ensemble Ehd,n,p augmente (au sens de l’inclusion) avec l’indice
d’itération p alors que l’ensemble Ehd,n,p décroı̂t. Les principaux avantages de cet algorithme sont
la robustesse et une facile adaptation à la dimension 3. De plus, contrairement à un algorithme
de suivi de front (voir [BEE+ 06] par exemple), l’algorithme 3 n’utilise aucune information issue
du pas de temps précédent pour déterminer le partionnement de l’ensemble E.
L’algorithme que nous proposons pour traiter les frontières libres garantit sur l’ensemble des
pas de temps que la charge hydraulique est nulle sur Ehw,n ,
Par définition de ωv , l’égalité (3.6) traduit que la vitesse normale est nulle sur la partie sèche
du drain Dhd,n et égale à la vitesse normale de la pluie sur la partie sèche de l’interface Ihd,n .
L’algorithme 3 assure également que la vitesse normale est supérieure ou égale à la fonction ωvn
sur la partie Ehw,n . La seule condition de complémentarité non vérifiée de manière exacte mais
uniquement de manière approchée concerne la charge hydraulique qui peut éventuellement être
positive sur la partie sèche Ehw,n .
46 Chapitre 3. Affleurement et drains
et ǫn est l’erreur dûe à la résolution du système non linéaire. Nous avons |ǫn | ≤ Cǫalg1 où ǫalg1
est la tolérance choisie dans l’algorithme 1 et C est une constante venant du fait que le critère
de convergence de l’algorithme 1 porte sur la charge hydraulique et non sur la teneur en eau.
La relation (3.7) implique facilement le résultat suivant de conservation du volume d’eau lorsque
la BDF1 est utilisée.
Propriété 1. Soit δV n le défaut de volume d’eau total pendant l’intervalle [(n − 1)δt, nδt] défini
comme
def n−1
δV n = Vgrnd
n
− Vgrnd − (FIn + FDn )δt.
Soit △V n le défaut de volume d’eau total pendant l’intervalle de temps [0, nδt] défini comme
n
X
n def
△V = δV i .
i=1
Alors
|△V n | ≤ nCǫalg1 ,
où C est une constante et ǫalg1 le seuil de tolérance utilisé dans l’algorithme 1.
et [(n − 2)δt, (n − 1)δt], il doit en être de même pour le membre de droite. Nous cherchons donc
à construire des flux de masse F̃In et F̃Dn sur l’intervalle [(n − 1)δt, nδt] de sorte que
Lors du premier pas de temps, le schéma DIRK3 ou de Crank–Nicolson permet d’obtenir F̃I1
et F̃D1 pour commencer la récurrence. Le résultat concernant la conservation du volume d’eau
lorsque la BDF2 est utilisée est le suivant.
Propriété 2. Soit δV n le défaut de volume d’eau total pendant l’intervalle [(n − 1)δt, nδt] défini
comme
def n−1
δV n = Vgrnd
n
− Vgrnd − F̃ n δt, (3.12)
def
avec F̃ n = F̃In + F̃Dn . Soit △V n le défaut de volume d’eau total pendant l’intervalle de temps
[0, nδt] défini comme précédemment. Alors
1
|△V n | ≤ |δV 1 | + nCǫalg1 ,
2
où C est une constante et ǫalg1 le seuil de tolérance utilisé dans l’algorithme 1.
d’où
X n
1 1
|△V | ≤ |δV 1 | +
n
|ǫi | ≤ |δV 1 | + nCǫalg1 .
2 2
i=1
Des raisonnements similaires avec des formules de récurrence plus longues peuvent être
considérés pour des BDF d’ordre quelconque.
48 Chapitre 3. Affleurement et drains
3.2 Résultats
Dans les deux cas tests présentés, l’un avec affleurement et l’autre avec un drain, nous
effectuons un bilan de masse afin de vérifier que celle-ci est bien conservée. Pour déterminer si
les conditions de complémentarité sur l’ensemble E sont vérifiées, nous représentons les couples
{ψ, vn } sur l’ensemble admissible A1 à différents instants de la simulation. L’évolution du type
de condition aux limites imposé sur l’ensemble Eh est suivie en indiquant, en fonction du temps,
le nombre de faces dans Eh ayant une condition de Dirichlet. Pour repérer l’état de stationnarité
du système, nous introduisons l’indicateur Is défini par
3.2.1 Affleurement
Nous reprenons un cas test aboutissant à l’affleurement d’une nappe proposé par Beaugendre
et al. [BEE+ 06] et inspiré d’un article d’Abdul et Gilham [AG84]. Une caractéristique du domaine
est la pente de la surface qui vaut 14% (voir la figure 3.3). La géométrie, le temps final de la
simulation et le pas de temps utilisés sont
La condition initiale est hydrostatique avec une nappe située à 0.7m et des conditions aux limites
d’imperméabilité sont imposées sur les parois et le fond,
ψ 0 = −z + 0.7m dans Ω,
vN = 0 sur (W ∪ B) × [0, T ].
Une pluie d’intensité constante égale à 10% de la conductivité à saturation est imposée sur le
sol pendant toute la simulation,
1m
0.8m
vN = 0
•
(0, 0) 1.4m
Le sol est constitué d’un sable dont les propriétés hydrodynamiques sont données par les relations
de Van-Genuchten (2.30) avec les valeurs suivantes pour les différents paramètres,
4
0.52
0.5 3
0.48
2
0.46
1
0.44
0.42 0
−40 −35 −30 −25 −20 −15 −10 −5 0 −40 −35 −30 −25 −20 −15 −10 −5 0
Charge hydraulique ψ [cm] Charge hydraulique ψ [cm]
Une partie du domaine étant saturée à l’instant initial, nous utilisons le schéma de Crank–
Nicolson lors du premier pas de temps. De plus, nous choisissons un maillage de 2006 mailles qui
est plus fin dans la partie supérieure pour mieux capter l’affleurement de la nappe : l’ensemble
I comporte 78 arêtes alors que le fond B n’en compte que 17 (voir la figure 3.5).
Sur la figure 3.6, la position de la nappe est dessinée en trait continu à plusieurs instants
significatifs au cours de la simulation. La partie saturée du sous-sol (qui se trouve en-dessous de la
nappe) est de couleur foncée alors que la partie insaturée est de couleur claire. A l’initialisation,
la répartition de la pression est hydrostatique avec une nappe horizontale située à une altitude
de 0.7m. La pluie imposée sur la partie supérieure du domaine entraine une augmentation
progressive de la pression ce qui se traduit par une montée de la nappe notamment sur la partie
droite du domaine où le milieu est le moins insaturé en surface (1h). Après 1h30 de pluie,
la quantité d’eau infiltrée est suffisante pour saturer la partie droite du domaine et provoquer
l’affleurement de la nappe. La dernière position correspond à l’état d’équilibre entre l’infiltration
et l’exfiltration dans le système. La zone de saturation du milieu est alors maximale.
0h 1h
100
0
0 100 140 0 100 140
1h30 6h
100 100
0 0
0 100 140 0 100 140
Fig. 3.6: Position de la nappe à différents instants.
3.2. Résultats 51
La figure 3.7 montre les différents couples charge hydraulique/vitesse normale sur A1 pendant
la simulation. Pour différents instants et pour chaque arête de l’ensemble I, le couple (ψhn , vhn ·nΩ )
est représenté par une croix (la valeur moyenne de ψhn est utilisée sur chaque maille). Les phases
précédemment décrites sont illustrées par les positions successives des nuages de points :
- A 10min, la charge hydraulique est négative et la vitesse normale est égale à l’intensité de
la pluie. Le nuage de points est seulement sur la branche {vn = vr ·nΩ } traduisant une infiltration
complète de la pluie dans le sol.
- A 1h, les couples sont toujours sur la branche {vn = vr · nΩ } mais avec une augmentation
de la charge hydraulique puisque le milieu se sature.
- A 1h40, la vitesse normale est supérieure à l’intensité de la pluie pour quelques faces.
Le nuage de points est situé sur les deux branches et la nappe commence à affleurer.
- A 6h, le nuage se trouve sur les deux branches avec une augmentation de la vitesse normale
traduisant une saturation du milieu. Quelques points s’éloignent quelque peu de l’ensemble
admissible A1 illustrant le fait que la condition de non-positivité de la charge hydraulique sur
la partie sèche de la frontière libre n’est pas assurée de manière exacte par notre algorithme.
10min 1h
0.0005 0.0005
0.0004 0.0004
0.0003 0.0003
0.0002 0.0002
0.0001 0.0001
0 0
-0.0001 -0.0001
-0.0002 -0.0002
-30 -25 -20 -15 -10 -5 0 5 -30 -25 -20 -15 -10 -5 0 5
1h40 6h
0.0005 0.0005
0.0004 0.0004
0.0003 0.0003
0.0002 0.0002
0.0001 0.0001
0 0
-0.0001 -0.0001
-0.0002 -0.0002
-30 -25 -20 -15 -10 -5 0 5 -30 -25 -20 -15 -10 -5 0 5
Dans ce cas test d’affleurement, l’ensemble E est réduit à l’interface I puisqu’il n’y a pas de
drain. Le flux de masse F̃In à l’interface à l’instant discret nδt peut être décomposé en un flux
n
d’infiltration F̃I,in n
et un flux d’exfiltration F̃I,out ,
Z Z
def n def def
F̃In = n
F̃I,in + n
F̃I,out , avec F̃I,in = − vh⋆,n,− et n
F̃I,out = − vh⋆,n,+ ,
I I
où vh⋆,n,− < 0 et vh⋆,n,+ > 0 désignent respectivement la partie négative et la partie positive de
vh⋆,n sur I. Pour obtenir le bilan de masse sur toute la durée de la simulation, nous utilisons la
décomposition du flux sur l’interface, multiplions l’équation (3.12) par la masse volumique de
l’eau ρ et sommons sur l’intervalle [0, T ],
n
X n
X n
X
n 0 i i
ρ(Vgrnd − Vgrnd ) = ρδtF̃I,in + ρδtF̃I,out + ρδV i .
| {z } i=1 | {z } i=1 | {z } i=1 | {z }
∆Mn
grnd Miin Miout Ei
La grandeur ∆Mngrnd est la variation totale de masse pendant l’intervalle [0, nδt]. Les deux
quantités Miin et Miout représentent respectivement la masse d’eau qui entre et qui sort du système
entre (i − 1)δt et iδt. L’erreur commise
P P entre (i − 1)δt et iδt est
sur la conservation de la masse
notée Ei . Les quantités ∆Mngrnd , ni=1 Miin (notée ΣMnin en abrégé), ni=1 Miout (notée ΣMnout en
P
abrégé) et ni=1 Ei (notée ΣEn en abrégé) sont représentées sur la figure 3.8. Le bilan de masse
confirme que la pluie s’infiltre totalement pendant les deux premières heures de la simulation
car la variation de masse est égale à la masse apportée par la pluie. L’eau commence à s’exfiltrer
pendant la phase suivante (entre 2 et 4 heures). Cela correspond à la période où l’affleurement
de la nappe se produit. Après 4h de simulation, la masse d’eau dans le sol devient constante et la
masse d’eau exfiltrée est égale à la masse d’eau infiltrée. Ce bilan montre également que l’erreur
commise sur la conservation de la masse est faible par rapport à la masse d’eau initialement
contenue dans le sol (écart relatif de l’ordre de 0.02%).
12
ΣMn
in
8
∆Mn
grnd
Masse [kg]
0
ΣEn
−4 ΣMn
out
0 1 2 3 4 5 6
Temps [h]
Le nombre d’arêtes ayant une condition de Dirichlet est indiqué sur la figure 3.9 et confirme
les trois phases distinguées dans le bilan de masse. La simulation commence avec une condition
de Neumann sur l’ensemble de l’interface I. Ensuite l’affleurement de la nappe débute avec
3.2. Résultats 53
l’apparition de la première arête avec une condition de Dirichlet. Le nombre d’arêtes de Dirichlet
se stabilise lorsque le système atteint l’équilibre.
L’indicateur Is reporté sur la figure 3.10 traduit correctement l’état de stationnarité du
système en devenant très faible (inférieur à 10−8 ) quand la nappe atteint sa position finale.
−2
50 10
Nombre d’arêtes de Dirichlet
40
−4
10
Indicateur Is
30
−6
10
20
10 −8
10
0
−10
10
0 1 2 3 4 5 6 0 1 2 3 4 5 6
Temps [h] Temps [h]
Le champ de vitesse (en mm.h−1 ) à proximité de la surface du domaine est tracé sur la fi-
gure 3.11 à plusieurs instants. Nous avons choisi de l’évaluer en utilisant le gradient par morceaux
sur chaque élément (noté ∇h) de la charge hydraulique (sans utiliser la correction H(div, Ω)-
conforme introduite par Ern, Nicaise et Vohralı́k [ENV07]),
où K(ψh) est la fonction constante par morceaux égale à la valeur moyenne de K(ψh) sur chaque
élément (évaluée par une quadrature à 7 points de degré 5). Comme la charge hydraulique est P1 ,
les deux composantes de la vitesse sont constantes sur chaque triangle. Dans le champ de vitesse
en haut de la figure 3.11, les triangles proches de la surface ont des vecteurs vitesse verticaux et
dirigés vers le bas ce qui traduit l’infiltration totale de la pluie dans le sol. Les autres éléments
du domaine ont un vecteur vitesse d’intensité faible. Le deuxième champ de vitesse représenté à
1 heure présente également des vecteurs vitesse verticaux descendants sur les triangles proches
de la surface. La montée de la nappe est visible sur les éléments de la partie gauche du domaine
proches de la cote 0.7m puisque les vecteurs vitesse de ces éléments sont verticaux ascendants.
Le troisième champ est représenté à 1h40 et se différencie du champ précédent par l’exfiltration
qui est visible sur la dernière arête de l’interface I. L’exfiltration apparait clairement dans le
dernier champ sur les triangles où les vecteurs vitesse sont dirigés vers l’extérieur du domaine.
Il est clair que les lignes de courant proches de la surface ont leurs deux extrémités sur l’interface
I confirmant ainsi que l’eau infiltrée sur une partie de la surface s’exfiltre par une autre partie.
Une vitesse normale quasiment nulle sur quelques éléments est visible puisque certains vecteurs
vitesses sont pratiquement parallèles à l’interface.
54 Chapitre 3. Affleurement et drains
10min
100
95
90
85
80
75
70
0 20 40 60 80 100 120 140
1h
100
95
90
85
80
75
70
0 20 40 60 80 100 120 140
1h40
100
95
90
85
80
75
70
0 20 40 60 80 100 120 140
6h
100
95
90
85
80
75
70
0 20 40 60 80 100 120 140
3.2.2 Drain
Nous proposons un cas test avec un drain situé à l’intérieur d’un domaine rectangulaire (voir
la figure 3.12). La géométrie, le temps final de la simulation et le pas de temps utilisés sont
La condition initiale est hydrostatique avec une nappe coı̈ncidant avec le fond du domaine et
des conditions aux limites d’imperméabilité sont imposées sur les parois et le fond,
ψ 0 = −z dans Ω,
vN = 0 sur (W ∪ B) × [0, T ].
Une pluie d’intensité constante égale à 10% de la conductivité hydraulique à saturation est
imposée sur le sol pendant les 6 premières heures uniquement,
Le centre cd du drain de rayon rd est situé proche du fond sur l’axe de symétrie du domaine
parallèle aux parois latérales,
1m
2cm
vN = 0 ▽
(0, 0) •
2m 2m
Le sol est constitué du sable utilisé dans le cas test d’Haverkamp qui est caractérisé par les
relations (2.27) et dont nous rappelons les valeurs des trois paramètres physiques θr , θs et Ks
La teneur en eau et la conductivité hydraulique sont représentées sur la figure 2.5 de la sec-
tion 2.2.3. Les paramètres de ce cas test (géométrie et intensité de la pluie) ne conduisent pas à
l’affleurement de la nappe. L’ensemble E sur lequel s’appliquent des conditions de Signorini est
ainsi réduit au seul drain D.
56 Chapitre 3. Affleurement et drains
Concernant la prise en compte du drainage, il est important de noter que les dimensions
des drains sont de l’ordre du centimètre alors que les dimensions du domaine sont ici de l’ordre
du mètre. Le rayon du drain étant de 1cm, son périmètre est d’environ 6cm et, en utilisant
une vingtaine de segments pour représenter convenablement sa circonférence, la taille des arêtes
est de l’ordre de 0.3cm. Un maillage structuré du domaine, constitué par exemple de triangles
rectangles ayant leur plus petit côté égal à 0.3cm aurait 720000 éléments (300 × 1200 × 2).
Pour réduire ce nombre nous effectuons un raffinement local autour du drain comme le montrent
les figures 3.13 et 3.14. En prenant des triangles de 10cm de côté sur les bords extérieurs du
domaine, le maillage localement raffiné contient 2009 éléments. La présence du drain multiplie
ainsi par 2.5 le nombre d’éléments du maillage structuré constitué de triangles rectangles de
10cm de côté puisque ce maillage contient 800 éléments (10 × 40 × 2). Le gain par rapport au
maillage fin uniforme décrit précédemment est considérable car le nombre d’éléments est divisé
par 360. La simulation a été faite sur l’ensemble du domaine mais aurait pu être réalisée sur une
seule moitié pour d’évidentes raisons de symétrie.
Sur la figure 3.15, la position de la nappe est repérée en trait continu à plusieurs instants.
Comme dans le cas test précédent, la partie saturée du sous-sol est de couleur foncée alors que la
partie insaturée est de couleur claire. La nappe est horizontale pendant la saturation de la partie
du sous-sol située sous le drain. Ensuite l’effet du drain est clairement visible sur la position de
la nappe à 3 et 6 heures. En effet, en l’absence de drain, les iso-valeurs de la charge hydraulique
seraient horizontales ce qui correspondrait à une répartition de pression hydrostatique. Après
l’arrêt de la pluie il y a un abaissement progressif de la nappe. Cela est observé pour la position
finale de la nappe (18 heures) qui est presque horizontale au niveau du drain. Comme attendu,
toutes les positions de la nappe présentent une symétrie par rapport à l’axe vertical passant par
le centre du drain.
2h
100
0
3h
100
0
6h
100
0
18h
100
0
Fig. 3.15: Position de la nappe à différents instants.
58 Chapitre 3. Affleurement et drains
La figure 3.16 montre les différents couples charge hydraulique/vitesse normale sur l’ensemble
admissible A1 pendant la simulation. Les phases précédemment décrites sont illustrées par les
positions successives des nuages de points :
- A 2h, la charge hydraulique est négative et la vitesse normale est nulle. Le nuage de points
est seulement sur la branche {vn = 0} traduisant une saturation du sol autour du drain.
- A 3h, les couples sont sur la branche {ψ = 0} puisque le sol est complètement saturé au
niveau du drain.
- A 6h, les couples se sont déplacés vers le haut de la branche verticale traduisant une
augmentation de la vitesse normale.
- A 18h, le nuage se trouve sur les deux branches puisque le drain a permis de désaturer le
sol après l’arrêt de la pluie.
2h 3h
0.07 0.07
0.06 0.06
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-5 -4 -3 -2 -1 0 1 -5 -4 -3 -2 -1 0 1
6h 18h
0.07 0.07
0.06 0.06
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-5 -4 -3 -2 -1 0 1 -5 -4 -3 -2 -1 0 1
Le bilan de masse est obtenu comme précédemment en multipliant l’équation (3.12) par la
masse volumique de l’eau ρ et en effectuant la somme sur l’intervalle [0, nδt],
n
X Xn Xn
n 0
ρ(Vgrnd − Vgrnd )= ρδtF̃Ii + ρδtF̃Di + ρδV i .
| {z } i=1 | {z } i=1 | {z } i=1 | {z }
∆Mn
grnd Miin Miout Ei
800
ΣMn
in
400
Masse [kg]
∆Mn
grnd
0
ΣEn
−400
−800 ΣMn
out
0 2 4 6 8 10 12 14 16 18
Temps [h]
Le nombre d’arêtes ayant une condition de Dirichlet est indiqué sur la figure 3.18 et corrobore
les trois phases identifiées dans le bilan de masse. Le milieu étant insaturé à l’initialisation, la
simulation débute avec une condition de Neumann homogène imposée sur toutes les arêtes du
drain. Puis, la zone au niveau du drain devient saturée avec le passage à des conditions aux
limites de Dirichlet sur toutes les arêtes situées au bord du drain. Quand la pluie cesse, le sol se
draine lentement et les conditions aux limites imposées sur certaines arêtes du drain redeviennent
peu à peu (après 10h de simulation) de type Neumann.
L’indicateur de stationnarité Is tracé sur la figure 3.19 montre que le système tend exponen-
tiellement vers une solution stationnaire sur l’intervalle [3h, 6h]. Cette première solution station-
naire correspond à un état d’équilibre entre la pluie et le drainage. Dans cette phase, l’indicateur
est divisé par 10 en 2.4 heures environ : le temps caractéristique de la convergence exponentielle
vers l’état stationnaire est de l’ordre d’une heure. Après l’arrêt de la pluie, le système tend vers
60 Chapitre 3. Affleurement et drains
une seconde solution stationnaire qui est l’état d’équilibre hydrostatique ayant une position de
nappe située au niveau du drain. Dans cette deuxième phase, les variations de l’indicateur de
stationnarité coı̈ncident avec les changements de conditions aux limites autour du drain.
−2
10
20
Nombre d’arêtes de Dirichlet
−4
10
15
Indicateur Is
−6
10 10
5 −8
10
0
−10
10
0 2 4 6 8 10 12 14 16 18 0 2 4 6 8 10 12 14 16 18
Temps [h] Temps [h]
La figure 3.20 présente à plusieurs instants le champ de vitesse constant par morceaux (en cm
par 30s) sur quelques éléments près du drain. Le premier champ confirme que les vitesses sont
pratiquement nulles au début de la simulation. Le champ tracé à 3h est radial et converge vers le
centre du drain car toutes les conditions sont de Dirichlet. Les vecteurs vitesse sur les triangles
proches du fond sont horizontaux, confirmant que la vitesse normale est nulle sur ce bord où
une condition de Neumann homogène est imposée. Il en est de même du champ tracé à 6h dans
lequel la norme des vecteurs vitesse est plus importante qu’à 3h car le milieu est plus saturé.
Le champ de vitesse à 12h est différent des deux précédents. La pluie ayant cessé, les vecteurs
sont d’amplitude moindre. De plus, les deux types de conditions aux limites sont imposés sur
le drain : Neumann sur les 6 arêtes supérieures et Dirichlet sur les 14 arêtes inférieures où
l’exfiltration se produit.
2h 3h
6 6
5 5
4 4
3 3
2 2
1 1
0 0
197 198 199 200 201 202 203 197 198 199 200 201 202 203
3.2. Résultats 61
6h 12h
6 6
5 5
4 4
3 3
2 2
1 1
0 0
197 198 199 200 201 202 203 197 198 199 200 201 202 203
Les conditions aux limites de Signorini considérées dans le chapitre précédent ne prennent
pas en compte le ruissellement surfacique lorsqu’une nappe affleure ou en cas de ruissellement
hortonien lorsque les capacités d’infiltration du sol sont insuffisantes. Nous proposons dans ce
chapitre d’inclure ce phénomène dans un modèle de couplage entre l’équation de Richards et
l’approximation de l’onde cinématique.
Nous commençons par présenter le problème continu avec notamment la condition de cou-
plage imposée sur la charge hydraulique et la hauteur d’eau à la surface du sol. L’équation
de l’onde cinématique est ensuite discrétisée par une méthode de volumes finis et le schéma
d’Euler explicite. Nous proposons deux algorithmes de couplage pour lesquels une propriété de
conservation de la masse est établie. Nous détaillons ensuite les quatre cas tests considérés pour
évaluer nos algorithmes dans diverses situations : assèchement d’un sol, ruissellement dû à la
pluie, ruissellement dû à l’injection d’eau dans le sol et ruissellement hortonien. Enfin, nous
présentons une application concrète en hydrologie relative au fonctionnement d’un petit bassin
versant drainé.
Le système régissant les écoulements superficiels est composé de l’équation de l’onde cinématique
avec un terme source −vr · nΩ provenant de la pluie et un terme source (ou puits) v(ψ) · nΩ
provenant du couplage, ainsi qu’une condition initiale h0 et une condition à la limite hA ,
∂t h + ∂x ϕ(h, S) = (v(ψ) − vr ) · nΩ sur I × [0, T ],
h(·, 0) = h0 sur I, (4.2)
h(A, ·) = hA en A × [0, T ].
64 Chapitre 4. Couplage avec le ruissellement surfacique
A•
h
I •B
I d,t
I w,t
W Ω W
nΩ
L’équation de l’onde cinématique est une loi de conservation scalaire hyperbolique si la hauteur
d’eau est strictement positive. De plus, une condition à la limite est uniquement nécessaire au
point A puisque les ondes se propagent de la gauche vers la droite.
Le problème modèle pour le couplage des écoulements souterrains et surfaciques consiste à
trouver le couple de fonctions (ψ, h) tels que
ψ est solution de (4.1) dans Ω × [0, T ],
h est solution de (4.2) sur I × [0, T ], (4.3)
(ψ, h) ∈ A2 sur I × [0, T ],
dans lequel A2 est l’ensemble des états (ψ, h) admissibles. L’ensemble admissible A2 , tracé à la
figure 4.2, comprend deux branches. La branche {h = 0} est associée à un sol insaturé (ou sec
par abus de langage) sur lequel la charge hydraulique est inférieure ou égale à zèro. La branche
{h = ψ} est associée à un sol saturé (ou mouillé par abus de langage) où la charge hydraulique
et la hauteur d’eau sont égales. Ainsi, l’ensemble A2 est défini par
def
A2 = {(ψ, h) ∈ R2 , h = ψ + }, (4.4)
sol mouillé
sol sec
× ψ
(0, 0)
Nous considérons principalement des situations dans lesquelles l’écoulement de surface est
produit par l’exfiltration. Des déplacements du couple (ψ, h) sur l’ensemble admissible sont donc
attendus. Des situations plus drastiques, comme des ondes ruisselantes sur des sols insaturés dans
le cas du ruissellement hortonien, peuvent conduire dans certains cas à un départ de l’ensemble
admissible en particulier si le sol est très sec. Dans ces cas limites, d’autres modèles de couplage
peuvent être plus pertinents. Si les processus d’infiltration sont très lents par exemple, un modèle
où les écoulements de surface se produisent sur un sol insaturé peut être envisagé. Ces études
sortent du cadre de ce travail et ne seront donc pas abordées.
Les deux ensembles I w,t et I d,t dépendent du temps et sont définis en fonction de la quantité
h de la façon suivante,
def def
I w,t = {x ∈ I, h(x, t) > 0} et I d,t = {x ∈ I, h(x, t) = 0}.
Il sera utile de reformuler le système (4.1) en désignant par ωψ et ωv les données à imposer
respectivement dans les conditions aux limites de Dirichlet et de Neumann sur I × [0, T ],
∂t [θ(ψ)] + ∇ · v(ψ) = 0 dans Ω × [0, T ],
v(ψ) = −K(ψ)∇(ψ + z) dans Ω × [0, T ],
v(ψ) · nΩ = vN sur (W ∪ B) × [0, T ],
(4.5)
v(ψ) · n Ω = ω v sur {(x, t), x ∈ I d,t },
ψ = ωψ sur {(x, t), x ∈ I w,t },
ψ(·, 0) = ψ 0 sur Ω.
t t + n′ δt′
Ecoulement superficiel δt′
| | | | |
δt temps
Ecoulement souterrain
t t + δt
Fig. 4.3: Multi-pas de temps pour les écoulements souterrain et superficiel.
66 Chapitre 4. Couplage avec le ruissellement surfacique
Soient xi , li , xi− 1 et xi+ 1 définis respectivement comme le centre, la longueur, et les sommets
2 2
gauche et droit d’une arête quelconque ei du maillage de l’interface I (voir figure 4.4). La pente
de ei est notée Si . Comme la fonction de flux ϕ est convexe et que la hauteur d’eau est positive
ou nulle, le flux de Godunov coı̈ncide avec le flux amont, ce qui conduit au schéma
∀k ∈ {0 · · · n′ − 1}, ∀i ∈ {1 · · · NI }, (4.6)
′ Z
δt
hin−1,k = hin−1,k−1 + n−1,k−1
ϕ(hi−1 , Si−1 ) − ϕ(hin−1,k−1 , Si ) − li vrn−1,k−1 · nΩ + vh⋆,n ,
li ei
h−1
hi−1
hi
x− 1
•
2 xi− 3 hi+1
•
2 xi− 1
A •
2 xi xi+ 1
× •
2 xi+ 3
2
•
li
Le schéma (4.6) nécessite la connaissance de la hauteur d’eau à t = 0 et sur une arête fictive
à gauche de la première arête pour tous les temps discrets,
′
∀i ∈ {1 · · · NI }, h0i = h0 (xi ) et ∀n ∈ {1 · · · NT }, ∀k ∈ {0 · · · n′ − 1}, hn,k nδt+kδt
−1 = hA .
En utilisant la définition de la fonction flux ϕ et en supposant que hmax est une borne supérieure
a priori de la hauteur d’eau h sur I × [0, T ], cette condition est assurée si
3 −1
δt′ ≤ 2/3
· min li Si 2 .
5Khmax 1≤i≤N I
Sans la pluie et le terme de couplage, la condition CFL garantit un principe du maximum discret.
finis pour l’onde cinématique décrit au paragraphe précédent, nous obtenons ainsi un schéma
global pour approcher le système couplé (4.3) à condition que l’évolution dans le temps des
variables de couplage {Ihd,n , Ihw,n , ωvn , ωψn } pour n ∈ {1 · · · NT } soit connue. L’évolution dans
le temps du système couplé est conçu pour satisfaire deux objectifs. Le premier est d’assurer
que l’approximation de (ψ, h) reste dans l’ensemble admissible A2 pour tous les pas de temps.
Le second est de garantir une conservation de la masse du système total (écoulements souterrain
et surfacique). L’algorithme proposé s’appelle « algorithme de couplage à un pas » puisqu’il utilise
la BDF1 dans laquelle seule la solution au pas de temps précédent est nécessaire. Pour simplifier
la présentation de cet algorithme, nous définissons
– ψhn ← Richards BDF1(Ihd,n , Ihw,n , ωvn , ωψn , ψhn−1 ) comme la résolution par l’algorithme 1 de
l’équation de Richards sur un pas de temps par la méthode SIPG, la BDF1 et les conditions
aux limites sur I données par {Ihd,n , Ihw,n , ωvn , ωψn },
– vh⋆,n ← Vitesse Normale(Ihd,n , Ihw,n , ωvn , ωψn , ψhn ) comme l’évaluation de la vitesse normale
vh⋆,n sur I, donnée comme précédemment par
(
def ωvn |F si F ∈ Ihd,n ,
vh⋆,n |F = w,n
v(ψhn )|F · nΩ + ηKs d−1 n n
F (ψh |F − ωψ |F ) si F ∈ Ih .
Nous détaillons maintenant le principe de l’algorithme 4. Tout d’abord, une prédiction h̃nh
de la hauteur d’eau est calculée sans terme de couplage (vh⋆,n = 0). Cette prédiction sert de
condition aux limites de Dirichlet dans l’équation de Richards. Comme le schéma de Godunov
satisfait un principe du maximum discret, h̃ni ≥ 0 pour tout i ∈ {1 · · · NI }, ce qui implique que
Ihd,n,1 = ∅ et Ihw,n,1 = I. Nous commençons ainsi par supposer que l’interface I est complètement
saturée, la détermination de ωvn,1 étant alors inutile. Ensuite, l’équation de Richards est avancée
d’un pas de temps et la vitesse normale vh⋆,n,p correspondante est calculée puis utilisée afin
d’évaluer la hauteur d’eau hn,p n,p
h . Le signe de hh est alors examiné sur toutes les faces de I.
Si la hauteur d’eau est positive ou nulle sur l’ensemble des faces de I, les évaluations de la
charge hydraulique et de la hauteur d’eau sont acceptées comme solutions du système couplé
à l’instant nδt. Dans le cas contraire, une nouvelle partition de I est déterminée en imposant
une condition de Neumann sur les faces où la hauteur d’eau est négative. Cette condition est
évaluée pour assècher complètement l’eau de surface ce qui traduit une infiltration complète dans
le sol. Un nouveau couple charge hydraulique/hauteur d’eau est calculé jusqu’à convergence.
La convergence de ce processus itératif est assurée car l’ensemble Ihd,n,p croı̂t avec p (au sens de
l’inclusion) alors que l’ensemble Ihw,n,p décroı̂t. Un point essentiel est que l’algorithme 4 garantit
une hauteur d’eau positive ou nulle. De plus, sur la partie mouillée de l’interface, nous avons
l’égalité
∀n ∈ {1 · · · NT }, ∀F ∈ Ihw,n , ψhn |F = h̃nh |F ,
puisque la valeur imposée sur la condition de Dirichlet ωψn,p est fixée pendant la boucle en p dans
l’algorithme 4. Ce n’est pas la condition ψ = h imposée par l’ensemble admissible A2 mais une
approximation en O(δt) de celle-ci. Sur la partie sèche de l’interface, la hauteur d’eau est nulle
et l’inégalité suivante est vérifiée,
∀n ∈ {1 · · · NT }, ∀F ∈ Ihd,n , ψhn |F ≤ h̃nh |F .
Il s’agit aussi d’une approximation en O(δt) de la condition ψ ≤ 0 qui est imposée sur l’ensemble
admissible. Enfin, dans le cas particulier où la hauteur d’eau hin−1 et les flux amont sont nuls sur
une arête quelconque ei pendant l’intervalle de temps [(n − 1)δt, nδt], la condition de Neumann
dans l’équation de Richards sur cette arête est égale à l’intensité de la pluie. Nous retrouvons
ainsi la situation traitée au chapitre 3.
Nous étudions maintenant la conservation du volume d’eau. En prenant une fonction test
φ égale à 1 dans la méthode SIPG, en sommant sur l’ensemble des éléments du maillage et en
utilisant la BDF1 pour approcher le terme instationnaire, nous obtenons
n n−1
Vgrnd − Vgrnd = FIn + FWB
n
δt + ǫn , (4.7)
dans laquelle FIn (resp. FWBn ) est le flux pendant l’intervalle [(n − 1)δt, nδt] à travers l’interface I
n def
FABr = FAn + FBn + Frn ,
n′ n ′ n Z
′
def δt′ X n−1,k def δt′ X n−1,k def δt′ X
FAn = ϕ(hA ), FBn = − ϕ(hN ) et Frn = − vrn−1,k · nΩ .
δt δt I δt I
k=1 k=1 k=1
Le volume total d’eau contenu dans le système couplé est la somme du volume d’eau de chaque
sous-système,
def
V n = Vgrnd
n n
+ Vover .
Lorsque les équations (4.7) et (4.9) sont sommées, les flux à l’interface s’annulent, ce qui donne
V n − V n−1 = FWB
n n
+ FABr δt + ǫn .
Propriété 3. Soit δV n le défaut de volume d’eau total pendant l’intervalle [(n − 1)δt, nδt] défini
comme
def
δV n = V n − V n−1 − (FWBn n
+ FABr )δt.
Soit △V n le défaut de volume d’eau total pendant l’intervalle de temps [0, nδt] défini par
n
X
def
△V n = δV i .
i=1
Alors
|△V n | ≤ nCǫalg1 ,
où C est une constante et ǫalg1 le seuil de tolérance utilisé dans l’algorithme 1.
Pour identifier l’expression de ΦnI , nous observons que l’utilisation de la BDF2 change (4.7) en
3 n n−1 1 n−2
Vgrnd − 2Vgrnd + Vgrnd = FIn + FWB
n
δt + ǫn ,
2 2
qui peut se réécrire
3 n n−1
1 n−1 n−2
n
Vgrnd − Vgrnd − Vgrnd − Vgrnd − FWB δt = FIn δt + ǫn , (4.11)
2 2
70 Chapitre 4. Couplage avec le ruissellement surfacique
3 n n−1
1 n−1 n−2
3 n 1 n−1 3
n 1 n−1
V − Vover − Vover − Vover + − FABr + FABr δt = − ΦI + ΦI δt. (4.12)
2 over 2 2 2 2 2
Le nouveau flux à l’interface ΦnI est déterminé pour que le flux de masse FIn soit exactement
équilibré par le flux à l’interface dans (4.11), d’où
3 1 2 1
FIn = ΦnI − ΦIn−1 ⇒ ΦnI = FIn + ΦIn−1 .
2 2 3 3
Au premier pas de temps où un schéma implicite à un pas est utilisé, la conservation du volume
def
d’eau est directement imposée en prenant Φ1I = FI1 .
L’algorithme de couplage à deux pas est présenté par l’algorithme 5 dans lequel seules les
différences avec l’algorithme 4 sont indiquées. La modification principale concerne l’évaluation
de hn,pi . La donnée de Neumann ωv
n,p
est aussi modifiée afin que la condition de Neumann
conduise à un assèchement de la face correspondante. Enfin, l’approximation discrète ψhn−2 au
temps (n − 2)δt est ajoutée aux données d’entrée et la vitesse à l’interface vh⋆,n est ajoutée
aux données de sortie à chaque pas de temps puisque celle-ci est utilisée dans le pas de temps
suivant. Le résultat concernant la conservation du volume d’eau pour l’algorithme à deux pas
est le suivant.
Propriété 4. Soit δV n le défaut de volume d’eau total pendant l’intervalle [(n − 1)δt, nδt] défini
comme
def
δV n = V n − V n−1 − (F̃WB
n n
+ FABr )δt, (4.13)
def n + 1 F̃ n−1 . Soit △V n le défaut de volume d’eau total pendant l’intervalle de
n
où F̃WB = 23 FWB 3 WB
temps [0, nδt] défini comme précédemment. Alors
1
|△V n | ≤ |δV 1 | + nCǫalg1 , (4.14)
2
où C est une constante et ǫalg1 la tolérance dans l’algorithme 1.
4.2. Résultats 71
Démonstration. Grâce à la relation (4.12), les termes de couplage sont éliminés lorsque les
équations (4.10) et (4.11) sont sommées,
3 n 1 n 3 n 1 n−1
V − V n−1 − V n−1 − V n−2 − FWB δt + FABr δt − F δt = ǫn .
2 2 2 2 ABr
En utilisant la définition de δV n , nous obtenons la relation de récurrence
1 2
δV n = δV n−1 + ǫn .
3 3
La fin de la démonstration est identique à celle du chapitre 3.
4.2 Résultats
L’algorithme 5 que nous avons implémenté dans le logiciel R+R DG est évalué sur quatre
cas tests : assèchement d’un sol, ruissellement dû à la pluie, ruissellement dû à l’injection d’eau
dans le sol et ruissellement hortonien. Les propriétés hydrodynamiques du sol sont données par
les relations d’Haverkamp (2.27) avec les valeurs suivantes pour les différents paramètres,
0.5
1
Teneur en eau θ [−]
0.4 0.8
0.6
0.3
0.4
0.2
0.2
0.1 0
−50 −40 −30 −20 −10 0 −80 −70 −60 −50 −40 −30 −20 −10 0
Charge hydraulique ψ [cm] Charge hydraulique ψ [cm]
Fig. 4.5: Propriétés hydrodynamiques du sol dans les cas tests de couplage.
4.2.1 Assèchement
Dans ce premier cas test, l’écoulement surfacique et le drainage du sol sont provoqués par le
débit à l’exutoire situé au point B. L’interface I comprend trois parties (voir la figure 4.6),
dont les pentes respectives sont J1 , J2 et J3 . La géométrie et le temps final de simulation sont
où L est la longueur du domaine et P l’altitude du point le plus bas de I. La condition initiale
est hydrostatique avec une nappe située à la cote z = 30.25cm et des conditions aux limites
d’imperméabilité sont imposées sur les parois et le fond,
ψ 0 = −z + 0.3025m dans Ω,
vN = 0 sur (W ∪ B) × [0, T ].
Pour l’écoulement surfacique, la condition initiale est une surface libre horizontale et la hauteur
d’eau imposée en A est nulle,
h0 = (−z + 0.3025)+ m sur I,
hA = 0 en A × [0, T ].
0.3025m 0.3m
vN = 0
•
(0, 0) 3m
Les calculs sont effectués avec un maillage quasi-uniforme de 2063 triangles (ce qui correspond
à une taille de maille d’environ 3.5cm) et des pas de temps δt = 2.5s et δt′ = 0.25s. Nous avons
vérifié que les vitesses normales à l’interface obtenues lorsque δt = δt′ = 0.25s peuvent être
superposées à celles obtenues lorsque δt = 2.5s et δt′ = 0.25s. Le gain obtenu au niveau du
temps CPU est considérable puiqu’il est d’un facteur 9.
La figure 4.7 présente le long de l’interface, à trois temps caractéristiques (10s, 100s et 300s),
la surface libre de l’écoulement surfacique (hnh + topographie) ainsi que la valeur moyenne de la
vitesse normale utilisée dans l’algorithme 5 et notée v̄h⋆⋆,n ,
Z
⋆⋆,n 1
∀i ∈ {1 · · · NI }, v̄h |ei = (2v ⋆,n,p + vh⋆,n−1 )/3. (4.15)
li ei h
La surface libre étant constante par morceaux, elle est représentée sur chaque arête par une ligne
horizontale. La vitesse normale moyenne est dessinée avec un cercle bleu (sombre) si la face est
mouillée (F ∈ Ihw,n ) et avec un cercle vert (clair) si la face est sèche (F ∈ Ihd,n ).
La figure 4.8 confirme que les couples charge hydraulique/hauteur d’eau restent sur l’en-
semble admissible A2 . Pour les mêmes temps que ceux de la figure 4.7 et pour l’ensemble des
arêtes de Ih, chaque couple (ψhn , hnh ) est representé par une croix (la valeur moyenne de ψhn est
considérée sur chaque arête). L’ensemble admissible A2 est aussi tracé en trait continu.
Sur la figure 4.7, le ressaut hydraulique qui apparaı̂t dans l’écoulement à surface libre au
début de la simulation est clairement visible au niveau du changement de pente. De plus, une
exfiltration apparaı̂t sur quelques arêtes situées sur I2 et I3 . Pendant la simulation, une condition
aux limites de type Neumann est imposée progressivement sur les faces où la hauteur d’eau
devient nulle. Cela est visible sur la figure 4.8 dans laquelle le nombre de couples (ψhn , hnh ) sur
la branche {h = 0} augmente au cours du temps.
4.2. Résultats 73
10 s 10 s
0.6 0.25
0.5
0.2
0.4
0.15
0.3
0.2 0.1
0.1
0.05
0
0
-0.1
-0.2 -0.05
0 50 100 150 200 250 300 -0.2 -0.1 0 0.1 0.2
100 s 100 s
0.6 0.25
0.5
0.2
0.4
0.15
0.3
0.2 0.1
0.1
0.05
0
0
-0.1
-0.2 -0.05
0 50 100 150 200 250 300 -0.2 -0.1 0 0.1 0.2
300 s 300 s
0.6 0.25
0.5
0.2
0.4
0.15
0.3
0.2 0.1
0.1
0.05
0
0
-0.1
-0.2 -0.05
0 50 100 150 200 250 300 -0.2 -0.1 0 0.1 0.2
où L est la longueur du domaine, P l’altitude du point le plus bas de I et J la pente constante
de l’interface I (voir la figure 4.9). La condition initiale est hydrostatique avec une nappe située
à la cote z = 0.85cm et les parois et le fond sont imperméables,
ψ 0 = −z + 0.85m dans Ω,
vN = 0 sur (W ∪ B) × [0, T ].
h0 = 0 sur I,
hA = 0 en A × [0, T ].
Une pluie d’intensité constante égale à 10% de la conductivité hydraulique à saturation est
imposée pendant les trois premières minutes,
J = 0.5%
▽
1m
0.85m
vN = 0
•
(0, 0) 6m
Un maillage quasi-uniforme avec 2049 triangles (ce qui correspond à des arêtes d’environ
10cm de longueur) et des pas de temps δt = δt′ = 1s sont utilisés. Nous avons vérifié que les
vitesses normales à l’interface obtenues avec un maillage plus fin (8763 elements) et des pas de
temps plus petits (δt = δt′ = 0.5s) peuvent être superposées à celles présentées ci-après. Le pas
de temps δt = 1s est proche de la condition CFL, impliquant que le pas de temps utilisé dans
l’équation de Richards pour des raisons de précision est proche de celui imposé dans l’équation
de l’onde cinématique par la condition CFL pour des raisons de stabilité.
La figure 4.10 présente la surface libre de l’écoulement surfacique et la vitesse normale le long
de l’interface à quatre temps caractéristiques de la simulation (10s, 60s, 180s et 360s). Les mêmes
notations sont utilisées que sur la figure 4.7.
4.2. Résultats 75
10 s 60 s
4 4
3 3
2 2
1 1
0 0
-1 -1
-2 -2
0 100 200 300 400 500 600 0 100 200 300 400 500 600
180 s 360 s
4 4
3 3
2 2
1 1
0 0
-1 -1
-2 -2
0 100 200 300 400 500 600 0 100 200 300 400 500 600
La figure 4.11 montre que les différents couples charge hydraulique/hauteur d’eau restent sur
l’ensemble admissible A2 au cours de la simulation. Pour les mêmes temps que ceux de la figure
4.10 et pour chaque arête de l’interface, chaque couple (ψhn , hnh ) est représenté comme à la figure
4.8. Des échelles différentes étant utilisées pour h et pour ψ, la branche {h = ψ} semble presque
verticale. Les quatre phases précédemment décrites sont clairement illustrées par les positions
successives des nuages de points :
1 - A 10s, la charge hydraulique est négative et la hauteur d’eau est nulle. Le nuage de points
est seulement sur la branche {h = 0} traduisant un état totalement sec.
2 - A 60s, la charge hydraulique est égale à la hauteur d’eau pour quelques faces. Le nuage
de points est situé sur les deux branches puisque le sol contient des zones saturée et insaturée.
3 - A 180s, la charge hydraulique est égale à la hauteur d’eau pour toutes les faces. Le nuage
de points est seulement sur la branche {h = ψ} traduisant un état totalement mouillé.
4 - A 360s, la charge hydraulique redevient négative sur les faces où la hauteur d’eau est
nulle. De nouveau, le nuage de points est situé sur les deux branches.
10 s 60 s
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-0.01 -0.01
-14 -12 -10 -8 -6 -4 -2 0 -14 -12 -10 -8 -6 -4 -2 0
180 s 360 s
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-0.01 -0.01
-14 -12 -10 -8 -6 -4 -2 0 -14 -12 -10 -8 -6 -4 -2 0
La grandeur ∆Mn est la variation totale de masse pendant l’intervalle [0, nδt] et se décompose
naturellement comme la somme de la variation totale d’eau située dans le sol ∆Mngrnd et de la
P P
variation totale d’eau située sur le sol ∆Mnover . Les deux quantités ni=1 Miin et ni=1 Miout sont
respectivement la masse d’eau qui entre et qui sort du système. L’erreur commise sur la P conserva-
n . Les cinq quantités ∆Mn , ∆Mn n , n i
tionPnde la masse entre 0 et nδt est notée E grnd , ∆M over i=1 Min
et i=1 Mnout sont représentées sur la partie gauche de la figure 4.12. Cette figure confirme les
quatre phases de la simulation. En effet la pluie est totalement absorbée par le sol jusqu’à 50s
car ∆Mn = ∆Mngrnd . Ensuite l’augmentation de ∆Mngrnd fléchit et ∆Mnover devient positive suite
à la saturation d’une partie de la surface du sol. De 90s à la fin de la simulation, les variations
de ∆Mn et ∆Mnover sont les mêmes, indiquant une saturation complète du sol. Enfin, pendant
la dernière phase, la masse d’eau infiltrée est constante (la pluie ayant cessé) et la masse d’eau
sortant du système est égale à la variation totale de masse.
Afin d’étudier l’influence du nouveau flux à l’interface dans l’algorithme de couplage à deux
pas, nous définissons l’algorithme 4′ comme l’algorithme 4 dans lequel l’équation de Richards
est résolue par la BDF2. Les erreurs commises sur la conservation de la masse obtenues avec
l’algorithme 5 et l’algorithme 4′ sont comparées sur la figure 4.13. L’algorithme 5 est bien plus
précis que l’algorithme 4′ . L’erreur produite par l’algorithme 4′ est néanmoins faible en compa-
raison des quantités totales d’eau puisqu’elle est seulement de l’ordre de quelques pourcents des
masses d’eau considérées.
12
ΣMn
in 0.03
8
∆Mn
Masse [kg]
Masse [kg]
0.02
4
∆Mn
grnd
0 ∆Mn
over 0.01
−4
0
ΣMn
out
−8
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350
Temps [s] Temps [s]
Fig. 4.13: Erreur En
Fig. 4.12: Bilan de masse. pour l’algorithme 4′ ( --- ) et 5 ( ---- ).
Enfin, la figure 4.14 détaille les flux de masse dans l’équation de l’onde cinématique. Comme
dans la section 3.2.1, le flux de masse F̃In est décomposé en un flux d’exfiltration F̃I,out n et un
n n n n n
flux d’infiltration F̃I,in . Les quatres quantités ρδtF̃I,in , ρδtF̃I,out , Min et Mout sont tracées sur la
figure 4.14 en fonction du temps. Nous voyons sur cette figure que les quantités d’eau infiltrée
et exfiltrée deviennent faibles lorsque le milieu est saturé.
78 Chapitre 4. Couplage avec le ruissellement surfacique
0.08
Mn
in
0.04
Masse [kg]
n
ρδtF̃I,in
0
n
ρδtF̃I,out
−0.04
Mn
out
−0.08
0 50 100 150 200 250 300 350
Temps [s]
avec les notations usuelles (voir figure 4.15). La condition initiale est hydrostatique avec une
nappe située à 0.1m et la condition aux limites sur les parois est un flux nul,
ψ 0 = −z + 0.1m dans Ω,
vN = 0 sur W × [0, T ].
Un flux d’infiltration avec un profil parabolique et une valeur moyenne v̄N égale à 5% de la
conductivité hydraulique à saturation est imposé pendant 2 minutes sur la partie gauche du
fond (Bl = {(x, z) ∈ B, x ∈ [0, 1]}). Ce flux d’infiltration est linéaire en temps pendant les 10
premières secondes puis constant,
−1 si (x, t) ∈ B × [0, 10],
x(x − 1) 0.03Ks t cm.s
l
vN (x, t) = x(x − 1) 0.3Ks cm.s −1 si (x, t) ∈ Bl × [10, 120],
0, sinon.
h0 = 0 sur I,
hA = 0 en A × [0, T ].
Un maillage quasi-uniforme avec 1874 triangles (ce qui correspond à des arêtes d’environ 2.5cm
de longueur) et des pas de temps δt = δt′ = 1s sont utilisés. Nous avons vérifié que les vitesses
normales à l’interface obtenues avec un maillage plus fin (7310 éléments) et des pas de temps
plus petits (δt = δt′ = 0.5s) peuvent être superposées à celles présentées ci-après. A nouveau, la
précision de la résolution de l’équation de Richards impose un pas de temps comparable à celui
donné par la condition CFL.
4.2. Résultats 79
J = 0.2%
v̄N = 0.05Ks
▽
0.2m
0.1m
vN = 0
•
(0, 0) 1m 1m
Fig. 4.15: Géométrie et condition initiale.
5s 35 s
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-0.1 -0.1
-0.2 -0.2
0 50 100 150 200 0 50 100 150 200
50 s 100 s
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-0.1 -0.1
-0.2 -0.2
0 50 100 150 200 0 50 100 150 200
150 s 360 s
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
-0.1 -0.1
-0.2 -0.2
0 50 100 150 200 0 50 100 150 200
5s 35 s
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-0.01 -0.01
-10 -8 -6 -4 -2 0 -10 -8 -6 -4 -2 0
50 s 100 s
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-0.01 -0.01
-10 -8 -6 -4 -2 0 -10 -8 -6 -4 -2 0
150 s 360 s
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
-0.01 -0.01
-10 -8 -6 -4 -2 0 -10 -8 -6 -4 -2 0
−3
x 10
0.6
ΣMn
in
3
0.4
∆Mn
Masse [kg]
Masse [kg]
0.2 2
∆Mn
grnd
0 ∆Mn
over
1
−0.2
−0.4
ΣMn
out
0
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350
Temps [s] Temps [s]
Fig. 4.19: Erreur En
Fig. 4.18: Bilan de masse. pour l’algorithme 4′ ( --- ) et 5 ( ---- ).
−3
x 10
6 0.05
Mn
in
4
0.04
150s
Hauteur d’eau hh [cm]
n 100s
ρδtF̃I,in
2
Masse [kg]
0.03
0
n
ρδtF̃I,out
0.02 50s
−2
Mn
out
−6 0
0 50 100 150 200 250 300 350 0 40 80 120 160 200
Temps [s] x [cm]
Fig. 4.20: Flux dans l’onde cinématique. Fig. 4.21: Hauteur d’eau à différents temps.
Les résultats des figures 4.18, 4.19 et 4.20 sont proches de ceux du cas précédents, en par-
ticulier dans la comparaison des erreurs sur la conservation de la masse avec (algorithme 5) et
sans (algorithme 4′ ) les modifications du flux à l’interface. La figure 4.21 présente la hauteur
d’eau hnh à différents instants et permet d’avoir une vision plus directe de l’évolution temporelle
de la surface libre.
avec les notations usuelles (voir la figure 4.22). La condition initiale est hydrostatique avec une
nappe située à la cote z = 0.6cm et des conditions aux limites d’imperméabilité sont imposées
4.2. Résultats 83
ψ 0 = −z + 0.6m dans Ω,
vN = 0 sur (W ∪ B) × [0, T ].
Pour l’écoulement surfacique, la condition initiale est une hauteur d’eau nulle,
h0 = 0 sur I.
Il est clair que cette simulation ne peut pas être prolongée en temps long sans violer l’hypothèse
de validité de l’équation de Richards selon laquelle la pression de l’air dans le milieu poreux
est constante et n’influence pas l’écoulement de l’eau. En effet, en continant la simulation, la
pression de l’eau augmenterait sans rencontrer la résistance de l’air emprisonné ce qui est en
contradiction avec l’expérience physique.
2cm
▽
1m
0.6m
vN = 0
•
(0, 0) 6m
Le domaine Ω étant similaire à celui utilisé dans le cas test d’exfiltration par la pluie de
la section 4.2.2, nous utilisons le même maillage. Le pas de temps dans la discrétisation de
l’équation de Richards est δt = 1s et celui dans l’équation de l’onde cinématique est de δt′ = 0.1s,
la condition CFL imposant un pas de temps inférieur à 0.15s.
La figure 4.23 présente la surface libre de l’écoulement et la vitesse normale le long de
l’interface à 5s, 15s et 60s. A 5s, l’eau introduite au début de la simulation par la condition à
la limite en A est totalement absorbée par le sol puisque les conditions sont de type Neumann
sur l’ensemble de l’interface. Les conditions aux limites deviennent progressivement de type
Dirichlet en suivant l’avancement du front de ruissellement. Il est clairement visible à 15s que
les conditions aux limites sont de type Neumann homogène sur les arêtes où la hauteur d’eau est
nulle. Les nuages de points de la figure 4.24 montrent que la continuité de la pression est vérifiée
sur l’ensemble de la zone mouillée de l’interface même si certains couples situés probablement
au niveau du front s’éloignent légérement de l’ensemble admissible A2 .
Le bilan de masse de la figure 4.25 permet de vérifier que le débit entrant augmente pendant
les dix premières secondes puis est constant dans le temps car la somme de la masse apportée au
84 Chapitre 4. Couplage avec le ruissellement surfacique
5s 5s
6 2
5
1.5
4
3
1
0.5
1
0
0
-1
-2 -0.5
0 100 200 300 400 500 600 -45 -40 -35 -30 -25 -20 -15 -10 -5 0
15 s 15 s
6 2
5
1.5
4
3
1
0.5
1
0
0
-1
-2 -0.5
0 100 200 300 400 500 600 -45 -40 -35 -30 -25 -20 -15 -10 -5 0
60 s 60 s
6 2
5
1.5
4
3
1
0.5
1
0
0
-1
-2 -0.5
0 100 200 300 400 500 600 -45 -40 -35 -30 -25 -20 -15 -10 -5 0
système est linéaire après 10s. Nous voyons également qu’il y a égalité entre la masse apportée
et l’augmentation de masse dans le système tant que l’eau en surface n’a pas atteint l’exutoire.
A partir de cet instant (42s), le débit à l’exutoire devient non nul et la masse d’eau totale dans
le système augmente moins vite. Comme précédemment, la comparaison des erreurs sur le bilan
de masse (figure 4.26) montre l’amélioration obtenue sur la conservation de la masse en adaptant
le terme de couplage dans l’équation de l’onde cinématique.
0.2
ΣMn
in
300 0
∆Mn −0.2
200
Masse [kg]
Masse [kg]
−0.4
∆Mn
grnd
−0.6
100
∆Mn
over −0.8
−1
0
−1.2
ΣMn
out
−100
0 10 20 30 40 50 60 0 10 20 30 40 50 60
Temps [s] Temps [s]
Fig. 4.26: Erreur En
Fig. 4.25: Bilan de masse. pour l’algorithme 4′ ( --- ) et 5 ( ---- ).
La figure 4.27 concernant la répartition des masses au niveau de l’onde cinématique montre
qu’environ la moitié de l’eau apportée ruisselle et que l’autre moitié s’infiltre. Lorsque l’eau de
surface atteint l’exutoire, cette répartition change : la masse d’eau infiltrée dans le sol augmente
moins vite puisque le milieu se sature et la masse d’eau de surface augmente peu puisqu’il y a
ruissellement sur toutes les mailles de surface. Comme nous l’avons déjà précisé, cette situation
n’est pas extrapolable en temps long sans rompre l’hypothèse que l’air dans le sol est à la pression
atmosphérique. La figure 4.28 présente la hauteur d’eau hnh à differents instants.
2
6 Mn
in
60s
1.6
30s
4 15s
Hauteur d’eau hh [cm]
Masse [kg]
2 1.2
ρδtF̃I+
0
0.8
−2 ρδtF̃I−
0.4
−4
Mn
out
−6 0
0 10 20 30 40 50 60 0 100 200 300 400 500 600
Temps [s] x [cm]
Fig. 4.27: Flux dans l’onde cinématique. Fig. 4.28: Hauteur d’eau à différents temps.
86 Chapitre 4. Couplage avec le ruissellement surfacique
∂t [θ(ψ)] + ∇ · v(ψ) = 0 dans Ω × [0, T ],
v(ψ) = −K(ψ)∇(ψ + z) dans Ω × [0, T ],
v(ψ) · nΩ = 0 sur (W ∪ B) × [0, T ],
v(ψ) · nΩ = ωv sur {(x, t), x ∈ I d,t },
ψ = ωψ sur {(x, t), x ∈ I w,t },
v(ψ) · nΩ = 0 sur {(x, t), x ∈ Dd,t },
ψ=0 sur {(x, t), x ∈ Dw,t },
(4.16)
ψ(·, 0) = ψ 0 sur Ω,
∂t h + ∂x ϕ(h, S) = (v(ψ) − vr ) · nΩ sur I × [0, T ],
h(·, 0) = h0
sur I,
h(A, ·) = 0 en A × [0, T ],
(ψ, v(ψ) · nΩ ) ∈ A1
sur D × [0, T ],
(ψ, h) ∈ A2 sur I × [0, T ].
4.3. Application concrète en hydrologie 87
Les données sont les conditions initiales ψ 0 et h0 , le flux normal imposé sur les parois et le
fond du domaine ainsi que la hauteur d’eau imposée en A. L’algorithme mis en place est une
généralisation de l’algorithme 5 en prenant en compte les éléments de l’algorithme 3 concernant
la prise en compte de drains.
vN = 0 sur (W ∪ B) × [0, T ].
d1 d2 d3 d4 1m
(0, 0) •
4m 8m 8m 8m 4m
Fig. 4.29: Géométrie de la parcelle drainée.
Les propriétes hydrodynamiques du sol, représentées sur la figure 4.30, sont données par le
modèle de Van Genuchten modifié (voir les équations (2.31)) avec les valeurs suivantes pour les
différents paramètres,
−4
x 10
0.44
2
0.42
0.41
1
0.4
0.39 0
−100 −80 −60 −40 −20 0 20 −100 −80 −60 −40 −20 0 20
Charge hydraulique ψ [cm] Charge hydraulique ψ [cm]
Fig. 4.30: Propriétés hydrodynamiques du sol dans le cas test avec plusieurs drains.
Le coefficient de Strickler K est égal à 30m1/3 .s−1 ce qui correspond à un frottement sur une
zone cultivée (pâturage par exemple).
4.3.3 Résultats
Influence de la condition initiale
Deux simulations sont réalisées en faisant varier uniquement la condition intiale afin d’étudier
l’influence de cette donnée sur le comportement hydrologique du système. Le temps final de la
simulation et les pas de temps utilisés sont
Une forte intensité de pluie est imposée sur la surface du sol pendant toute la simulation
Pour la première condition initiale (CI 1), les isovaleurs de la charge hydraulique sont horizontales
alors que pour la seconde condition initiale (CI 2), les isovaleurs de la charge hydraulique sont
parallèles à la surface du sol,
CI 1 : ψ 0 = −0.5z dans Ω,
0
CI 2 : ψ = −0.5(z + 0.005(32 − x)) dans Ω.
Les figures 4.31, 4.32 et 4.33 représentent les positions de la nappe obtenues avec les deux
conditions initiales à 1h, 1h15 et 1h30. Ces figures illustrent la sensibilité du système à la
condition initiale. L’échelle verticale est multipliée par 4 pour avoir une meilleure visualisation.
4.3. Application concrète en hydrologie 89
Les positions de la nappe entre les drains à 1h de simulation sont clairement dans la continuité
de la condition initiale : la nappe est horizontale dans le premier cas alors qu’elle est parallèle
au fond dans le second cas. A 1h15, la nappe monte plus vite au niveau du bord droit que dans
le reste du domaine avec la première condition initiale car la conductivité à la surface du sol
est d’autant plus importante que l’on s’éloigne du bord gauche. En revanche, la montée de la
nappe s’effectue de la même façon dans la totalité du domaine avec la seconde condition initiale.
La différence de position entre les deux nappes est accentuée à 1h30.
3600s
400
300
200
100
0
3600s
400
300
200
100
0
Le principe pour obtenir le bilan de masse est identique aux cas tests précédents. Un terme
supplémentaire dû à la la présence des drains dans le domaine apparaı̂t dans ce bilan qui s’écrit
n
X Xn
n 0 n 0
ρ(Vgrnd − Vgrnd ) + ρ(Vover − Vover )= ρδtFri + ρδt(FBi + FDi ) + ρ∆V n .
| {z } | {z } i=1 | {z } i=1 | {z } | {z }
∆Mn
grnd ∆Mn
over Miin Miout En
90 Chapitre 4. Couplage avec le ruissellement surfacique
La quantité FDi est le flux d’eau passant par les quatre drains sur l’intervalle [(i − 1)δt, iδt],
def
FDi = FDi 1 + FDi 2 + FDi 3 + FDi 4 .
Nous rappelons que le flux d’eau autour d’un drain est calculé en intégrant −vh⋆ sur le bord du
drain. Les figures 4.34 et 4.35 représentent pour la CI 1 et la CI 2 respectivementP les quantités
4
MnDj = ρδtFDnj (j étant un entier compris entre 1 et 4) et la quantité MnD = n
j=1 ρδtFDj .
La figure 4.34 montre que les drains n’arrivent pas au régime maximum en même temps avec
la première condition initiale. Un décalage d’environ 8 minutes est observé entre chaque drain,
conduisant à un écart d’environ 24 minutes entre les drains d1 et d4 . Ce résultat confirme
que le sol commence à saturer dans la partie droite puisque le drain d4 est actif le premier.
Par contre, la figure 4.35 montre une synchronisation parfaite entre les différents drains avec la
seconde condition initiale. Nous voyons également sur cette figure que la masse totale qui sort
du système est identique avec les deux conditions initiales quand les quatre drains fonctionnent
au régime maximum.
0.14 0.14
D D
0.12 0.12
Masse [kg]
0.08 0.08
0.06 0.06
0.04 D4 D 3 D 2 D 1 0.04 D4 D3 D2 D1
0.02 0.02
0 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Temps [h] Temps [h]
Fig. 4.34: Bilan de masse aux drains - CI 1. Fig. 4.35: Bilan de masse aux drains - CI 2.
La figure 4.36 représente la position de la nappe à deux instants où l’affleurement commence
à se produire (t = 6000s pour la figure du haut et t = 6150s pour la figure du bas). Comme
attendu, la figure du haut montre que la nappe débute l’affleurement aux milieux des zones
inter-drains. La symétrie des zones insaturées situées au dessus des drains et à proximité de
la surface est également visible sur cette figure. Nous pouvons observer que cette symétrie est
brisée sur la figure du bas ce qui met en évidence le ruissellement sur la surface du sol. Il est
clair que ce phénomène ne peut pas être détecté sans un modèle de couplage prenant en compte
le ruissellement surfacique.
4.3. Application concrète en hydrologie 91
6000s
400
300
200
100
0
6150s
400
300
200
100
0
Fig. 4.36: Affleurement de la nappe aux instants t = 6000s (haut) et t = 6150s (bas).
La figure 4.37 représente la hauteur d’eau à la surface du sol à quatre instants de la simula-
tion : t1 = 1h40, t2 = t1 + 7min30s, t3 = t1 + 15min et t4 = t1 + 1h20min. Notons que t1 désigne
le temps où l’affleurement débute et que t4 correspond au temps final de simulation (t4 = T ).
Sur cette figure, la hauteur entre le bord inférieur et le bord supérieur du cadre représente
1.5mm. La hauteur d’eau à t1 est en alternance strictement positive puis nulle. Cela montre que
les conditions aux limites sur la surface du sol sont successivement de type Dirichlet (h > 0)
et de type Neumann (h = 0) et que notre algorithme est capable de traiter plusieurs zones
d’affleurement en même temps. La hauteur d’eau à t2 est positive sur l’ensemble de l’interface
ce qui montre que notre algorithme a su traiter la connexion de plusieurs zones de ruissellement.
La hauteur d’eau à t3 montre qu’après un quart d’heure d’exfiltration en surface, la hauteur
d’eau est assez proche de sa valeur finale.
t4
t3
t2
t1
La figure 4.38 représente les quantités ΣMnin , ∆Mn , ∆Mngrnd , ∆Mnover et ΣMnout . Cette figure
confirme que l’eau introduite par la pluie est totalement absorbée par le sol au début de la simu-
lation puisque la masse d’eau à la surface du sol est nulle jusqu’à environ 1h45min. Pendant cette
phase de saturation du sol, le drainage augmente puisque l’écart entre la courbe représentant
∆Mngrnd et celle représentant ΣMnin s’accroı̂t progressivement (nous retrouvons le résultat de la
figure 4.35). Une phase de transition d’environ un quart d’heure est ensuite visible : l’affleure-
ment se produit puisque la quantité ∆Mnover devient non nulle. Le sol se sature complètement
car la quantité ∆Mngrnd devient constante. Enfin, un état stationnaire semble être atteint après
deux heures de simulation : les masses d’eau contenues dans le sol et dans la lame ruisselante
semblent être constantes, la masse d’eau introduite dans le système étant alors égale à la masse
d’eau sortant du système.
92 Chapitre 4. Couplage avec le ruissellement surfacique
50
500 ΣMn
in
400 0
Masse [kg]
Masse [kg]
∆Mn
grnd −50
200
100 ΣMn
out,B
∆Mn
over
−100
0
−100 −150
−200 ΣMn
out ΣMn
out
−200
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Temps [h] Temps [h]
Fig. 4.39: Décomposition de la masse d’eau
Fig. 4.38: Bilan de masse. sortant du système
La masse totale sortie du sytème à l’instant nδt, notée ΣMnout , peut être décomposée suivant
que l’eau est sortie au niveau de l’exutoire ou au niveau des drains :
n
X n
X n
X
ΣMnout = ρδt(FBi + FDi ) = ρδtFBi + ρδtFDi
i=1 i=1 i=1
| {z } | {z }
ΣMn
out,B ΣMn
out,D
La figure 4.39 représente les trois quantités ΣMnout , ΣMnout,B et ΣMnout,D . Comme attendu, cette
figure montre que la masse d’eau qui sort du système est uniquement évacuée par les drains
jusqu’à une 1h45min de simulation. A partir de cet instant, l’eau commence à sortir à l’exutoire
et il est intéressant de noter que la pente de la droite représentant la quantité ΣMnout,B est
beaucoup plus importante que celle représentant la quantité ΣMnout,D . Cela signifie que les drains
retardent l’affleurement de la nappe mais que le débit à l’exutoire est supérieur au débit des
drains lorsque le milieu est complètement saturé. In fine, la masse totale qui transite à l’exutoire
est environ le triple de celle évacuée par les drains.
0
0 5 10 15 20 25 30 35 40 45
Temps [h]
La figure 4.41 représente les hauteurs de la nappe à l’interdrain obtenues par les mesures
expérimentales et la simulation. Cette comparaison montre que le sol modélisé est plus conduc-
teur que le sol réel car les affleurements de la nappe issus du calcul se produisent avant les
affleurements mesurés.
120
110
100
90
hexp
80
70
60
hnum
50
0 5 10 15 20 25 30 35 40 45
Temps
Conclusion
Ce cas test de simulation d’une partie d’un bassin versant expérimental nous a permis de
tester notre logiciel R+R DG dans une configuration réaliste. La première série de résultats
est encourageante puisque nous pouvons traiter l’apparition simultanée de plusieurs zones de
ruissellement ainsi que leurs connexions. La propagation du ruissellement est clairement visible
par la dissymétrie de la zone saturée sur la figure 4.36. D’un point de vue hydrologique, ce cas
test a montré l’impact de la condition initiale sur le comportement du système ce qui motive la
nécessité de bien connaı̂tre cette donnée d’entrée du modèle. Par ailleurs, le fait que la masse
d’eau totale passant par l’exutoire est supérieure à la masse d’eau évacuée par les drains traduit
que le réseau de drainage étudié retarde l’affleurement mais n’empêche pas le ruissellement
de grandes quantités d’eau. La seconde série de résultats montre des écarts significatifs entre
les résultats numériques et les mesures expérimentales. Toutefois, la réalisation d’autres tests en
changeant certains éléments de modélisation (rugosité du terrain et propriétés hydrodynamiques
du sol) devrait nous permettre d’obtenir des résultats plus proches de ceux mesurés sur le terrain.
Conclusion
Dans cette thèse, nous avons étudié la modélisation numérique des écoulements souterrains
en utilisant une méthode de Galerkine discontinue pour la discrétisation en espace et des schémas
BDF pour la discrétisation en temps. Nous avons aussi étudié le couplage de ces écoulements
avec le ruissellement surfacique discrétisé par une méthode de volumes finis et un schéma d’Euler
explicite.
Dans le premier chapitre, l’intérêt de ces travaux en hydrologie a été situé. Les différents
modèles disponibles pour décrire les deux types d’écoulements ont été présentés. Nous avons
choisi l’équation de Richards pour modéliser les écoulements dans le sol et l’approximation de
l’onde cinématique pour décrire les écoulements à la surface du sol.
Le chapitre 2 a permis de présenter la méthode de Galerkine discontinue SIPG, les schémas de
discrétisation temporelle ainsi que la linéarisation effectuée. Deux schémas ont été implémentés
pour l’initialisation : un schéma de Runge–Kutta diagonalement implicite utilisé lorsque le milieu
est insaturé à l’instant initial et le schéma de Crank–Nicolson utilisé dans le cas contraire.
Le chapitre 3 a introduit les problèmes de frontière libre en résolvant l’équation de Richards
avec des conditions aux limites de Signorini. Un algorithme pour traiter ce type de condition a
été proposé puis testé dans le cas de l’affleurement d’une nappe et dans le cas de la présence d’un
drain dans le sol. Des propriétés de conservation de la masse d’eau dans le sol ont été établies
pour les schémas BDF d’ordre 1 et 2.
Le chapitre 4 a constitué la partie principale de nos travaux : le couplage des écoulements
souterrains et superficiels en prenant en compte la propagation d’une lame d’eau ruisselante.
Deux algorithmes de couplage, l’un avec la BDF1 et l’autre avec la BDF2, ont été proposés.
Le choix du terme de couplage dans les écoulements superficiels est important pour assurer la
conservation de la masse lorsque l’ordre de la BDF est supérieur ou égal à deux.
La méthodologie numérique développée dans cette thèse est précise (conservation de la masse
d’eau, taux de convergence optimaux sur solution analytique), flexible (possibilité d’utiliser des
maillages non-coı̈ncidants, extension aisée à la dimension 3, pas de temps différents en subsur-
face et en surface) et robuste (possibilité de traiter les zones de transitions saturées/insaturées
du sol, l’affleurement de nappe, le ruissellement hortonien et le drainage, une pluie variable
dans le temps). Le logiciel R+R DG que nous avons développé est ainsi à même d’être utilisé
pour étudier, certes dans des configurations encore simplifiées, le fonctionnement hydrologique
de petits bassins versants drainés en milieu agricole. L’identification des conditions initiales d’un
système est un exemple d’études actuellement réalisables. La prise en compte de l’hétérogénéité
en espace et éventuellement en temps du terrain serait aussi un élément à étudier par un traite-
ment quantitatif des incertitudes par une approche déterministe (analyse de sensibilité, modèle
inverse) ou stochastique (polynômes de chaos).
Annexe A
Le but de cette annexe est de présenter le couplage des écoulements souterrains avec les
écoulements surfaciques qui sont modélisés par les équations de Saint–Venant. Nous commençons
par exposer la discrétisation de ce système : le flux numérique est précisé, quatre types de
conditions aux limites peuvent être imposés et une modification du flux numérique permet
d’obtenir un schéma équilibré. Nous présentons également un cas test de transfert d’eau entre
plusieurs sillons qui n’est pas réalisable avec l’approximation de l’onde cinématique. En effet
celle-ci ne peut pas modéliser la dynamique de l’écoulement dans un sillon puisque l’hypothèse
d’équilibre entre gravité et frottements implique un écoulement des zones ayant l’altitude la plus
grande vers les zones ayant une altitude plus faible.
def
La vitesse du fluide est définie comme le rapport u = q/h. Nous rappelons également que les
def √
valeurs propres de la matrice jacobienne ∂W F sont u±c où c = gh est la vitesse de propagation
des ondes.
La discrétisation des équations de Saint–Venant s’effectue avec les volumes finis mais une
méthode DG, comme le font par exemple Aizinger et Dawson [AD02], Tassi, Bokhove et Vion-
net [TBV07] ou Ern, Piperno et Djadel [EPD08], est également envisageable. Comme pour la
discrétisation de l’équation de l’onde cinématique, nous considérons la trace du maillage Th sur
l’interface I. Les notations sont également identiques : xi , li , xi− 1 et xi+ 1 sont respectivement
2 2
le centre, la longueur, et les sommets gauche et droit d’une arête quelconque ei du maillage de I
(voir la figure A.1). La quantité zi désigne la cote du centre de l’arête ei .
98 Annexe A. Couplage avec les équations de Saint–Venant
Le système régissant les écoulements superficiels est composé du système (A.1), d’une condi-
tion initiale W0 ainsi que de conditions aux limites qui peuvent être imposées aux extrémités A
et B de l’interface et qui seront précisées ci-après.
hi−1 hi+1
hi
x− 1 xi− 3 xNI + 1
•
2
•
2
(xi , zi ) xi+ 1 •
2
xi− 1 2
A 2
× • xi+ 3 B
• 2
•
li
• Flux HLLE
Plusieurs flux numériques peuvent être utilisés dans la discrétisation des équations de Saint–
Venant. Le flux de Harten, Lax et Van Leer [HLvL83], noté HLL, est une généralisation du
flux de Lax–Friedrichs. Le flux HLLE proposé par Einfeldt [Ein88] est une amélioration du flux
HLL puisqu’il préserve la positivité de la hauteur moyenne d’eau dans le cas monodimensionnel.
De même, le flux HLLC proposé par Toro [Tor01] assure la positivité de la hauteur moyenne en
dimension deux d’espace. Dans notre cas, nous choisissons le flux HLLE avec un schéma d’Euler
explicite en temps. Pour simplifier la présentation, nous supposons que le pas de temps δt′ utilisé
pour le système de Saint–Venant est égal au pas de temps δt choisi pour l’équation de Richards.
Sous cette hypothèse, le schéma s’écrit
Z
Wn+1 − Win n vh⋆,n+1
li i = F(Wi−1 , Win ) − F(Win , Wi+1
n
)+ , (A.2)
δt′ ei ghSi
où F(WL , WR ) est le flux HLLE. Ce flux est basé sur le fait que l’approximation de la solution
repose sur trois états appelés WL , W0 , WR et séparés par deux vitesses de propagation des
ondes c+ et c− . Partant des deux états
hL hR
WL = et WR = ,
qL qR
Les vitesses de propagation des ondes c+ et c− sont calculées à partir des quantités précédentes,
def def
c+ = max(0, uR + cR , u0 + c0 ) et c− = min(0, uL − cL , u0 − c0 ).
A.1. Discrétisation des équations de Saint–Venant 99
Le flux HLLE est la somme d’une partie centrée et d’une partie décentrée,
def 1 1
F = F(WL ) + F(WR ) + Q(WL − WR ),
2 2
où la matrice d’ordre deux Q est définie par
def c+ + c− 0 1 c+ c− 1 0
Q = −2 .
c+ − c − −u20 + gh0 2u0 c+ − c − 0 1
Pour éviter de modifier le pas de temps δt′ au cours de la simulation, nous remplaçons cette
condition CFL par
1
δt′ = √ · min li ,
2 ghmax 1≤i≤NI
dans laquelle hmax est une borne supérieure a priori de la hauteur d’eau h sur I × [0, T ].
soit
√ 1 √
WA = hA .
2 ghA + q1 /h1 − 2 gh1
Au point B avec une hauteur hB , on utilise l’état WNI , ce qui donne
1
WB = hB √ p .
−2 ghB + qNI /hNI + 2 ghNI
100 Annexe A. Couplage avec les équations de Saint–Venant
– Un débit est également imposé en conservant les invariants de Riemann. Au point A avec un
débit qA , on obtient
qA p q1 p p q p
1
− 2 ghA = − 2 gh1 ⇒ −2 ghA hA + hA − + 2 gh1 + qA = 0, (A.4)
hA h1 h1
soit
hA
WA = ,
qA
• Schéma équilibré
La préservation d’états d’équilibre est une caractéristique utile des schémas approchant les
équations de Saint–Venant. Nous considérons en particulier les états stationnaires au repos
caractérisés par les deux conditions
où z désigne la cote du fond et C est une constante. Les schémas préservant de tels états sont dits
équilibrés. La modification suivante du schéma (A.2), proposée par Audusse et al. [ABB+ 04],
permet qu’il soit équilibré en l’absence de terme de couplage,
Wn+1 − Win n n n n g 0
li i = F(Wi−1− , Wi−1+ ) − F(Wi+1− , Wi+1+ ) + , (A.5)
δt ′ 2 (hi+1− ) − (hni−1+ )2
n 2
hi−1+ hi+1−
où Wi−1+ = et Wi+1− = .
qi−1+ qi+1−
∀i ∈ {1 · · · NI }, ∀n ∈ {0 · · · NT }, ∀C ∈ R,
hni ≥ 0, hni + zi = C et qin = 0 ⇒ hn+1
i ≥ 0, hn+1
i + zi = C et qin+1 = 0.
A.2. Transfert entre sillons 101
L = 3m, PA = PB = 1m et T = 2h,
avec les notations usuelles (voir la figure A.2). La condition initiale est hydrostatique avec une
nappe située au niveau du fond des sillons, c’est-à-dire à 0.8m, et la condition aux limites sur
les parois et le fond est un flux nul,
ψ 0 = −z + 0.8m dans Ω,
vN = 0 sur (W ∪ B) × [0, T ].
Pour l’écoulement surfacique, la condition initiale est un remplissage du sillon situé à proximité
du point A, (
0 1m − z si x ≤ 1,
h =
0 si x ≥ 1,
et les conditions aux limites sont
hA = 0 en A × [0, T ],
hB = 0 en B × [0, T ].
Nous pourrons considérer que l’état initial correspond au remplissage instantané du premier
sillon au début de la simulation.
h0 + z = 1m
0.2m ▽
1m
vN = 0
•
(0, 0) 3m
Fig. A.2: Géométrie et condition initiale.
Les propriétés hydrodynamiques du sol sont identiques à celles utilisées dans le chapitre 4, elles
correspondent aux relations d’Haverkamp (2.27) avec les valeurs suivantes pour les différents
paramètres,
θs = 0.5 α̃ = 0.028cm−1 Ks = 10−2 cm.s−1 γ = 4
θr = 0.05 β = 4 Ã = 0.030cm−1
102 Annexe A. Couplage avec les équations de Saint–Venant
Un maillage quasi-uniforme avec 3283 triangles, représenté à la figure A.3, est considéré. Les
pas de temps δt = 2s et δt′ = 0.2s sont utilisés, ce qui conduit à 3600 résolutions de l’équation
de Richards pendant la simulation.
2s 2s
110 20
100 15
90 10
80 5
70 0
60 -5
0 50 100 150 200 250 300 -20 -15 -10 -5 0 5 10 15 20
200 s 200 s
110 20
100 15
90 10
80 5
70 0
60 -5
0 50 100 150 200 250 300 -20 -15 -10 -5 0 5 10 15 20
7200 s 7200 s
110 20
100 15
90 10
80 5
70 0
60 -5
0 50 100 150 200 250 300 -20 -15 -10 -5 0 5 10 15 20
2s
100
200s
100
2h
100
Fig. A.6: Etat de saturation du sol au niveau des sillons à différents instants.
Annexe B
Le logiciel R+R DG
L’objectif de cette annexe est de présenter le logiciel R+R DG que nous avons conçu,
développé en C et testé dans plusieurs configurations. Ce logiciel contient environ 9000 lignes
(hors routines d’algèbre linéaire) et nous décrivons quelques extraits de nos programmes in-
formatiques pour illustrer la façon dont nous avons implémenté la méthode SIPG ainsi que le
couplage. Nous commençons par l’organisation générale du code, puis nous détaillons l’assem-
blage de la matrice et du second membre. Nous terminons par préciser les éléments nécessaires
à la prise en compte du couplage au niveau de la surface du sol.
dinteg 0 2 4 5 8 12 16
N b2d 1 3 6 7 16 33 52
N b1d 1 2 3 4 5 7 9
Données d’entrée
• Propriétés hydrodynamiques du sol : ψ 7→ K(ψ), ψ 7→ θ(ψ) et ψ 7→ C(ψ)
c c
• Maillage Th du domaine Ω : construction avec Freefem
et renumérotation avec Matlab
• CI, CL, pluie et coefficient de Strikler : ψ 0 , vN (x, t), h0 , hA (t), vr (t) et K.
• Eléments de discrétisation : degrés p en espace et q en temps, η, T, δt, δt′ , ǫalg1
• Intégration numérique : degré, poids et points de Gauss (1d/2d), valeurs des fonctions de base
Résolution du problème
• Initialisation de la charge et de la hauteur d’eau : ψ 0 → ψh0 et h0 → h0h
• Algorithme de couplage à un pas avec le schéma DIRK3 ou CN2 (Algorithme 2)
• Boucle en temps (BDF2)
• Algorithme de couplage à deux pas (Algorithme 5)
◦ Résolution de l’onde cinématique sans terme de couplage
◦ Boucle itérative
− Partition de l’interface et des drains
− Détermination des conditions aux limites
− Résolution de l’équation de Richards (Algorithme 1)
− Evaluation de la vitesse normale vh⋆
− Calcul de la hauteur d’eau (prise en compte de vh⋆ )
◦ Mise à jour de la charge ψh et de la hauteur d’eau hh
c c
• Stockage et tracé de la charge et de la hauteur d’eau avec Gnuplot
et Plotmtv
• Calcul des masses d’eau et des flux pour obtenir le bilan de masse
Données de sortie
c
• Bilan de masse avec Matlab
• Animations : isovaleurs de la charge, hauteur d’eau et nuages de points
L’étape qui a nécessité le plus de développement dans la résolution du problème est l’algorithme
de couplage que nous avons détaillé dans la section 4.1.4. Cet algorithme comporte une boucle
itérative permettant de trouver les variables de couplage (partition de l’interface et données
pour les conditions aux limites) en résolvant l’équation de l’onde cinématique par le schéma de
Godunov ainsi que l’équation de Richards par l’algorithme itératif 1. Le stockage et le tracé
de la charge hydraulique et de la hauteur d’eau au cours de la simulation permettent de suivre
l’évolution du calcul en temps réel et de détecter s’il y a un problème (mauvaise CI ou CL, pas de
temps δt trop grand) sans attendre la fin de la simulation (cela permet un gain de temps pendant
le débuggage). La fréquence de tracé des diverses grandeurs (charge hydraulique, hauteur d’eau,
isovaleurs, nuages de points, champs de vitesse) est choisie par l’utilisateur. Le post-traitement
de nos données de sortie permettant d’obtenir notamment le bilan de masse s’effectue avec
c c
Matlab
mais pourrait aussi être réalisé avec Scilab
.
Ces éléments montrent qu’il est facile d’enrichir les fonctionnalités du logiciel en ajou-
tant des modèles pour décrire les propriétés hydrodynamiques du sol, en ajoutant des for-
mules d’intégration et en changeant éventuellement le logiciel de maillage et le logiciel de post-
traitement des résultats.
B.2. Assemblage de la matrice et du second membre 107
Nous avons implémenté une conductivité hydraulique tensorielle puisque les termes K11 , K12 ,
K21 et K22 apparaissent dans le listing ci-après. Ayant remarqué que le temps CPU augmente en
prenant en compte tous les termes du tenseur, nous avons introduit une condition (via le booléen
tenseur extradiag) pour calculer l’apport des termes extradiagonaux seulement si l’utilisateur
le souhaite. Cette condition n’a pas été utilisée dans les cas tests présentés dans cette thèse
mais l’anisotropie est un élément disponible dans le logiciel R+R DG. Des simulations avec des
sols hétérogènes pourraient également être réalisées puisque les quatre composantes du tenseur
dépendent de l’espace.
if ( terme_source )
{ b_K(et) += pds_2d(j) * fonc_base_2d[j] * source(x_pts_2d(et)(j),y_pts_2d(et)(j),tps) ; }
M_K(et)[i1] += pds_2d(j)
* ( fonc_dx_2d(j,i1) * fonc_dx_2d[j] * K11(psi,x_pts_2d(et)(j),y_pts_2d(et)(j)) +
fonc_dy_2d(j,i1) * fonc_dy_2d[j] * K22(psi,x_pts_2d(et)(j),y_pts_2d(et)(j)) ) ;
if ( tenseur_extradiag ) {
b_K(et) -= pds_2d(j) * K12(psi,x_pts_2d(et)(j),y_pts_2d(et)(j)) * val_fonc_dy_2d[j] ;
M_K(et)[i1] += pds_2d(j) *
( fonc_dx_2d(j,i1) * fonc_dy_2d[j] * K21(psi,x_pts_2d(et)(j),y_pts_2d(et)(j)) +
fonc_dy_2d(j,i1) * fonc_dx_2d[j] * K12(psi,x_pts_2d(et)(j),y_pts_2d(et)(j)) ) ; }
}
}
b(et) *= jacobien(et) ;
M_K(et) = M_K(et) * jacobien(et) ;
}
108 Annexe B. Le logiciel R+R DG
La contribution dans le second membre provenant d’un élément est composée d’une intégrale
comportant la conductivité hydraulique,
Z
bK = − K(ψh)ez ∇φ.
τ
Nous rappelons que cette intégrale fait suite à l’intégration par parties effectuée pour éviter
d’introduire la dérivée de la conductivité hydraulique (voir la section 2.2.5). Pour vérifier nos
programmes lors des tests de convergence avec des solutions analytiques (voir la section 2.2.2),
nous avons pris en compte un terme source volumique dépendant de l’espace et du temps (intitulé
source dans le listing et activé via le booléen terme source),
Z
sf = f φ.
τ
Dans le listing concernant la boucle sur les éléments, le tableau contenant les valeurs des fonctions
de base aux points d’intégration numérique est noté fonc base 2d et permet de calculer la
matrice M C. Les tableaux val fonc dx 2d et val fonc dy 2d contiennent les valeurs des dérivées
par rapport à x et y des fonctions de base aux points d’intégration numérique. En utilisant la
transformation linéaire, notée inv T K, qui permet de passer de l’élément courant et à l’élément
de référence, ces deux tableaux permettent de calculer la matrice M K.
def 1 + def
{ξ}F = (ξ + ξ − ) et [[ξ]]F = ξ − − ξ + ,
2
avec ξ + = ξ|τ + et ξ − = ξ|τ − . Pour une fonction vectorielle, les opérateurs de moyenne et de saut
sont définis composante par composante. De plus, nF désigne la normale unitaire à F dirigée de
τ − vers τ + . Les intégrales effectuées sur une face interne contiennent les opérateurs de moyenne
et de saut et font donc intervenir quatre matrices élémentaires A IP, B IP, C IP et D IP qui
traduisent respectivement la contribution de ψ + sur τ + , de ψ − sur τ − , de ψ + sur τ − et de ψ −
sur τ + .
Les différentes contributions dans la matrice provenant des intégrales contenant le flux sca-
laire Fψ sont
Z Z Z
1 1
K(ψh+ )∇φ · nτ + ({ψh} − ψh+ ) = K(ψh+ )∇φ · nτ + ψh− − K(ψh+ )∇φ · nτ + ψh+ ,
F 2 F 2 F
Z Z Z
1 1
K(ψh− )∇φ · nτ − ({ψh} − ψh− ) = K(ψh− )∇φ · nτ − ψh+ − K(ψh− )∇φ · nτ − ψh− .
F 2 F 2 F
Ces quatre intégrales sont prises en compte dans les différentes matrices élémentaires par
l’intermédiaire des vecteurs vec 1, vec 2, vec 3 et vec 4 calculés à partir des vecteurs
élémentaires a elem, b elem, c elem et d elem et des composantes de la normale Norm(·,·)(1)
et Norm(·,·)(2).
B.2. Assemblage de la matrice et du second membre 109
Les différentes contributions dans la matrice provenant des intégrales contenant la première
partie du flux Fu sont
Z Z Z
1 + + 1
{K(ψh)∇ψh}F · nτ + φ = K(ψh )∇ψh · nτ + φ + K(ψh− )∇ψh− · nτ + φ,
F 2 F 2 F
Z Z Z
1 1
{K(ψh)∇ψh}F · nτ − φ = K(ψh− )∇ψh− · nτ − φ + K(ψh+ )∇ψh+ · nτ − φ.
F 2 F 2 F
Notons que les deux intégrales qui traduisent la contribution de ψ + sur τ + ainsi que la contri-
bution de ψ − sur τ − sont directement calculées à partir des calculs précédents (en effectuant les
transposées des matrices A IP et B IP).
Les différentes contributions dans la matrice provenant des intégrales contenant la seconde
partie du flux Fu sont
Z Z Z
+
[[ψh]]F nF · nτ + φ = ψh φ − ψh− φ,
F F F
Z Z Z
[[ψh]]F nF · nτ − φ = ψh− φ − ψh+ φ.
F F F
Ces quatre intégrales sont ainsi calculées à partir de la matrice de masse élémentaire et muti-
pliées par des matrices de permutations. Elles sont notées tPk MM 1d Pk(·,·) dans le listing et
sont multipliées par le coefficient ηKs d−1
F . Six permutations (en deux dimensions d’espace) sont
possibles pour passer de la face courante à la face de référence définie sur l’élément de référence.
Les six matrices de permutations Pk(·) que l’on trouve dans les trois boucles effectuées sur les
faces sont rectangulaires, dépendent du degré d’approximation polynômial (taille 2×3 en P1 ,
3×6 en P2 ) et sont calculées en pré-traitement. Ces matrices permettent de sélectionner correc-
tement les degrés de liberté sollicités sur une face (suivant le numéro de permutation) parmi
l’ensemble des degrés de liberté du triangle. Pour éviter de calculer à chaque pas de temps et
lors des deux boucles itératives (boucle de l’algorithme de couplage et boucle de l’algorithme de
linéarisation) les matrices faisant intervenir ces matrices de permutations, nous avons également
déterminé en pré-traitement les trente-six matrices tPk MM 1d Pk(·,·) possibles.
Deux intégrales proviennent de l’intégration par parties réalisées dans le second membre pour
éviter d’introduire la dérivée de la conductivité hydraulique :
Z Z
+ −
I1 = K(ψh|τ + )ez φ nτ + et I1 = K(ψh|τ − )ez φ nτ − .
F F
Cette boucle portant sur les faces internes fait intervenir les gradients des fonctions de base
qu’il faut déterminer en fonction du numéro de la permutation de la face F suivant que l’on
se place sur le triangle τ + (fonc dx 1d et fonc dy 1d) ou sur le triangle τ − (fonc dx 1d v et
fonc dy 1d v). Dans ce listing, nous retrouvons également le booléen tenseur extradiag qui
permet de prendre en compte la contribution éventuelle des termes extra-diagonaux du tenseur
de conductivité.
110 Annexe B. Le logiciel R+R DG
if ( tenseur_extradiag ) {
a_elem += fonc_dy_1d_et[j] * K12(psi,pts_1d(1),pts_1d(2)) ;
b_elem += fonc_dx_1d_et[j] * K21(psi,pts_1d(1),pts_1d(2)) ;
c_elem += fonc_dy_1d_ev[j] * K12(psi_v,pts_1d_v(1),pts_1d_v(2)) ;
d_elem += fonc_dx_1d_ev[j] * K21(psi_v,pts_1d_v(1),pts_1d_v(2)) ;
A_IP += transpose(A_IP) ;
B_IP += transpose(B_IP) ;
La contribution dans le second membre provenant d’une face ayant une condition de Dirichlet
est composée de trois intégrales I1 , I2 et I3 ,
Z Z Z
I1 = K(ψh)∇φ · nΩ ψD , I2 = φψD et I3 = K(ψh)ez φnτ .
F F F
Nous rappelons que les intégrales I1 et I2 font intervenir les flux numériques et contiennent la
valeur de Dirichlet ψD que nous souhaitons imposer sur l’ensemble des faces FhD . Nous observons
dans le listing ci-dessous que cette donnée ψD est une fonction dépendant de l’espace et du temps.
Comme pour les faces internes, l’intégrale I3 provient de l’intégration par parties réalisée dans
le second membre pour éviter d’introduire la dérivée de la conductivité hydraulique.
if ( tenseur_extradiag ) {
I_3 += pds_1d(j) * fonc_base_1d[j] * K12(psi,x_pts_1d(ea)(j),y_pts_1d(ea)(j)) * Norm(ea,1)(2) ; }
}
E_IP += transpose(E_IP) ;
E_IP = ARETE[ea-1].longueur() * (-E_IP + eta / h_K(et) * tPk_MM_1d_Pk(i_et,i_et) ) ;
M_K(et) += E_IP ;
b(et) += ARETE[ea-1].longueur() * ( -I_1 + eta / d_K(et) * tPk(i_et) * I_2 + tPk(i_et) * I_3 ) ;
}
112 Annexe B. Le logiciel R+R DG
Le listing ci-dessous montre que la donnée dans la condition aux limites de Neumann vN dépend
de l’espace et du temps. Cette fonctionnalité est utilisée dans les deux cas tests où la pluie s’arrête
au cours de la simulation (cas test avec drain de la section 3.2.2 et cas test de ruissellement dû
à la pluie dans la section 4.2.2) ainsi que dans le cas test avec l’injection d’eau dans le sol de la
section 4.2.3.
A chaque pas de temps, les faces sont réparties en fonction de leur label dans les quatre
tableaux suivants
– arete dirichlet continu pour le label 1,
– arete dirichlet discret pour le label 2,
– arete neumann continu pour le label 3,
– arete neumann discret pour le label 4.
Nous précisons que l’initialisation de la première valeur de chaque tableau à NA+1 (NA étant le
nombre total de faces) est utile dans la suite des programmes pour savoir si tous les types de
conditions aux limites sont utilisés. Nous avons également implémenté un algorithme de tri pour
classer les deux groupes de faces situées sur la surface suivant le sens des x croissants : la fonction
tri a bulle est appliquée aux tableaux arete dirichlet discret et arete neumann discret
lorsqu’il existe au moins une face ayant le label 2 et 4 respectivement. Cet algorithme de tri,
dédié au cas monodimensionnel de l’onde cinématique, permet une écriture plus compacte du
flux de Godunov.
switch ( ARETE[ea-1].condition() )
{ case 1 : arete_dirichlet_continu(i_dirichlet_continu) = ea ;
i_dirichlet_continu += 1 ;
break ;
case 2 : arete_dirichlet_discret(i_dirichlet_discret) = ea ;
i_dirichlet_discret += 1 ;
break ;
case 3 : arete_neumann_continu(i_neumann_continu) = ea ;
i_neumann_continu += 1 ;
break ;
case 4 : arete_neumann_discret(i_neumann_discret) = ea ;
i_neumann_discret += 1 ;
break ;
default : cout << "Probleme dans les conditions aux limites " << endl ;
cout << "Programme Conditions_limites.h " << endl ;
exit(1) ;
}
}
nb_dirichlet_continu = i_dirichlet_continu -1 ;
nb_dirichlet_discret = i_dirichlet_discret -1 ;
nb_neumann_continu = i_neumann_continu -1 ;
nb_neumann_discret = i_neumann_discret -1 ;
if ( arete_dirichlet_discret(1) != NA+1 )
{ r = tri_a_bulle( arete_dirichlet_discret ) ; }
if ( arete_neumann_discret(1) != NA+1 )
{ r = tri_a_bulle( arete_neumann_discret ) ; }
114 Annexe B. Le logiciel R+R DG
if ( nb_dirichlet_discret != 0 )
{ vecteur <int> numero_dirichlet(nb_dirichlet_discret) ; }
if ( nb_neumann_discret != 0 )
{ vecteur <int> numero_neumann(nb_neumann_discret) ; }
int i1 = 1 , i2 = 1 ;
switch ( ARETE[ea-1].condition() )
{ case 2 : numero_dirichlet(i1) = i ;
i1 += 1 ;
break ;
case 4 : numero_neumann(i2) = i ;
i2 += 1 ;
break ;
[ABB+ 04] Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and
Benoı̂t Perthame. A fast and stable well-balanced scheme with hydrostatic recons-
truction for shallow water flows. SIAM J. Sci. Comput., 25(6) :2050–2065 (electro-
nic), 2004.
[ABCM01] D.N. Arnold, F. Brezzi, B. Cockburn, and L.D. Marini. Unified analysis of disconti-
nuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 39(5) :1749–
1779 (electronic), 2001.
[AD02] V. Aizinger and C. Dawson. A discontinuous Galerkin method for two-dimensional
flow and transport in shallow water. Advances in Water Resources, 25 :67–84, 2002.
[AG84] A.S. Abdul and R.W. Gillham. Laboratory studies of the effects of the capillary
fringe on streamflow generation. Water Resources Research, 20(6) :691–698, 1984.
[Arn82] D.N. Arnold. An interior penalty finite element method with discontinuous elements.
SIAM J. Numer. Anal., 19(4) :742–760, 1982.
[Bak77] G. A. Baker. Finite element methods for elliptic equations using nonconforming
elements. Math. Comp., 31(137) :45–59, 1977.
[Bas99] P. Bastian. Numerical computation of multiphase flows in porous media. Habilitation
Thesis, 1999.
[Bas03] P. Bastian. Higher order discontinuous Galerkin methods for flow and transport
in porous media. In Challenges in scientific computing—CISC 2002, volume 35 of
Lect. Notes Comput. Sci. Eng., pages 1–22. Springer, Berlin, 2003.
[BC64] R.H. Brooks and A.T. Corey. Hydraulic properties of porous media. Fort Collins,
Colorado : Colorado State University. Hydrology paper, 3 :27, 1964.
[BEE+ 06] H. Beaugendre, A. Ern, T. Esclaffer, E. Gaume, I. Ginzburg, and C.Kao. A seepage
face model for the interaction of shallow water tables with the ground water surface :
Application of the obstacle-type method. Journal of Hydrology, 329 :258–273, 2006.
[BIR+ 07] P. Bastian, O. Ippisch, F. Rezanezhad, H.J. Vogel, and K. Roth. Numerical simula-
tion and experimental studies of unsaturated water flow in heterogeneous systems.
In Reactive Flows, Diffusion and Transport, pages 579–597. Springer Berlin Heidel-
berg, 2007.
[BJ67] G. Beavers and D. Joseph. Boundary conditions at a naturally impermable wall. J.
Fluid. Mech., 30 :197–207, 1967.
[BK04] M. Bause and P. Knabner. Computation of variably saturated subsurface flow
by adaptive mixed hybrid finite element methods. Advances in Water Resources,
27 :565–581, 2004.
[BN82] M.D. Barcelo and J.L. Nieber. Influence of soil pipe networks on catchment hydro-
logy, 1982. paper presented at the Summer meeting, Am. Soc. Agric. Eng., Univ of
Wisc., Madison.
118 Bibliographie
[BR04] P. Bastian and B. Rivière. Discontinuous galerkin methods for two-phase flow in
porous media. Technical Report 2004–28, IWR (SFB 359), Universität Heidelberg,
2004.
[BT78] Z.P. Bazant and W. Thonguthai. Pore pressure and drying of concrete at high
temperature. J. Eng. Mechanics Division. ASCE, 104 :1059–1079, 1978.
[CB90] M.A. Celia and E.T. Bouloutas. A general mass-conservative numerical solution for
the unsaturated flow equation. Water Resources Research, 26(7) :1483–1496, 1990.
[Cia02] P.G. Ciarlet. The Finite Element Method for Elliptic Problems. SIAM, 2002.
[CJ] E. Cuthill and J.McKee. Reducing the bandwidth of sparse symmetric matrices.
Naval Ship Research and Development Center. Washington, D.C. 20007.
[CKG07] J.P. Carlier, C. Kao, and I. Ginzburg. Field-scale modeling of subsurface tile-drained
soils using an equivalent-medium approach. Journal of Hydrology, 341 :105–115,
2007.
[Cor54] A.T. Corey. The interrelation between gas and oil relative permeabilities. Producers
Monthly, 19(1) :38–41, 1954.
[CP88] R.F. Carsel and R.S. Parrish. Developing joint probability distributions of soil water
retention characteristics. Water Resources Research, 24(5) :755–769, 1988.
[CS98a] B. Cockburn and C.-W. Shu. The local discontinuous Galerkin method for time-
dependent convection-diffusion systems. SIAM Journal on Numerical Analysis,
35(6) :2440–2463 (electronic), 1998.
[CS98b] B. Cockburn and C-.W. Shu. The runge-kutta discontinuous galerkin finite element
method for conservation laws v : Multidimensional systems. J. Comput. Phys.,
141 :199–224, 1998.
[Dar56] H. Darcy. Les fontaines publiques de la ville de Dijon. V. Dalmont, Paris, 1856.
[Daw06] C. Dawson. Analysis of discontinuous finite element methods for ground wa-
ter/surface water coupling. SIAM Journal on Numerical Analysis, 44(4) :1375–1404,
2006.
[DB70] T. Dunne and R.D. Black. Partial area contributions to storm runoff in a small new
england watershed. Water resources research, 7 :1160–1172, 1970.
[DB02] L. Dormieux and E. Bourgeois. Introduction à la micromécanique des milieux poreux.
Presses de l’École Nationale des Ponts et Chaussées, 2002.
[DLA01] E.B. Diaw, F. Lehmann, and Ph. Ackerer. One-dimensional simulation of solute
transfert in saturated-unsaturated porous media using the discontinuous finite ele-
ment method. Journal of Contaminant Hydrology, 51 :197–213, 2001.
[DMQ02] M. Discacciati, E. Miglio, and A. Quarteroni. Mathematical and numerical models
for coupling surface and groundwater flows. Applied Numerical Mathematics, 43 :57–
74, 2002.
[DP99] H.-J.G Diersch and P. Perrochet. On the primary variable switching technique for
simulating unsaturated-saturated flows. Advances in Water Resources, 23 :271–301,
1999.
[DPEG08] D.A Di Pietro, A Ern, and J.L Guermond. Discontinuous Galerkin methods for
anisotropic semidefinite diffusion with advection. SIAM J. Numer. Anal., 46(2) :805–
831, 2008.
[dSV71] A.J.C de Saint-Venant. Théorie du mouvement non-permanent des eaux, avec ap-
plication aux crues des rivières et à l’introduction des marées dans leur lit. Compte-
Rendu à l’Académie des Sciences de Paris, 73 :147–154, 1871.
Bibliographie 119
[EFGV00] M. Esteves, X. Faucher, S. Galle, and M. Vauclin. Overland flow and infiltration
modelling for small plots during unsteady rain : numerical results versus observed
values. Journal of hydrology, 228 :265–282, 2000.
[EG06] A. Ern and J.L. Guermond. Discontinuous Galerkin methods for Friedrichs’ sys-
tems. II. Second-order elliptic PDEs. SIAM J. Numer. Anal., 44(6) :2363–2388
(electronic), 2006.
[Ein88] B. Einfeldt. On godunov-type methods for gas dynamics. Journal of Computational
Physics, 25 :294–318, 1988.
[ENV07] A. Ern, S. Nicaise, and M. Vohralı́k. An accurate H(div) flux reconstruction for
discontinuous Galerkin approximations of elliptic problems. C. R. Math. Acad. Sci.
Paris, 345(12) :709–712, 2007.
[EPD08] A. Ern, S. Piperno, and K. Djadel. A well-balanced Runge-Kutta discontinuous
Galerkin method for the shallow-water equations with flooding and drying. Internat.
J. Numer. Methods Fluids, 58(1) :1–25, 2008.
[Esc04] T. Esclaffer. Etude théorique de la formation des débits de crues à l’échelle du bassin
versant. Technical report, Ecole Nationale des Ponts et Chaussées, 2004.
[ESZ08] A. Ern, A.F. Stephansen, and P. Zunino. A discontinuous Galerkin method with
weighted averages for advection-diffusion equations with locally small and anisotro-
pic diffusivity. IMA J. Numer. Anal., 2008. to appear.
[FFRH04] S. Fagherazzi, D.J. Furbish, P. Rasetarinera, and M. Youssuff Hussaini. Application
of the discontinuous spectral Galerkin method to groundwater flow. Advances in
Water Resources, 27 :129–140, 2004.
[FK59] I. Fatt and W.A. Klikoff. Effect of fractional wettability on multiphase flow through
porous media. Trans., AIME (Am. Inst. Min. Metall. Eng.), 216 :426–432, 1959.
[Fre] Logiciel Freefem. http ://www.freefem.org/.
[FS86] G. Fipps and R.W. Skaggs. Drains as a boundary condition in finite elements. Water
resources research, 22(11) :1613–1621, 1986.
[GA11] W.H. Green and G. Ampt. Studies on soils physics,1. The flow of air and water
through soils. J. Agric. Sci., 4(1) :1–24, 1911.
[GCK04] Irina Ginzburg, Jean-Philippe Carlier, and Cyril Kao. Lattice boltzmann approach
to richards equation. In Proceedings of Computational Methods in Water Resources
2004 International Conference, 2004.
[Gen80] M.T. Van Genuchten. A closed-form equation for predicting the hydraulic conduc-
tivity of unsaturated soils. Soil Sci Soc Am J, 44 :892–898, 1980.
[Gin04] V. Ginting. Computational upscaled modeling of heterogeneous porous media flow
utilizing finite volume method. PhD thesis, Texas AM University, 2004.
[GKJR96] R.S. Govindaraju, M.L. Kavvas, S.E. Jones, and D.E. Rolston. Use of green–ampt
model for analyzing one-dimensional convective transport in unsaturated soils. Jour-
nal of hydrology, 178 :337–350, 1996.
[GL96] G.H. Golub and C.F. Van Loan. Matrix computations. The Johns Hopkins University
Press, London, 1996.
[GMMS02] M. Gabbouhy, A. Maslouhi, Z. Mghazli, and Z. Saâdi. Modélisation numérique du
transport de soluté dans la zone non saturée d’un sol trés sableux. Math-Recherche
et Applications, 4 :77–99, 2002.
120 Bibliographie
[GP01] J.-F. Gerbeau and B. Perthame. Derivation of viscous saint-venant system for la-
minar shallow water ; numerical validation. Discrete and Continuous Dynamical
Systems : Series B, 1(1) :89–102, 2001.
[GSK78] R. Gharaaty-Sani and L.G. King. Interceptor drains on sloping land, 1978. paper
presented at the Summer Meeting, Am. Soc. Agric. Eng., Logan, Utah.
[GY75] A.B. Gureghian and E.G. Youngs. The calculation of steady-state water-table
heights in drained soils by means of the finite element method. Journal of hydrology,
27 :15–32, 1975.
[Haz11] A. Hazen. Discussion : dams on sand fondations. Transactions of the American
Society of Civil Engineers, 73 :199–203, 1911.
[Hel97] R. Helmig. Multiphase Flow and Transport Processes in the Subsurface : A Contri-
bution to the Modelling of Hydrosystems. Springer-Verlag, 1997.
[Her03] J.M. Hervouet. Hydrodynamique des écoulements à surface libre : Modélisation
numérique avec la méthode des éléments finis. Presses de l’École Nationale des
Ponts et Chaussées, 2003.
[HG93] S.M Hassanizadeh and W.G. Gray. Thermodynamic basis of capillary pressure in
porous media. Water Resour. Res., 29 :3389–3405, 1993.
[HLvL83] A. Harten, P.D. Lax, and B. van Leer. On upstream differencing and godunov-type
schemes for hyperbolic conservation laws. SIAM Rev., 25 :35–61, 1983.
[Hor33] R.E. Horton. The role of infiltration in the hydrologic cycle. Trans. Am. Geophys.
Union, 14 :446–460, 1933.
[HR06] X. He and L. Ren. A multiscale finite element linearization scheme for the unsa-
turated flow problems in heterogeneous porous media. Water Resources Research,
42 :W08417, 2006.
[HVT+ 77] R. Haverkamp, M. Vauclin, J. Touma, P.J. Wierenga, and G. Vachaud. A compa-
rison of numerical simulation models for one-dimensionnal infiltration. Soil Science
Society of America Journal, 41 :285–294, 1977.
[JM00] W. Jäger and A. Mikelic. On the interface boundary condition of Beavers, Joseph
and Saffmann. SIAM J. Appl. Math., 60 :1111–1127, 2000.
[JW00] J E Jones and C S Woodward. Newton-krylovmultigrid solvers for large-scale, highly
heterogeneous, variably saturated flow problems. manuscript in preparation. In
Advances in Water Resources, pages 763–774, 2000.
[KKB+ 02] K. R. Kavanagh, C. T. Kelley, R. C. Berger, J. P. Hallberg, and S. E. Howing-
ton. Nonsmooth nonlinearities and temporal integration of richards equation. In
Proceedings of the XIV International Conference on Computational Methods in Wa-
ter Resources, S. Majid Hassanizadeh, Ruud J. Schotting, W. G. Gray, and G. F.
Pinder, editors, pp 947-954., 2002.
[KM06] S.J. Kollet and R.M. Maxwell. Integrated surface-groundwater flow modeling : A
free-surface overland flow boundary condition in a groundwater flow model. Ad-
vances in Water Resources, 29 :945–958, 2006.
[KN08] P.M. Kienzler and F. Naef. Subsurface storm flow formation at different hillslopes
and implications for the ’old water paradox’. Hydrological processes, 22 :104–116,
2008.
[Koz27] J. Kozeny. About capillaries conducting water in the earth. Commitee Report of the
Viennese Academy, 136(2a) :271–306, 1927.
Bibliographie 121
[Sto45] G. Stokes. On the theories of the internal friction of fluids motion and of the
equilibrium and motion of elastic solids. Trans. Cambridge Phil. Soc., 8 :287–305,
1845.
[TBV07] P.A. Tassi, O. Bokhove, and C.A. Vionnet. Space discontinuous Galerkin method
for shallow water flows-kinetic and HLLC flux, and potential vorticity generation.
Advances in Water Resources, 30 :998–1015, 2007.
[Td06] P. Tardif d’Hamonville. Modélisation et simulation du transport advectif et diffusif
an milieu poreux monophasique et diphasique. PhD thesis, Ecole Nationale des Ponts
et Chaussées, 2006.
[TdED07] P. Tardif d’Hamonville, A. Ern, and L. Dormieux. Finite element evaluation of diffu-
sion and dispersion tensors in periodic porous media with advection. Computational
Geosciences, 11(1) :43–58, 2007.
[Tor01] E.F. Toro. Shock-capturing methods for free-surface shallow flows, Wiley, New York,
2001.
[VCEL03] P.L. Viollet, J.P. Chabard, P. Esposito, and D. Laurence. Mécanique des fluides
appliquée. Presses de l’École Nationale des Ponts et Chaussées, 2003.
[VGC01] T. Vogel, M.Th. Van Genuchten, and M. Cislerova. Effect of the shape of the soil
hydraulic functions near saturation on variably-saturated flow predictions. Advances
in Water Resources, 24 :133–144, 2001.
[VL01] J.E. VanderKwaak and K. Loague. Hydrologic-response simulations for the R-5
catchment with a comprehensive physics-based model. Water Resources Research,
37(4) :999–1013, 2001.
[VT62] B.S. Vimoke and G.S. Taylor. Simulating water flow in soil with an electric resistance
network. Soil and Water Conserv. Res. Div., U.S. Agirc. Res. Serv.,Columbus,
Ohio., 32 :1–136, 1962. Rep. 41-65, 51pp.
[WD00] C.S. Woodward and C. Dawson. Analysis of expanded mixed finite element methods
for a nonlinear parabolic equation modeling flow into variably saturated porous
media. SIAM J. Numer. Anal., 37(3) :701–724 (electronic), 2000.
[Whe78] M.F. Wheeler. An elliptic collocation-finite element method with interior penalties.
SIAM J. Numer. Anal., 15(1) :152–161, 1978.
.