Modelisation Des Ecoulements
Modelisation Des Ecoulements
THESE
présentée
pour obtenir
Résumé
Dans ce travail, on considère la modélisation des écoulements dans des hydro-systèmes
comprenant des sols et des aquifères géologiquement complexes et hétérogènes. On considèrera par
exemple le cas d’un aquifère côtier soumis à l’intrusion saline, avec couplage densitaire (eau douce /
eau salée), phénomène auquel peuvent se greffer d’autres couplages (écoulements à saturation
variable, couplages surface / souterrain). On choisit une approche ayant les caractéristiques
suivantes :
♦ le modèle est spatialement distribué afin de représenter l’hétérogénéité du milieu ;
♦ le modèle est fortement couplé afin d’appréhender les écoulements dans leur complexité
physique.
On utilise dans ce but un modèle fortement intégré, à une seule équation générique de type EDP,
basée sur une loi de Darcy généralisée permettant de décrire différents « régimes » d’écoulements la
co-existant dans un même domaine, tout en conservant robustesse et efficacité.
Le travail est divisé en trois parties.
Dans une première partie on élabore un nouveau modèle numérique 3D, pour la modélisation des
écoulements en milieux poreux à densité variable dans l’hypothèse d’une interface abrupte. Ce
nouveau modèle est basé sur des relations ‘effectives’ non linéaires de saturation et de perméabilité,
dans une équation d’écoulement de type Richards modifiée.
La seconde partie correspond à l’élaboration et l’implémentation d’un modèle verticalement intégré
d’intrusion saline en aquifère côtier, permettant d’étudier l’effet de l’hétérogénéité stochastique de
l’aquifère. Le modèle, basé sur l’hypothèse ‘interface abrupte’, est implémenté comme un module
‘2D’ dans le code volumes finis BIGFLOW 2D/3D. Le nouveau module 2D est utilisé pour analyser
la variabilité de l’interface eau douce / eau salée par simulations stochastiques de type Monte Carlo à
échantillonnage spatial (réalisation unique). Ces résultats sont comparés à nouvelle théorie, où
l’interface aléatoire auto-corrélée est analysée par transformation de variable, combinée à une
méthode de perturbation et à une représentation spectrale (Fourier / Wiener-Khinchine).
Dans la troisième et dernière partie, on présente un modèle de couplage fortement « intégré » pour la
modélisation des écoulements de surface et souterrain en hypothèses d’écoulement plan,
verticalement hydrostatique. On s‘intéresse au cas d’une vallée fluviale avec cours d’eau, plaine
d’inondation, et nappe d’accompagnement. L’écoulement en surface est modélisé par l’équation
d’onde diffusante et l’écoulement souterrain par l’équation de Dupuit-Boussinesq. Ce modèle couplé
est appliqué à la vallée fluviale de la Garonne dans la région de Toulouse - Moissac (France). Cette
application a nécessité l’élaboration d’une méthode d’interpolation géostatistique adaptée à
l’élaboration d’un Modèle Intégré Numérique de Terrain (‘MINT’), de façon à inclure le fond de la
rivière au MNT topographique en haute résolution. Enfin, au-delà de cette application particulière, le
modèle d’écoulement couplé surface / souterrain est généralisé au cas d’un couplage densitaire eau
douce / eau salée, lorsque la nappe est sujette à l’intrusion saline au voisinage d’une embouchure ou
d’un estuaire.
Mots Clés
Intrusion saline Modelisation stochastique Couplage surface/souterrain
Milieu poreux et hydrogéologie Aquifère et nappe souterraine Volumes finis 2D / 3D
Modèle Intégré Numérique de Terrain (MINT, MNT) Loi de Darcy Ward
Equation de Richards Equation de Dupuit-Boussinesq Equation d’onde diffusante
ii
Abstract
In this work, we consider water flow modeling in hydro-systems that include geologically
complex and heterogeneous soils and aquifers, e.g., a coastal aquifer undergoing seawater intrusion,
with density coupling (freshwater / saltwater), along with other coupled phenomena (variable
saturation, surface / subsurface coupling). The selected approach has the following characteristics:
♦ the model is spatially distributed in order to represent the heterogeneity of the medium
♦ the model is strongly coupled in order to apprehend the physical complexity of flow systems
We use for this purpose a strongly integrated model, governed by a single generic equation (PDE)
based on generalized Darcy law, to describe different flow ‘regimes’ co-existing in the same
domain, while conserving robustness and efficiency.
The work is divided into three parts:
In the first part, we develop a new 3D numerical model for variable density flow in porous media
under the sharp interface approximation. This new model is based on non-linear ‘effective’
saturation and conductivity relations, in a modified Richards flow equation.
The second part corresponds to the development and implementation of a vertically integrated
saltwater intrusion model, to study the effect of stochastic heterogeneity in a coastal aquifer.
The model, based on the sharp interface hypothesis, is implemented as a 2D module in the finite
volumes code BIGFLOW 2D/3D. The new module is used for analyzing the variability of the
salt / fresh interface through Monte Carlo simulations with spatial sampling (single realization).
These results are compared to a new theory where the random field interface is analyzed via a
transformation combined to a perturbation method and a spectral representation (Fourier / Wiener-
Khinchine).
In the third and last part, we present a strongly integrated model to simulate coupled surface /
subsurface plane flows, such as a river valley with stream, floodplain, and free surface aquifer.
Surface flow is modeled via the diffusive wave equation, and subsurface flow is modeled using the
Dupuit-Boussinesq equation. This coupled model is applied to the Garonne river valley in the
Toulouse-Moissac region (France). This application has required the elaboration of a geostatistical
interpolation technique that produces an Integrated Digital Elevation Model (‘IDEM’). The IDEM
incorporates a high resolution representation of river channels into the topographic DEM. Finally,
beyond this specific application, the coupled surface / subsurface model is generalized to the case of
salt / fresh density coupling, where the aquifer is subject to saltwater intrusion near a river mouth or
an estuary.
Keywords
Saltwater intrusion Stochastic modeling Surface / subsurface coupling
Porous media & hydrogeology Groundwater & aquifers Finite volumes 2D / 3D
Integrated Digital Elevation Model (IDEM, DEM) Darcy-Ward law
Richards equation Dupuit-Boussinesq equation Diffusive wave equation
iii
Remerciements
Je tiens tout d’abord à remercier en premier lieu mon directeur de thèse, Rachid Ababou, pour son
soutien scientifique. Ses conseils promulgués durant ma thèse m’ont permis d’élargir mes horizons.
Je n’oublie pas aussi sa chaleur humaine qui m’a beaucoup touché.
Je remercie aussi les membres du jury devant lesquels j’ai eu l’honneur de présenter ma thèse.
Messieurs Philippe Renard, Philippe Ackerer, Claudio Paniconi ont accepté de juger ce travail en
tant que rapporteurs, et je les en remercie. Ils ont contribué par leurs nombreuses remarques et
suggestions à améliorer la qualité de ce mémoire, et je leur en suis très reconnaissant.
Mes remerciements vont également à Monsieur Michel Quintard pour m’avoir accueilli au sein de
l’équipe GEMP de l’IMFT, et pour m’avoir fait découvrir les karsts des gorges du Tarn ! Je n’oublie
pas les membres de l’équipe GEMP et les doctorants de l’IMFT.
Merci à l'équipe LEH Laboratoire d’Ecologie des Hydrosystèmes du CNRS pour les données du
site de Monbéqui.
Ce travail n’aurait pas pu aboutir sans la bourse de thèse du Ministère de l’Education et sans le
soutien financier du projet européen SWIMED.
Mes remerciements les plus affectueux vont à ma femme Hanaa, mes parents et mes trois sœurs
qui m’ont toujours soutenu dans mes choix.
iv
Chapitre I INTRODUCTION
Introduction 2
I – 1 Contexte _______________________________________________________________ 3
I - 1.1 Les défis ______________________________________________________________ 3
I - 1.2 Modélisation des transferts ______________________________________________ 4
I - 1.3 Hétérogénéité, milieux aléatoires, et géostatistique ___________________________ 4
I - 1.4 Modélisation stochastique _______________________________________________ 5
I - 1.5 Types d’écoulements, types de milieux, et lois phénoménologiques ______________ 6
I – 2 Objectifs de la thèse ______________________________________________________ 8
I – 3 Plan de la thèse__________________________________________________________ 9
Introduction 3
I-1 CONTEXTE
I - 1.1 Les défis
L’eau douce constitue 2.5 % de l’eau dans le globe, et 95% de l’eau douce utilisable est dans le
sous sol. Ces chiffres montrent l’importance de cette source d’eau. Elle est d’autant plus importante
si l’on considère la pénurie d’eau à laquelle le monde devra faire face dans le futur. Le nombre de
pays touchés par la surexploitation (à plus de 40%) des eaux souterraines va significativement
augmenter d’ici 2025 comme on le voit dans la figure 1.1. Les pays riches et moins riches devront
faire face à ce problème. Ces chiffres globaux sont encore plus alarmants lorsqu’on considère la
répartition inégale des ressources en eaux dans un même pays et à travers le temps, et si l’on
considère la concentration démographique dans les grandes villes et sur les côtes. Cette
surexploitation combinée aux changements climatiques majeurs conduira éventuellement à une
raréfaction et une dégradation de la qualité des ressources en eaux.
Fig. 1.1 Prévision des contraintes sur l’eau douce (WMO, 96).
Les défis majeurs à prendre en compte en gestion des ressources en eaux sont : la déplétion des
aquifères à cause de la surexploitation, l’abaissement du niveau de certaines rivières à cause des
changements climatiques, l’assèchement des lacs, la salinification des sols à cause de la
surexploitation agricole, la pollution des aquifères par les pesticides, les engrais et autres substances
chimiques nocives, et l’intrusion saline. Cette dernière sera accrue par la surélévation du niveau de la
mer et par la surexploitation des nappes côtières.
Ces différents problèmes font intervenir des hydro-systèmes complexes formés de
« composantes » divers comme les rivières, les lacs, les eaux souterraines et la végétation. Les eaux
souterraines peuvent elles-mêmes être constituées par un assemblage de différentes sous-
composantes telles que : nappes phréatiques d’eau douce, nappes confinées plus ou moins
profondes, zones côtières soumises à l’intrusion saline, sols non saturés, etc.
Ces composantes interagissent entre elles. Un abaissement du niveau des rivières conduira à un
rabattement du niveau des nappes en plaine alluviale. Une surexploitation des nappes conduira à un
abaissement du niveau des rivières. Ce type d’interaction dans les deux sens nécessite la
considération d’un couplage fort entre les deux composantes. Ce couplage est encore plus complexe
dans le cas d’un aquifère côtier connecté à la mer. Dans ce cas l’eau douce moins dense glisse en
Introduction 4
dessus de l’eau de mer plus dense qui envahit l’aquifère sous la forme d’un biseau salé. Ces deux
volumes d’eau interagissent aussi au niveau de la zone de transition qui se forme entre l’eau douce et
l’eau salée.
Plusieurs travaux de recherches ont été consacrés à l’élaboration d’une approche globale de
modélisation des écoulements ; citons parmi les plus récents : Putti et Paniconi (2004), Panday et
Huyakorn (2004), Gunduz et Aral (2005), Kollet et Maxwell (2006). D’autres projets ambitieux sont
en cours actuellement, comme le projet SEVE de couplage de codes hydrologiques. Ces approches,
combinées à des études de scénarios basées sur des règles d’évolutions socio-économiques,
permettront aux décideurs de gérer les risques de façon globale sur l’ensemble de l’écosystème afin
de prendre les meilleures décisions avec la plus faible incertitude.
denses. En pratique, on dispose de mesures qui ne recouvrent pas de façon suffisamment détaillée la
topographie - ou bathymétrie - de la rivière. Dans beaucoup de cas les profils le long de la rivière
sont très éloignés les uns des autres (plusieurs centaines de mètres). De telles contraintes exigent des
méthodes d’interpolation adaptées à ce type de problème en vue de l’obtention d’un Modèle
Numérique de Terrain (MNT). Lorsqu’un tel MNT résulte de l’intégration ou de la fusion de
données variées (altimétries sur grilles régulières, altimétries ponctuelles, profils en travers, profils
en longs), on pourrait parler plutôt de Modèle Intégré et Numérique de Terrain (MINT).
ensemble ou sur la totalité des réalisations générées. Enfin une reconstitution des moments de la
variable d’intérêt est effectuée à partir des résultats des simulations numériques. Plusieurs
milliers de réalisations sont nécessaires pour obtenir de bonnes estimations statistiques sur les
moments.
• Ecoulement saturé :
Ecoulement en milieu poreux saturé avec une porosité effective totalement remplie d’eau.
• Ecoulement non-saturé :
Ecoulement en milieu non saturé, la porosité étant seulement partiellement remplie d’eau et le
reste d’air. La résistance à l’écoulement de l’air est négligée et la pression de l’air est en
équilibre avec la pression de l’atmosphère en permanence. Les équations classiques régissant
l’écoulement non saturé (sans saturation variable) sont les équations de Richards (Richards
1931).
• Ecoulement variablement saturé :
L’eau s’écoule dans un milieu partiellement saturé, i.e., contenant des zones saturées et d’autres
zones non saturées. Dans le cas général, plusieurs zones saturées et non saturées peuvent co-
exister, et leurs distributions spatiales peut évoluer dans le temps, e.g.: écoulement partiellement
saturé avec front infiltration descendant vers la surface libre d’une nappe (Freeze 1971 ; Vauclin
et al. 1979). Afin de modéliser l’écoulement en milieu partiellement saturé dans un seul domaine
et avec une seule équation, l’équation de Richards est reformulée en variables mixtes et sous
forme conservative, avec deux variables d’état, la teneur en eau et la pression : voir Ababou et
al. (1988, 1992) et Celia et al. (1990), parmi d’autres.
• La loi de Darcy, la perméabilité et la charge hydraulique :
La loi de Darcy exprime la proportionnalité entre la densité de flux q [L T-1] à travers un milieu
poreux et le gradient de charge hydraulique dans le milieu poreux. Notons que la vitesse de
l’eau est donnée par V = q/θ, où θ est la porosité (ou la teneur en eau). La perméabilité [m2],
intrinsèque au milieu poreux, exprime l’inverse de la résistance visqueuse au flux. La
conductivité hydraulique [m/s] exprime la même chose, mais elle dépend du fluide (eau). Dans
la zone non saturée, la loi de Darcy est de la même forme –quasi linéaire – avec une
conductivité hydraulique fonction de la pression ou la teneur en eau. Parmi les références
historiques liées à ce paragraphe nous citons par exemple Darcy (1856); Buckingham (1907);
Richards (1931).
• Milieux macroporeux ‘dynamiques’ et ‘cinétiques’ :
Un milieu macroporeux peut être considéré comme un milieux ouvert, une couche superficielle
de sol couvert de végétation dense, ou un banc de galet, etc. Dans le code de calcul BIGFLOW
par exemple, ce type de milieu est représenté comme un milieu poreux à pores grossiers et à
Introduction 7
forte perméabilité, allant à l’extrême vers une perméabilité quasi-infinie et une porosité de
100% dans certains sous domaines (Ababou et al. 1996, 1998, 2002, 2006 ; Trégarot 2000).
Comme expliqué dans les références précitées, cette approche originale conduit à prendre en
compte explicitement deux effets distincts pour les écoulements en macroporeux à saturation
variable :
a) Effets dynamiques : à cause du fort nombre de Reynolds, la perte de charge n’est plus
linéaire, les effets inertiels deviennent importants, et la loi de Darcy est remplacée par la loi
quadratique en vitesse de Ergun-Ward-Forchheimer (e.g.,Ward 1964) ;
b) Effets cinétiques : la courbe de rétention teneur en eau / pression est une fonction escalier
(ou quasi-escalier, vue la discrétisation numérique), conduisant à une imbibition ou drainage
instantané de la porosité en tout endroit où une surface libre existe.
• Ecoulements plans de Dupuit-Boussinesq :
En considérant un écoulement verticalement hydrostatique quasi plan, l’équation de masse et la
loi de Darcy peuvent être intégrées verticalement pour obtenir l’équation de Dupuit-Boussinesq
en écoulement plan (x,y). Les variables d’état deviennent : la charge hydraulique moyenne H
[m], et le flux spécifique en 2D Q [L2 T-1]. La seconde variable provient de l’intégration de la
vitesse de Darcy. Dans ce cas, la zone non saturée est négligée (réponse instantanée de la
surface libre au processus d’imbibition et de drainage). Le coefficient de drainage est la porosité
effective. Cette approche se base sur les travaux de Dupuit (1863) et Boussinesq (1904).
• Ecoulements plans de Saint-Venant :
Les équations 2D de Saint-Venant (Saint Venant, 1871) sont une approximation verticalement
hydrostatique et verticalement intégrées de Navier-Stokes (équations de conservation de masse
et de quantité de mouvement) dans le cas d’écoulement hydraulique à surface libre en canaux,
rivières, plaines d’inondation, etc. Les variables considérées sont le vecteur de vitesse de l’eau
en 2D, V [L1 T-1], et la hauteur d’eau “H” ou la lame d’eau “h”. La formulation 1D des
équations de Saint-Venant est largement utilisée en hydrologie de surface pour décrire les
écoulements transversalement intégrés le long des rivières et des canaux.
• Ecoulements plans d’Onde Diffusante et couplage nappe-rivière :
L’équation 2D d’onde diffusante (cinématique et diffusive) est une simplification
supplémentaire de l’équation de Saint-Venant (et donc de Navier-Stokes), valide seulement
pour des écoulements à vitesse lentement variée et suffisamment faible pour que les termes
inertiels soient négligeables. Il existe une version 1D pour les écoulements en canal, tandis que
la version 2D est utilisable pour décrire les ruissellements et inondations de surfaces en (x,y) –
ce qui est d’intérêt dans notre travail sur le couplage surface/souterrain. En effet, l’idée est de
coupler un écoulement de surface avec un écoulement souterrain peu profond en utilisant l’onde
diffusante pour l’eau de surface et l’équation de Boussinesq pour l’eau souterraine. Voir Chow
et al. (1988) ou Bedient et al. (2002) pour une présentation « basique » du modèle d’onde
diffusante ; ces auteurs présentent aussi des tables de valeurs des coefficients de rugosités
correspondants.
• Anisotropie :
Les milieux poreux et macroporeux naturels peuvent être fortement anisotropes. Cependant, la
définition de l’anisotropie dépend aussi de la résolution spatiale du modèle. Par exemple,
l’anisotropie d’un milieu stratifié n’est pas importante si le modèle d’écoulement est implémenté
avec une résolution plus fine que l’épaisseur des strates. La loi de Darcy peut être formulée avec
une perméabilité anisotrope avec un tenseur symétrique du second ordre.
Dans le modèle numérique BIGFLOW, le tenseur est supposé diagonal, mais une rotation peut
être appliquée au domaine de calcul par rapport au système de référence « horizontal/vertical ».
Le modèle équationnel de BIGFLOW prend en compte l’anisotropie dans plusieurs formes de la
Introduction 8
loi de Darcy : ainsi dans le modèle 2D plan, la transmissivité et le coefficient de rugosité sont
anisotropes ; pour les écoulements 3D macroporeux à grand Reynolds, la loi quadratique en
vitesse (Ward-Forchheimer) peut aussi être anisotrope, comme indiqué dans Trégarot (2000),
d’après Knupp et Lage (1995).
• Ecoulements à densité variable en milieux poreux (et hypothèse d’interface abrupte) :
Un écoulement en milieu poreux est à densité variable lorsque la densité dépend de la position et
du temps, notamment à travers d’autres variables d’état (telles que concentration, température,
etc.). On s’intéresse ici au cas de l’intrusion saline en nappe côtière, où la concentration de l’eau
en sel affecte l’écoulement. L’eau est le fluide porteur (solvant), et le sel est le soluté. Ce
problème pourrait être décrit par deux équations : une équation d’écoulement d’eau et une
équation de transport advectif / diffusif de sel. Cependant, on peut mettre en oeuvre une version
simplifiée en utilisant l’hypothèse « interface abrupte ».
• Ecoulements à densité variable en milieux poreux avec hypothèse d’interface abrupte :
Dans le modèle d’interface abrupte, on considère un écoulement à densité variable défini par
deux zones fluides de densités distinctes séparées par une interface abrupte. En effet, à grande
échelle et pour de faibles contrastes de densité (aquifères côtiers), et étant donné la variabilité
temporelle des forçages hydrologiques (marées, précipitations, etc.), l’épaisseur de la zone de
transition due à la diffusion pure de sel entre l’eau de mer et l’eau douce est relativement petite
par rapport à la taille du domaine, et par rapport à la variabilité importante due à l’hétérogénéité
intrinsèque du milieu. C’est cette approche « interface abrupte » qui est utilisée dans le cadre de
cette thèse, l’une verticalement intégrée, et l’autre, relativement novatrice, tridimensionnelle.
ainsi qu’avec les équations plus générales des milieux macroporeux à perte de charge quadratique
(Ward-Forchheimer). Puis une procédure de couplage surface/souterrain est développée. Le modèle
couplé est ensuite appliqué, en premier lieu, à des cas simplifiés (méandres rectangulaires, etc.).
Enfin la méthode est appliquée à l’écoulement dans une partie de la vallée alluviale de la Garonne au
niveau de Monbéqui, entre les villes de Toulouse et Moissac. Pour faire cette simulation, un modèle
numérique intégrant le fond de la rivière et la topographie était indispensable. Ainsi une méthode
pour l’obtention d’un modèle intégré numérique de terrain (MINT) à partir d’un nombre limité de
sections de rivières est également élaborée et appliqué au site « Garonne ».
Dans le chapitre 6 les conclusions majeures, les résultats et les perspectives découlant de ce
travail sont présentés.
Modélisation des écoulements à densité variable 11
II - 1 NOTIONS DE BASE
Dans les paragraphes suivants les quelques notions de base qui interviennent dans la description
des écoulements en milieux poreux seront présentées.
II - 1.2 Définitions
• Porosité φ : rapport du volume des vides au volume total Vt du sol
Vvide
φ= (2.1)
Vt
• Teneur en eau volumique θ (m3 /m3) : rapport du volume de l'eau Veau contenue dans les pores
(ou vides) du sol, au volume total Vt du sol considéré :
Veau
0 <θ = <1 (2.2)
Vt
• Saturation S : rapport du volume d'eau au volume des vides dans le sol :
S=θ φ (2.3)
• La conductivité hydraulique et la perméabilité
La conductivité à saturation peut s'exprimer sous la forme d'un produit de deux facteurs, l'un
dépendant des caractéristiques du fluide (l'eau), l'autre dépendant uniquement des caractéristiques de
la phase solide et traduisant la facilité avec laquelle le sol se laisse traverser par les fluides, appelée
perméabilité intrinsèque k (m2 ) :
ρg g
Ks = k= k (2.4)
μ ν
où ν et μ sont respectivement les viscosités cinématique (m2 /s) et dynamique (kg/m/s) de l'eau.
Modélisation des écoulements à densité variable 14
eau adsorbé
air
eau
particule
Fig. 2.1 Représentation d'un milieu poreux
− 2 σ cos(ϕ ) − A
h= = (2.7)
ρgr r
où
σ: Tension superficielle de l'eau, inversement proportionnelle à la température
(σ ≈ 0.0725 Nm-2 à 20 C) ;
ϕ: Angle de contact entre l'interface eau / air et la paroi du capillaire, supposé nul dans
les sols (en supposant une mouillabilité parfaite de la phase solide vis à vis de l'eau) ;
ρ: Masse volumique de l'eau (en négligeant celle de l'air) ;
r: Rayon du pore assimilé à un tube capillaire, égal au rayon du ménisque si ϕ = 0.
En résumé, en zone non saturée, l'eau est en sous pression par rapport à la pression atmosphérique
(h < 0). Une pression très fortement négative h → -∞ (ou forte succion) signifie une prédominance
des forces d'adsorption face aux forces de capillarité, les rayons de courbures des interfaces eau / air
étant petits. Une pression faiblement négative (faible succion) ou proche de zéro (pression proche de
la pression atmosphérique) signifie au contraire que les rayons des interfaces eau/air sont grands,
et/ou que les pores sont presque saturés.
A l’échelle macroscopique (supérieure au VER), l'interface entre la zone saturée et la zone non
saturée est définie, en terme de pression, par h = 0, c’est-à-dire, p = pATM.
L’iso-surface h(x,y,z,t) = 0 représente, lorsqu’elle existe, une surface libre interne au milieu poreux
(cas d’un domaine poreux « partiellement saturé »).
Pour dériver le terme de droite de l'équation (2.9), on considère la matrice solide comme rigide et
immobile mais faiblement déformable dans la direction verticale, z. On considère aussi que la
densité varie en fonction de la pression et de la concentration en soluté (Ackerer et al. 1999). On
néglige la variation de la densité en fonction de la température. Ainsi on obtient :
∂Mass ⎡ ∂p θ ∂ρ ∂p θ ∂ρ ∂C θ ∂Δz ∂p ⎤
=ρ ⎢θ ∂t + ρ ∂p ∂t + ρ ∂C ∂t + Δz ∂p ∂t ⎥ Δx Δy Δz (2.11)
∂t ⎣ ⎦
où p la pression de l'eau en [F/L2] et C est la fraction massique du soluté.
On définit α le coefficient de compressibilité élastique du milieux poreux en [L2/F], et β le
coefficient de compressibilité élastique de l'eau [L2/F] (Bear 1972).
1 ∂Δz 1 ∂ρ
α= et β = (2.12) et (2.13)
Δz ∂p ρ ∂p
En remplaçant (2.12) et (2.13) dans (2.11) on obtient :
∂Mass ⎡ ∂p θ ∂ρ ∂C ∂ θ ⎤
∂t
=ρ ⎢θ (α + β ) ∂t + ρ ∂C ∂t + ∂t ⎥ Δx Δy Δz (2.14)
⎣ ⎦
où :
qD
Re = (2.18)
υ
avec q la vitesse [L/T]
υ la viscosité cinématique (μ/ρ)
D la longueur représentatif de l’écoulement [L]
Pour les écoulements en milieux poreux, q est considéré comme la densité de flux, et D est égale
au diamètre effectif des grains d10 (la taille du filtre à travers duquel 10% en masse des grains passe).
Afin que la loi de Darcy soit applicable il faut que Re ≤ 1. Muskat (1937) montre plusieurs études
selon lesquelles l'écoulement est transitoire entre des valeurs de Reynolds comprises entre 1 et 12. A
titre indicatif pour un sable de diamètre moyen de 0.5 mm la vitesse d'écoulement limite pour l'eau à
température ambiante est de 0.2 cm/s. Cette valeur peut être dépassée dans le cas d'écoulement en
milieux fracturés ou karstiques. Dans ce cas une loi de perte de charge quadratique (Darcy-
Forchheimer) est utilisé (voir Chapitre 5).
ou encore :
∂θ (h, x )
= ∇.[K ( h, x ) ∇h] − ∇.[K ( h, x ) g B ] + Qs (2.20)
∂t
On remarque que le terme de stockage élastique 'M' est négligé en comparaison au stockage dû à
la variation de la teneur en eau (Freeze et Cherry, 1979). Le terme M peut être considéré dans les
zones saturées au cours d'une simulation d'un milieu poreux variablement saturé, surtout lorsque la
zone saturée est sujette à de fortes variations de pression et dans les zones argileuses sujettes à des
phénomènes de gonflement et de retrait. Le terme C de variation de la densité en fonction de la
concentration est nul car on considère que la densité est constante.
L'éq. (2.20) aux Dérivées Partielles (EDP) est une forme conservative mixte en h/θ de l'équation
de Richards exprimée à l'origine en h. Elle est de type parabolique en milieu non saturé (et non
linéaire) et elliptique en milieu saturé (si le terme de sportivité spécifique traduisant les effets de
compressibilité est nul). Un terme M non nul préserve le caractère parabolique de l'équation, évitant
ainsi certaines difficultés de convergence qui se présentent lorsque l'équation devient elliptique dans
les zones saturées (Trégarot 2002). Elle est applicable à tout instant t et en tout point x à l'intérieur
du domaine de calcul, la loi de Darcy étant utilisée pour traiter les frontières du domaine. La
résolution de la forme mixte de l'équation de Richards (2.20) permet donc de simuler les
écoulements variablement saturés (h), en milieux hétérogènes (x) et anisotropes.
θe = θ − θd (2.21)
1 Hystérèse. Dans la nature, un sol est soumis à des successions de phases de drainage et d'imbibition plus
ou moins complètes, suivant des courbes θ(h) non uniques, comprenant des courbes primaires et des
courbes de passage imbibition/drainage. Ce phénomène d'hystérèse peut s’expliquer en partie par les
variations de l'angle de contact ϕ formé par les interfaces eau / air et eau / solide. Ainsi, à une pression h
donnée, θ est plus importante lors d'une phase de drainage du sol que lors d'une phase d'imbibition
Cependant, dans la pratique de la modélisation hydrologique à l'échelle d'un versant, ce phénomène
d'hystérèse est négligé devant les effets des hétérogénéités, et une relation θ(h) univoque est considérée
pour définir localement les propriétés du sol.
Modélisation des écoulements à densité variable 19
S − Sd θ − θd θe
Se = = = (2.22)
S s − S d θ s − θ d θ es
Le Tableau 2.1 présente une liste des modèles de θ(h) les plus utilisés:
Tab. 2.1 Modèles de saturation effective (humidité) Se (h), d'après (Trégarot 2000).
Auteur(s) Modèle θ(h) Remarques
Le modèle K(Se) de Mualem (1976) associé au modèle Se(h) de Van Genuchten (1980) donne:
K (h) =
[1 − ( −αh ) [1 + ( −αh ) ]
n −1 n −m
]
2
[1 + ( −αh ) ]
r
(2.27)
n m/2
Le même modèle de K(Se) (Mualem 1976) associé cette fois-ci au modèle Se(h) de Brooks et Corey
(1964) donne:
2 + 2 ,5 / b
⎛h ⎞
K r (h) = ⎜ b ⎟ (2.28)
⎝ h ⎠
L'équation (2.29) est linéaire et parabolique (elliptique si le terme de storativité spécifique Ss = 0).
La loi de comportement de Darcy exprime le débit spécifique Qs (en m3/s/m), ou bien la densité
de flux q (en m3/s/m2 ), comme suit :
Qs = ηq = − K s η ∇Zs (2.30)
L'équation de conservation de masse s’écrit :
∂Z s
φe = −div[Q s ] = −div[ηq] (2.31)
∂t
avec :
T= ∫K
Zinf
s dz , avec cette fois, Ks = Ks (x,y,z) (2.35)
Lorsque les variations temporelles de la surface libre Zs sont négligeables par rapport à la valeur
moyenne de l'épaisseur η = Zs - Zinf , ou lorsque la répartition verticale de Ks est telle qu'elle entraîne
de faibles variations temporelles de T, alors T = T(x,y) et l'équation (2) devient linéaire.
Modélisation des écoulements à densité variable 23
Dans le cas général d’une nappe confinée en aquifère hétérogène et à plancher et toit variables,
nous avons : S = S(x,y), H = H(x,y,t), Ks = Ks (x,y), Zsup = Zsup (x,y), Zinf = Zinf (x,y).
L'équation (3) est de type parabolique, linéaire en raison de la transmissivité constante
ρ = f (Ci, p, T ) (2.38)
Cette interdépendance lorsqu'elle ne peut pas être négligée, nous mène à considérer le couplage à
densité variable. Ces problèmes appliqués aux milieux poreux sont rencontrés dans différents
systèmes naturels et industriels en hydrogéologie (ou hydrologie souterraine), en géophysique, en
génie de réservoir, en génie nucléaire, et en génie de matière. Les applications varient du transport
des polluants denses, à l'intrusion saline en aquifères côtiers, à l'infiltration de lixiviat dans les
décharges et le stockage des déchets industriels, à la conception des systèmes de chauffage
géothermique, à la convection dans les couches de neige, et tant d'autres (Dierch et Kolditz, 2002).
La même analyse peut être faite sur la viscosité cinématique qui varie en fonction de la
concentration, mais son effet est négligeable dans les problèmes qui nous intéressent (II - 3.2.2).
L'un des problèmes majeur en hydrogéologie sujet à l'écoulement à densité variable est l'intrusion
saline. En fait l’intrusion de l’eau de mer dans les aquifères d’eau douce est un phénomène naturel
qui se produit dans les zones côtières. Dans une configuration classique d’intrusion saline l’eau
douce glisse sur l’eau salée plus lourde dont la densité est supérieure à 1022 kg/m3. L’eau salée
forme ainsi un biseau dans l’aquifère d’eau douce. A ce phénomène vient s’ajouter l’effet du
pompage dans les aquifères côtiers. Dans cette configuration une zone salée de forme conique ce
forme au niveau du puits (Diersch et al. 1984, Diersch et Nillert 1990, Reilly et Goodman 1987 et
Holzbecher 1995).
Deux familles de modèles dont chacune intègres plusieurs variantes sont utilisées pour la
modélisation des écoulements à densité variable appliquée à l'intrusion saline :
La première repose sur les équations couplées du transport et d'écoulement des fluides miscibles.
Parmi les travaux de recherche qui ont utilisé cette méthode on peut citer : Segol et al. 1975,
Huyakorn et al. 1987, Frind 1982, Voss (1984,1999), Voss et Souza (1987), Putti et Paniconi (1995),
Diersch (1988), Kolditz et al. (1998). Cette méthode a vu un fort développement avec l'augmentation
de la performance et de la capacité des calculateurs depuis une 20 année. Diersh et Kolditz (2002)
présentent un état de l’art de cette méthode.
La seconde considère deux fluides non-miscibles avec une interface abrupte les séparant.
Plusieurs formulations existent pour cette approche. Certaines sont basées sur les modèles de suivis
d’interface. D’autres ressemblent aux écoulements multiphasiques en milieux poreux. Dans le cas où
la zone salée est quasi statique, l'hypothèse de Badon-Ghyben (1888) et de Herzberg (1901) est
appliquée. Nous utiliserons cette approche pour développer notre modèle 3D d'intrusion saline. Les
deux approches seront présentées dans les deux sections suivantes.
C Advection
C Advection et diffusion
X
Advection, diffusion
et dispersion
C
∂ ( ρθ )
C = −C∇.( ρq ) + Cρ S QS (2.42)
∂t
Enfin on soustrait l'équation (2.42) de l'équation (2.41) :
∂C
θρ = − ρφ V ⋅ ∇ C + ∇ ⋅ ( ρφ D ⋅ ∇ C ) + ρQ (C s − C ) (2.43)
∂t
Les coefficients de dispersion du tenseur de diffusion sont obtenus par une des méthodes suivantes :
• Model géométrique (Taylor et Aris) :
Dans ce modèle, les coefficients sont obtenus pour une configuration simple, déterministe. Par
exemple, dans un milieu constitué de cylindres, le coefficient de dispersion Dii est donné par :
o écoulement parallèle aux cylindres :
D xx
= 0 .002 Pe 2 (2.44)
Dm
o écoulement perpendiculaire aux cylindres :
D xx
= 0 .07 Pe 1.7 (2.45)
Dm
où Pe est le nombre adimensionnel de Peclet qui exprime le rapport entre la convection
forcée (advection) et la diffusion moléculaire. Le nombre de Peclet est donné par :
V L
Pe = (2.46)
Dm
où V est la vitesse moyenne d’écoulement [L/T] et L une longueur caractéristique [L]
• Méthode stochastique de changement d’échelle (ou Method of Volume Averaging) :
Dans les méthodes de changement d’échelle, les propriétés macroscopiques sont calculées pour
plusieurs configurations d’hétérogénéité à partir de changement d’échelle (Dagan 1982, Gelhar et
Axness 1983). Les résultats obtenus prennent la forme du modèle de puissance :
D xx
= a Pe b (2.47)
Dm
tel que 1 < b < 2. Pour des cylindres aléatoirement organisés en couches : a=0.7 et b=1.2
• Modèle statistique (Bear 1961, Scheiddegger 1961) :
Le tenseur de dispersion est donné par :
VV
D = D m I + (α L − α T ) + αT V I (2.48)
V
avec Dm la diffusion moléculaire effective après prise en compte de la tortuosité en [L2/T], I le
tenseur unité et αL et αT les coefficients de dispersivité intrinsèque, longitudinale et transversale
en [L]. Le coefficient αL est défini dans la direction principale de l'écoulement. Il varie de
plusieurs ordre suivant le degré d'hétérogénéité du domaine et la longueur de l'écoulement
(Gelhar 1982). Le coefficient αT est défini suivant la direction transversale à la direction
principale de l'écoulement, Il est plus petit. En pratique on prend αT = 0.1 à 0.01 αL.
Les composantes du tenseur de diffusion regroupant la diffusion moléculaire et la dispersion
obtenues par la méthode statistique sont données par (Bear 1979):
Modélisation des écoulements à densité variable 27
2 2 2
1 vx 1 vy 1 vz
Dxx = α L + αT + αT + Dm τ xx
θ v θ v θ v
2 2 2
1 vy 1 vx 1 vz
D yy = α L + αT + αT + Dm τ yy
θ v θ v θ v
2 2 2
1 vz 1 vx 1 vy
Dzz = α L + αT + αT + Dm τ zz
θ v θ v θ v (2.49)
1vv
Dxy = D yx = (α L − α T ) x y
θ v
Dxz = Dzx = (α L − α T )
1 vx vz
θ v
1 v y vz
D yz = Dzy = (α L − α T )
θ v
Dm représente la diffusion moléculaire et ( τ xx , τ yy , τ zz ) sont les composantes principales du
tenseur de tortuosité. Ce modèle a été largement étudié dans les travaux de Lever et Jackson
(1985), Hassanizadeh (1986), Kolditz et al. (1997) parmi d'autre.
⎡ ρ 0 ε ∂C ∂θ ⎤ 1
∂p
⎢θ S p ∂t + θ ρ C + = ∇[ρ K (∇p − ρ g )] + ρ s
∂t ⎥⎦ ρ
QS (2.54)
⎣ max ∂t ρ
En considérant la charge hydraulique de l’eau douce comme variable d'état, on obtient enfin :
⎡ ∂h f ρ 0 ε ∂C ∂θ ( h f ) ⎤ 1 ⎡ 2 ⎛⎜ ρ − ρf ⎞⎤ ρ s
⎟ ⎥ + Qs
θ
⎢ s ∂t + θ + = ∇ ⋅ ⎢ ρ ∇ + ∇
∂t ⎥⎦ ρ
S K ⎜ fh z ⎟ (2.55)
⎣ ρ Cmax ∂t ⎢⎣ ⎝ ρf ⎠⎥⎦ ρ
♦ Effet rétroactif de la viscosité
La viscosité cinématique varie en fonction de la concentration en sel et de la température. Dans
notre cas, on considère la viscosité à température ambiante. Les modèles qui donnent la variation de
la viscosité en fonction de la concentration « c » sont de nature empiriques. Herbert et al. (1988)
proposent le modèle suivant :
μ = 1 .002 × 10 −3 (1 + 0 .4819 c − 0 .2774 c 2 + 0 .7814 c 3 ) (2.56)
La figure (2.4) montre l’évolution de la viscosité cinématique en fonction de la concentration à
température ambiante. Cette variation est faible pour des faibles valeurs de concentration. Ainsi elle
est négligée dans la plus part des situations.
2,00E-03
1,90E-03
1,80E-03
1,70E-03
viscosité (m /s)
2
1,60E-03
1,50E-03
1,40E-03
1,30E-03
1,20E-03
1,10E-03
1,00E-03
0 0,2 0,4 0,6 0,8 1
concentration
II - 3.3.1 Introduction
Dans l’approche « interface abrupte », on considère l’eau salée et l’eau douce comme deux fluides
non-miscibles séparés par une interface : voir figure (2.5).
L’interface eau salée / eau douce est une surface imperméable en équilibre de pression. Autrement
dit, la pression est continue de part et d’autre de l’interface. Par contre, la densité du fluide est
discontinue de part et d’autre de l’interface. Il s’agit donc d’un modèle à deux fluides non miscibles,
sans diffusion de sel.
Modélisation des écoulements à densité variable 30
te p
zone
abru
eau salée
ce
zone
r fa
te eau douce
in
Deux méthodes sont possibles pour la résolution de l’intrusion saline avec une approche de type
interface abrupte (voir figure 2.6).
- La première est une méthode de suivi d’interface (surface tracking) qui consiste à diviser le
domaine en deux régions distinctes dont chacune est associée à une équation d’écoulement.
Ensuite, la solution des deux équations est utilisée pour retrouve la position de la surface qui
représente l’interface eau douce / eau salée (Bear, 1999). La méthode de suivi d’interface peut
être appliquée aussi dans le cas 2D ou quasi 3D.
- La deuxième méthode, dite ‘multi-phasique’, consiste à considérer l'eau salée et l'eau douce
comme deux fluides non-miscibles séparés par une interface, l’écoulement de chaque phase
fluide étant cependant résolu dans l ‘ensemble du domaine, de telle façon que la position de
l’interface est obtenue implicitement à la fin de la résolution du problème. Au niveau de
l'interface il n'y a pas de mélange (pas d’échange de masse) entre les deux phases fluides.
Théoriquement aucune diffusion n'est permise au niveau de l'interface abrupte.
Modélisation des écoulements à densité variable 31
Intégration
verticale
Formulation quasi 3D
(multicouche)
Formulation en 2D
Hypothèse
de G-H
Formulation en 2D Formulation en 3D Formulation en 3D
(G-H) (G-H) (G-H)
τ ( x, y , t ) = H S (1 + ε ) − H f ε (2.59)
ρf
où ε = représente le contraste de densité entre l’eau douce et l’eau salée.
ρs − ρ f
La condition de flux nul à travers l’interface pour la région d’eau douce est donnée par :
⎛qf ⎞
⎜⎜ − u ⎟⎟ • n = 0 (2.60)
⎝ φ ⎠
En insérant les différentes expressions dans l’équation précédente, on obtient :
∂H f ∂H s
φε − φ (1 + ε ) = [K sf ( x )∇H f ]⋅ [∇z − (1 + ε )∇H s + ε∇H f ] (2.61)
∂t ∂t
En appliquant le même calcul pour la zone d’eau salée, on obtient de même :
∂H f ∂H s
φε − φ (1 + ε ) = [K ss ( x )∇H s ]⋅ [∇z − (1 + ε )∇H s + ε∇H f ] (2.62)
∂t ∂t
pc = p1 − p2 = f ( Θ 2 )
(2.65)
avec 1 = phase non aqueuse ; 2 = phase aqueuse
o Contrainte de conservation du volume :
Cette contrainte relie la saturation des deux phases à la porosité effective du milieu. Elle
exprime la conservation de masse au sein d’un VER.
Θ1 Θ2
θ= + (2.66)
ρ1 ρ2
♦ Méthode de résolution
Les équations diphasiques peuvent être formulées en pression, en saturation, ou en formulation
mixte (pression pour la phase 1 et saturation pour la phase 2). Wu et Forsyth (2001) font une étude
comparative pour le choix de la meilleure variable appropriée pour les équations de Richards
(milieux non saturés), les équations diphasiques et tri-phasiques en milieux fortement hétérogènes.
Ils recommandent l’utilisation de la formulation en saturation pour l’équation de Richards, une
formulation mixte en pression-saturation pour le système diphasique, et une formulation saturation-
pression-saturation pour un système tri-phasique.
L’équation de Richards correspond bien à un problème diphasique avec deux phases fluides,
eau / air, mais avec une seule des deux phases dynamiquement active, l’eau.
Modélisation des écoulements à densité variable 34
ρ d g(h + H) = ρ s gH (2.67)
Cette relation est obtenue en imposant la relation de continuité de pression de part et d’autre de
l’interface douce/salée, et en appliquant les hypothèses quasi-hydrostatiques précédentes dans les
deux nappes (douce et salée). On suppose aussi que la pression atmosphérique est constante.
x L
h
mer eau douce
interface
H q
eau salée
plancher
Fig.2.7 Aquifère côtier soumis à l’intrusion saline (la mer est à droite).
On en déduit la relation entre la profondeur de l’eau salée (H) et l’épaisseur de la lentille douce (h) :
h
H= (2.68)
ε
ρs − ρ0
où ε est le contraste de densité [-] : ε = avec, comme valeur indicative ε ≈ 1 .
ρ0 40
La relation H = h/ε sera modifiée plus loin, dans le Chapitre IV, de façon à prendre en compte
une épaisseur verticale non nulle de l’exutoire de la nappe d’eau douce à la mer, ce qui est
évidemment plus réaliste qu’un exutoire de section infinitésimale. En effet, dans ce dernier cas,
l’exutoire est un point triple singulier triple ; la vitesse de sortie de l’eau douce à la mer y est infinie
même si le débit de sortie de l’eau douce reste fini.
Finalement, en appliquant l’hypothèse de Ghyben-Herzberg à l’un des modèles numériques
d’interface abrupte (modèle de suivi d’interface ou modèle diphasique), on voit que l’on peut réduire
le problème à deux équations en un problème à une seule équation...moyennant quand même
l’approximation assez restrictive d’une zone salée quasi-immobile.
Modélisation des écoulements à densité variable 35
Cependant, on verra que la zone salée peut être évolutive (quasi-équilibre évolutif), et aussi, que
l’hypothèse d’écoulements plans dans l’eau douce peut être levée si la continuité de pression à
l’interface (cf. Chapitre III) est appliquée de façon locale en 3D.
Critères OPTIONS
2La convergence au sens des grilles numériques signifie que la solution n’est pas améliorée
Tab 2.3 Exemple de classification des modèles hydrogéologiques utilisé dans SALTRANS.
Code AD EF S/U NC NF CV HT DIM SD
CODEBRIGHT * p S/U 3C 2/3 FE
d3f * p S/U 1C * 2/3 FV
FEFLOW * h S/U 1C * (*) * 2/3 FE
HYDRUS-2D * h S/U MC * 2 FE
MOCDENSE * p S 2C 2 AE
MOCDENS3D * h S 1C 2/3 AE
PSE2D * p S 1C 2 AE
RETRASO * p S/U MC * 2 FE
ROCKFLOW * p S/U 1C * 2/3 FE
SALTFLOW * h S 1C 2/3 FE
SEAWAT * h S MC 2/3 FD/AE
SUTRA/SUTRA3D * h/p S/U 1C * * 2/3 FE
TRANSDENSE * h S/U 1C * 2/3 FE
AD – Advection-Dispersion model
EF – Equation Formulation in terms of h = head, p = pressure
S/U – Saturated/Unsaturated
NC – Number of Components, MC = Multicomponent
NF – Non-Fickian dispersion
CV – Consistent Velocity approximation
HT – Heat Transport
DIM – Dimensions
SD – Spatial Discretization, FE = Finite Element, FD = Finite Difference, FV = Finite Volume, AE =
Analytical Element
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 37
III - 1 INTRODUCTION
Dans ce chapitre un modèle d’intrusion saline à interface abrupte est proposé. Ce modèle est basé
sur l’approche multiphasique à un fluide évoquée dans le chapitre 2. Contrairement à l’approche
multiphasique à deux fluides avec une pseudo capillarité proposée par Huyakorn (1996) seul
l’écoulement dans l’eau douce est modélisé. L’approche consiste à adapter le modèle d’écoulement
non-saturé basé sur les équations de Richards aux problèmes de l’intrusion saline. Dans ce modèle le
domaine est divisé en trois zones : la zone saturée en eau douce, la zone non-saturée, et la zone
saturée en eau salée. La division du domaine se fait à l’aide des courbes de rétention adaptées.
L’écoulement en zone salée est considéré comme quasi-statique. La zone salée est considérée
comme une zone dépourvue d’eau douce et qui s’adapte automatiquement aux variations de
l’écoulement dans la zone d’eau douce. L’adaptation se fait en appliquant la condition de flux nul et
de continuité de pression à l’interface entre l’eau douce et l’eau salée. Cette approche a été proposée
par Larabi et De Smedt (1997), Sbai (1999), et Aharmouch (2003) avec des courbes de retentions
qui présente un problème d’anti-diffusion. Dans notre modèle le problème de l’anti-diffusion qui
apparaît au niveau de la transition entre la zone eau douce / eau salée est traité en proposant des
courbes de retentions adaptées.
∂Θ α
+ ∇ ⋅ ∇qα = Sα avec α=1,2 (3.1)
∂t
Les deux fluides non miscibles (eau douce et eau salée) sont assimilés à deux phases. En plus la
troisième phase (air) est considérée de façon implicite dans chacune des deux phases. Cette
troisième phase est considérée comme non-active.
P1 = P 2 en Σ( x, t ) (3.3)
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 40
Zone
d’eau douce
P1
P2
Σ ( x, t)−
Σ ( x, t ) +
Zone
d’eau salée
Fig. 3.1 Représentation de l’interface abrupte entre la zone d’eau salée et la zone
d’eau douce.
A partir de cette définition on peut définir les pressions critiques pour chaque fluide pour lesquelles
la condition dynamique est vérifiée. La pression critique exprime un état statique ou l’eau douce et
l’eau salée sont en équilibre. Dans ce cas seule la diffusion moléculaire existe. Mais cette diffusion
faible par rapport à la taille du domaine est négligée.
• Pour le fluide 1 (eau douce) la valeur de la pression critique est donnée par :
ρ 2 p 2 − patm
h1crit = = (1 + ε ) ( H 2 ( x, y, z , t ) − z ) en Σ( x, t ) + (3.7)
ρ1 ρ 2 g
• Pour le fluide 2 (eau salée) on a :
ρ 1 p1 − patm
h2 crit = = ( H1 ( x, y, z, t ) − z ) /(1 + ε ) en Σ( x, t ) − (3.8)
ρ 2 ρ1 g
Si la pression dans un des fluides dépasse la pression critique, il peut potentiellement envahir l’autre
fluide. On peut ainsi définir le potentiel d’invasion. Le potentiel d’invasion de l’eau douce dans
l’eau salée en chaque point est donné par :
Le potentiel d’invasion de l’eau salée dans la zone d’eau douce est donné par :
Si en un point de l’interface, hpf ou hps ne sont pas nuls, alors ce point de l’interface n’est pas en
équilibre dynamique. Ainsi les potentiels d’invasion peuvent être utilisés comme paramètres dans les
courbes de rétention pour définir de façon implicite la transition entre la zone d’eau douce et d’eau
salée. La résolution itérative ou simultanée des deux équations avec les courbes de rétention ainsi
définies donnera la position finale de l’interface sans avoir recours à des méthodes sophistiquées de
suivi d’interface et d’imposition de condition limite interne avec remaillage.
Cette relation est identique à la formulation de Ghyben-Herzberg dans la zone salée mais elle est
différente dans la zone d’eau douce qui n’est pas verticalement intégrée. Aucune hypothèse n’est
formulée pour la zone d’eau douce, i.e., l’écoulement dans la zone d’eau douce est 3D non-
hydrostatique.
• la zone 2 d’eau douce définie par une pression positive plus faible que la pression de la
zone salée hydrostatique.
• la zone 3 d’eau salée définie par une pression supérieure à la pression hydrostatique
critique. Cette zone est saturée en eau salée. Dans ce modèle la zone salée est une zone
imperméable à l’eau douce donc qui ne contient pas de l’eau douce. L’eau salée n’est
présente que conceptuellement. Cette zone a une teneur nulle en eau douce.
θ − θd
= e β ( h−h ) b
θ β (h − hb )
pour h > hcrit =e pf
avec hpf=hcrit-h
θs
La courbe de conductivité pour le modèle de Van Genuchten Mualem est donnée par :
2
K − Kr ⎛ ⎡ ⎤ ⎞⎟
m
=
1 ⎜1 − 1 − 1
pour h < 0 Ks − Kr (
1 + (− α h )
n
)
−m / 2
⎝ ⎣
(
⎜ ⎢ 1 + (− α h )n ⎥ ⎟
⎦ ⎠)
pour 0 ≤ h ≤ hcrit K = K sat (3.17)
⎛ ⎡ ⎤
m
K − Kr 1 ⎜ 1
⎜1 − ⎢1 − ⎥
pour h > hcrit Ks − Kr
=
(
1 + (− α h pf ))
n − m '/ 2
⎝
(
⎜ ⎢⎣ 1 + (− α h pf )
n'
)
⎥
⎦
avec hpf=hcrit-h
K − Kr α ( h −h )
pour h > hcrit = e pf b avec hpf=hcrit-h
Ks − Kr
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 43
On voit dans les figures 3.2 et 3.3 les courbes de teneur en eau et de conductivité hydraulique
respectivement et pour les modèles Van Genuchten Mualem et exponentiel. Les différentes zones
(1,2 et 3) obtenues sont indiquées sur la figure 3.6.
0,40
0,35
h=0.0 hcri
0,25
θ(h)
0,20
0,10
0,05
0,00
-10 -5 0 5 10 15 20 25 30 35 40
h (m)
1,2
h=0.0 hcri
Mualem
1,0
Exponentielle
0,8
Zone 3
K(h)
Zone 2
0,6 Zone 1
0,4
0,2
0,0
-10 -5 0 5 10 15 20 25 30 35 40
h (m)
mn α (θ s − θ d ) ( − α h pf ) n −1
C (h) = −
pour h > hcrit [1 + ( −α h pf ) n ] m +1
avec hpf=hcrit-h
0,25
h=0.0 hcri
0,20
C(h) Mualem
0,15
0,10
0,05
C(h)=dθ/dh
0,00
-0,05
-0,10
-0,15
-0,20
-0,25
-10 -5 0 5 10 15 20 25 30 35 40
h
Fig. 3.4 Capacité obtenue pour le premier modèle de Van Genuchten Mualem
proposé.
Etant donné que la capacité est négative sur la transition de la zone 2 vers la zone 3 et de la zone 3
vers la zone 2 , la diffusivité sera aussi négative (voir figure 3.5). Ce résultat est fortement
contraignant car il implique des valeurs diagonales négatives dans la matrice du système linéaire.
Ainsi la matrice symétrique du système n’est plus définie positive et la solution de la matrice n’est
plus stable.
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 45
6,0E+02
D(θ) Mualem
4,0E+02
transition zone 2 --> zone 1
2,0E+02
D(θ)
0,0E+00
-2,0E+02
transition zone 2 --> zone 3
-4,0E+02
-6,0E+02
0 0,05 0,1 0,15 0,2 0,25 0,3 0,35
θ
Fig. 3.5 Diffusivité obtenue pour le premier modèle de Van Genuchten Mualem
proposé.
Etant donné que hcrit dépend de la position et plus spécifiquement de l’élévation, il existe plusieurs
courbes de rétention. La figure 3.6 présente les positions des courbes obtenues par rapport à la
configuration de l’aquifère. La figure 3.7 montre les courbes de θ(h) pour plusieurs positions en
élévation ou profondeur dans un aquifère de 30 m de profondeur
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 46
Mer
c Zone 2 : zone saturée en eau douce
d
Zone 3 : zone saturée en eau
Fig. 3.6 Localisation des positions pour lesquelles les courbes de rétention, de
conductivité et diffusivité sont montrées dans les figures 3.7, 3.8 et 3.9
respectivement.
(a) (b)
(c) (d)
(a) (b)
(c) (d)
La diffusion théorique obtenue à partir des courbes précédentes est montrée dans la figure 3.9.
On remarque que la diffusion ne présente plus des valeurs négatives.
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 48
(a) (b)
(c) (d)
x L
h
mer eau douce
interface
H q
eau salée
plancher
2εqx
h( x ) = −
K (ε + 1)
(3.23) et (3.24)
2qx
H ( x) = −
Kε (ε + 1)
Longueur de pénétration du coin salé (L) :
On peut déterminer la longueur L du coin salé (L est la distance par rapport à la mer de
l’intersection entre l’interface salée et le plancher de l’aquifère).
− Kε (ε + 1)H 0 2
En x=L, H(L)=H0 ⇒ L =
2q
• Profils pour x > L (au-delà de la zone du coin salé, vers l’intérieur des terres) :
Si on veut déterminer la forme de la nappe pour x>L, il faut réécrire la loi de Darcy intégrée
verticalement comme suit :
∂h
q = −K (h + H 0 )
∂x
De plus on a la condition limite en x=L : h(L)=εH0
En intégrant, on obtient pour x > L :
2q
h( x) = H 0 ²(1 + ε )² − (x − L ) − H 0 (3.25)
K
Cette solution correspond à une certaine interprétation du problème de Glover (1964), cité par
exemple dans Cheng et Ouazar (2003).
Développement d’un modèle 3D d’intrusion saline avec interface abrupte et zone salée quasi-statique 50
teneur en eau
Solution analytique
BFSWIM 3D (contours)
L = 40.8
On observe que l’interface abrupte de notre modèle correspond assez bien à une courbe iso-
concentration médiane (C = 0.5 en adimensionnel). Cependant, le biseau salé sur le substratum
imperméable n’a pas la même forme pour le modèle interface abrupte (absence d’inflexion).
Cette forme est modifiée dans le problème de Henry lorsqu’on tient compte de la dispersion
hydrodynamique (en plus de la diffusion moléculaire) (Abarca et al., 2006). Or les résultats de la
figure 3.15 montrent que les nouvelles iso-concentrations obtenues de cette façon ont une forme plus
proche de notre modèle interface abrupte.
2m
q=0 & dc/dz = 0
q = cte - c = 0
1m BIGSWIM3D
c =1
q = 0 & dc/dz=0
III - 7 CONCLUSIONS
Dans ce chapitre nous avons adopté l’approche d’interface abrupte à deux fluides non miscibles
pour la modélisation de l’intrusion saline. Nous avons formulé les équations régissant le
comportement de l’interface abrupte séparant l’eau douce et l’eau salée. Nous avons développé une
méthode numérique, qui ne nécessite pas un suivi lagrangien d’interface ni un remaillage
dynamique.
Notre méthode présente un avantage par rapport à l’approche avec diffusion de sel par le fait que
seules les équations concernant l’écoulement d’eau douce, dans les zones saturée et non saturée, sont
résolues. Dans cette approche, l’eau salée est considérée comme étant quasi-hydrostatique (vitesses
quasi-nulles), mais cet état d’équilibre de l’eau salée peut évoluer en temps. La zone salée n’est donc
immobile, elle varie en temps.
L’écoulement en zone non saturée, est décrit à l’aide de courbes caractéristiques de rétention Θ(h)
et de conductivité hydraulique K(h), qui sont ré-interprétées comme des fonctions caractéristiques
‘effectives’ pour le système composé de l’eau douce, l’eau salée, et la zone non saturée. Nous avons
testé une version continues des fonctions Θ(h) et K(h) proposé par d’autre auteur. Cette version a
présenté des problèmes d’anti-diffusion. Ensuite on a proposé des nouvelles fonctions présentant un
Dirac au niveau de l’interface. Cette astuce numérique a résolu le problème de l’anti-diffusion.
Deux tests on été effectué pour comparé l’approche. Les tests effectués montrent que la technique
proposée est capable de déterminer avec précision la position de l’interface eau douce/eau salée. La
validation, en régime permanent, des résultats numériques avec ceux de la solution analytique de
Glover est concluante. Une deuxième comparaison est effectuée avec les résultats numériques de la
simulation du problème de Henry (1964) qui utilise l’approche diffusive. La comparaison montre
qualitativement que l’interface modélise par notre approche est proche de l’iso surface C=0.5.
Modélisation stochastique de l'intrusion saline en 2D plan 55
Chapitre IV MODELISATION
STOCHASTIQUE DE L'INTRUSION
SALINE EN 2D PLAN
Modélisation stochastique de l'intrusion saline en 2D plan 56
Résumé ____________________________________________________________________ 57
1. INTRODUCTION _______________________________________________________ 71
2. PROBLEM FORMULATION AND NUMERICAL MODEL ___________________ 71
3. CONTROL OF NUMERICAL ACCURACY AND ROBUSTNESS ______________ 72
3.1 Continuation method for highly heterogeneous media _______________________________ 72
3.2 Numerical accuracy for high contrast simulations___________________________________ 74
3.3 Robustness of the linear and non-linear system solvers_______________________________ 74
4. NUMERICAL RESULTS AND ANALYSES _________________________________ 75
5. CONCLUSIONS AND OUTLOOK _________________________________________ 76
6. Acknowledgements_______________________________________________________ 77
7. REFERENCES __________________________________________________________ 77
Modélisation stochastique de l'intrusion saline en 2D plan 57
Résumé
Ce chapitre est constitué de deux publications en anglais : un chapitre de livre (Albitar et Ababou,
2005) et un article de conférence (Albitar et Ababou, 2006). Dans ce chapitre on présente des
résultats numériques et analytiques de modélisation stochastique de l’intrusion saline dans des
domaines 2D plan hétérogènes jusqu'à fortement hétérogènes. L’objectif principale est de déduire
une relation entre l’incertitude de la position de l’interface eau douce/eau salé et l’hétérogénéité du
milieu poreux représenté à travers de la transmissivité hydraulique. Dans la première partie du
chapitre l’accent est mis sur les équations mathématiques, la solution numérique et la solution
analytique stochastique du problème. Tandis que dans la seconde partie, les méthodes numériques
utilisées pour la modélisation de l’écoulement en milieu fortement hétérogène sont présentées.
L’annexe D complète ce chapitre on pressentent la méthode spectrale, et les résultats numériques
complets.
Modélisation stochastique de l'intrusion saline en 2D plan 58
Seawater intrusion in coastal aquifers is a growing concern in Mediterranean regions, due to over-
population and over-exploitation of coastal groundwater resources. Under these circumstances, it is
essential to model the extent of seawater intrusion and to locate the saltwater-freshwater interface
taking into account heterogeneity and parameter uncertainty. There are different ways to couple salt
transport and freshwater flow in groundwater models. We choose here the vertically integrated sharp
interface approach, with two immiscible fluid regions (freshwater and seawater). We use this model
to analyze the effects of aquifer variability on the saltwater wedge in plane view, based on large
numerical simulations of 2D seawater intrusion in randomly heterogeneous unconfined aquifers.
We consider an unconfined coastal aquifer with an impervious bedrock at z = ZINF(x,y) and a fresh
water table of elevation z = ZS(x,y). In addition, because we use a plane flow model, the vertically
averaged freshwater hydraulic head H(x,y) coincides with the free surface elevation, i.e. :
H(x,y)≈ZS(x,y). In this 2D framework, all variables and parameters are spatially distributed in (x,y).
A schematic representation of the coastal aquifer and its salt wedge is shown in Fig. 1.
We assume that seawater and freshwater are separated by a sharp interface. More precisely, we
rely on the Ghyben-Herzberg approximation(s), that is:
• the seawater and freshwater fluids are assumed immiscible (sharp interface);
• the subsurface seawater wedge is assumed quasi-hydrostatic;
• the freshwater is assumed vertically hydrostatic (negligible vertical velocities).
Modélisation stochastique de l'intrusion saline en 2D plan 59
Fig. 1 Schematic view of seawater intrusion (sea level ZSEA shown at right) into a free surface aquifer (shown at left),
with saltwater interface ZSALT(x,y) and substratum ZINF(x,y)
Now, let ZSALT(x,y) be the elevation of the salt/fresh water interface. Applying the hydrostatic
assumptions and the pressure continuity condition at the interface, and modifying the Badon-
Ghyben-Herzberg configuration to account for a finite outflow face of height ΔZ located undersea,
we obtain:
ρF g (H − ZSALT ) = ρS g (ZSEA − ZSALT − ΔZ ) + ρF g ΔZ (26)
This gives finally the desired closure relation:
Z SALT = Z SEA −
(H − Z SEA ) − δZ
ε (27)
In these equations, ρF is freshwater density, ρS is saltwater density, and ε is the saltwater-to-
freshwater density contrast:
ρS − ρF 1
ε= ≈ (28)
ρF 40
Parameter ΔZ is the vertical depth of the freshwater outflow face at the shoreline, assumed much
smaller than aquifer thickness. It can be obtained from exact solution of seawater intrusion in a
vertical slice (x,z) of a homogeneous confined aquifer, without depth-averaging. Here, ΔZ is about
0.77 m, compared to 30 m aquifer thickness.
For freshwater flow, we use the Dupuit-Boussinesq plane flow approximation. The freshwater
thickness is defined as:
η ( x, y) = H ( x, y) − Zinf ( x, y) or η(x, y) = H(x, y) − Zsalt(x, y) (29)
depending on spatial location (x,y), within the salt wedge or not. The freshwater transmissivity
T(x,y) is then inferred from freshwater thickness η(x,y) as follows :
T ( x, y ) = K ( x , y ) × η ( x, y ) (30)
Note that T is spatially variable via K and ZINF , and also, nonlinear via the unknown variables H
and ZSALT on which it depends. Finally, we obtain the following system of vertically integrated flow
equations (steady state case):
1. Steady-state mass conservation (freshwater):
∂Θ
= −div (Q) where Θ is the water content (31)
∂t
2. Darcy’s law (vertically integrated):
Q = -T (H,Z SALT ,Z INF ,x,y ) grad(H ) (32)
Modélisation stochastique de l'intrusion saline en 2D plan 60
3. Freshwater transmissivity:
K(x, y) ×(H − ZINF(x, y)) if ZSALT < ZINF
T= (33)
K(x, y) ×(H − ZSALT(x, y)) if ZSALT ≥ ZINF
In this paper, we choose to study uncertainty without regard for specific data. That is, we choose to
simulate seawater intrusion in large unconditional single-replicates of the heterogeneous aquifer. We
use the XIMUL code to generate isotropic log-normal random fields K(x,y) on 1 million node grids
(1000º×º1000).
The XIMUL code deals more generally with Bayesian estimation and conditional simulation of
1,2,3-D random functions of space or time [Ababou et al 1994]. The unconditional generator uses
the Fourier Turning Band method based on a representation theorem of Matheron (1973): see
[Tompson et al. 1989] and references therein.
Numerical simulations of seawater intrusion are carried out using the BigFlow code BF 2000
[Ababou and Trégarot 2002]. It solves a generalized model equation for flow in heterogeneous,
anisotropic, partially saturated media. It can efficiently follow multiple interacting free surfaces in
3D, and it can represent “open” or “macroporous” media [Trégarot 2000]. A vertically integrated 2D
flow module is also available, including Boussinesq-Dupuit aquifer flows, free surface hydraulics
based on kinematic-diffusive wave, Darcy-Forchheimer flow in rough fractures [Spiller 2004], and
the seawater intrusion module SWIM2D used here.
The BF 2000 code is based on implicit 3D finite volume formulation of flux divergence equations
in conservative form (mixed form). It solves fully coupled transient and steady flow problems, using
a single infinite time step for steady state. It uses Preconditioned Conjugate Gradients for matrix
solution, and modified Picard iterations for nonlinear solution. The matrix-vector data structure is
very sparse. For more details on the numerics, see [Ababou et al 1992; Ababou and Bagtzoglou
1993; Ababou 1996].
Table 1 Summary of statistical parameters for two sets of simulations (small and large)
Parameters Set 1 Set 2
ni (number of nodes) 300 × 300 1000 × 1000
Fig. 2 and Fig. 3 display perspective views of simulated seawater intrusion for a highly variable
permeability (σlnK = ln10 ≈ 2.30). Two surfaces are displayed in each figure : ZSALT(x,y), the
salt/fresh interface level (mapped with color-coded logK values), and ZS(x,y), the freshwater
piezometric surface (or hydraulic head), also mapped with the same color-coded or grey-scale logK
values.
On Fig. 2, one can clearly observe the sharp local gradients of the saltwater interface occuring in
low permeability zones, which act as barriers to seawater (it should be kept in mind, however, that
K(x,y) is the depth-averaged permeability).
Fig. 4 depicts the effect of heterogeneity level on the mean penetration of the salt wedge, for a
300 × 300 grid. The mean ZSALT(x) profile is plotted versus distance from sea (x), after averaging
ZSALT(x,y) along the shorewise direction (y). The three profiles correspond to: σlnK = 0
(homogeneous), σlnK = 1 (moderate heterogeneity) and σlnK = ln10 (high heterogeneity). As the level
of variability σlnK increases, the mean elevation ZSALT(x) increases and the mean salt wedge
penetrates farther inland. The extra penetration of the mean wedge due to heterogeneity is about
200 m, for the most heterogeneous case.
A similar result (not shown here) was obtained for the larger 1000 × 1000 grid with heterogeneity
levels σlnK = 0, 1.0, √2, 1.6, and ln10. It confirms the monotonic increase of the mean penetration
length of the salt wedge as σlnK increases, compared to a homogeneous aquifer with geometric mean
permeability.
Modélisation stochastique de l'intrusion saline en 2D plan 62
As a first step towards uncertainty analysis (next section), let us develop further the statistical
analysis of the simulated salt wedge, based on single replicate unconditional simulations obtained on
the largest grid (1000 × 1000 cells).
The salt wedge is characterized by the shape of the saltwater interface elevation ZSALT(x,y) and
its horizontal extension inland. We consider ZSALT(x,y) as a random field and we analyze it
statistically. We focus in particular on the first and second order moments of ZSALT(x,y), including
its mean and its standard deviation. This analysis is applied to the 1000x1000 grid, with large
variability (σlnK = ln10).
Given the symmetries of the problem and the statistical stationarity of K(x,y), we expect the
surface ZSALT(x,y) to be stationary (statistically homogenous) along the y direction parallel to the
seashore. However, it will not be stationary along the x direction parallel to flow (transverse to
seashore).
Indeed, Fig. 5 shows 100 transects ZSALT(x,yn) along with the average profile, all plotted as
functions of the x-coordinate (perpendicular to sea shore). The profiles ZSALT(x) are clearly non-
stationary.
It is clear from both Fig. 4 (mean ZSALT) and Fig. 5 (random ZSALT) that the interface elevation
ZSALT(x,y) follows a nonlinear trend along x (for fixed y) and cannot be a stationary random function
of x. This observation has two consequences:
1. Given a single replicate of the coastal aquifer in (x,y), we can only sample in the shorewise
direction (y) to produce a statistical description of the salt wedge.
2. For theoretical purposes, we may seek a convenient transformation ZSALT → Φ to obtain an
approximate stationary field Φ from the non-stationary field ZSALT.
Following this idea, consider first the analytical solution of the homogenous problem (σ = 0) using
the Badon-Ghyben-Herzberg assumption, modified to include a submarine outflow face of height
ΔZ at the seashore :
(ε + 1) ⎛ ε Δ Z ⎞
2
⎛ ε ΔZ ⎞ x H 12 − Z SEA
2
Z SALT ( x ) = Z SEA − ⎜ ⎟ − ΔZ − +⎜ ⎟
⎝ ε +1 ⎠ ε LX ε (ε + 1) ⎝ ε +1 ⎠ (34)
Equation (34) holds for 0 ≤ x ≤ LSALT, where x=LSALT is the intersection of the saltwater interface
with the substratum. Thus, if the bedrock is at z = 0, the value of LSALT is defined by ZSALT (x) = 0.
Other variables in equation (34) are defined below:
• LX is the domain size in the x-direction, between the two fixed head boundaries H=H0 (sea at left)
and H=H1 (freshwater at right);
• ZSEA = H0 is the elevation and depth of the sea level above the substratum, at the sea boundary
x=0;
• H1 is the depth of the freshwater level above the substratum at the inland boundary x=LX ;
• ZSALT is the elevation of saltwater/freshwater interface above the substratum;
• LSALT is the x-wise penetration length of the salt wedge inland, on the bedrock;
• ΔZ is the vertical length of the submarine freshwater outflow face into the sea.
Modélisation stochastique de l'intrusion saline en 2D plan 63
The term (ε ΔZ)/(ε+1) can sometimes be neglected in equation (9); we have here:
(ε ΔZ)/(ε+1) = 0.018 m ,
which is indeed small compared to ZSEA = 30 m and to ΔH = H1 -ZSEA = 1 m.
Also, to simplify the above expression Zsalt(x), we define a new parameter h0 :
H12 − ZSEA 2 (ε + 1) (35)
h0 =
ε ⋅ (ε + 1)
This parameter, h0 , is a length scale on the same order as the thickness of the freshwater lens
imposed at the inland boundary (upstream). The solution of the homogenous problem σ = 0 can now
be expressed as:
h0 x
Z SEA − δZ − Z SALT ( x) ≈ ⋅ (36)
ε LX
This simple analytical expression (36) shows that there exists, for the homogeneous case, a
quadratic transform which makes the saltwater profile exactly linear in x. For the heterogeneous
case, with random field permeability, this suggests applying the same quadratic transform to the
nonlinear random function ZSALT(x,y). The transformed field is a new random “potential” field
φSALT(x,y) with:
Φ SALT = ( Z SEA − δZ − Z SALT ) 2 (37)
We may expect that the random ΦSALT(x,y) has a roughly linear trend. Furthermore, it is possible to
derive analytically the mean and variance of ZSALT from the moments of the random field
ΦSALT(x,y). Let us first normalize ZSALT and ΦSALT by ZSEA as follows: Z=(ZSALT-δZ)/ZSEA;
φ = ΦSALT/ZSEA. The Φ-transform is now:
φ (x, y ) = (1 − Z (x, y ))2 (38)
with φ = 0 (exactly) on the sea boundary x = 0, and φ = 1 at some fixed distance L1, the characteristic
length of penetration of the salt wedge. The latter is given, to order O(σ), by the analytical solution
for a homogeneous aquifer:
L1 = LSALT (σ ) ≈ LSALT (0 ) × (1 + O(σ ))
(39)
Thus, we may write the (approximate) boundary condition of the random case as:
x = 0 : φ = 0; x = L1 : φ ≈ 1 + O(σ ) (40)
The main idea, here, is that we prefer to solve for the Φ-field because it is more easily amenable
to statistical analysis than the Z-field (more on this below). With this goal in mind, let us define the
random fluctuations of φ and Z:
ϕ ( x, y ) = φ ( x, y ) − φ ( x) and z ( x, y ) = Z ( x, y ) − Z ( x) (41)
where the mean potential is given by: φ ( x) = 〈(1 − Z )²〉 .
The brackets <•> represent either the shorewise spatial average (spatial mean of a single replicate
along direction “y”), or the mathematical expectation E(•) over an ensemble of replicates : the two
are equivalent if ergodicity is assumed.
Now, substituting the random fluctuations in equation (38) and taking averages, we obtain:
( ) 2
φ = 1 − Z + z ² and σ Z 2 = φ − 1 − Z ( ) 2
(42)
where the mean <Z> remains to be determined. On the other hand, from equation (38) :
Z = 1 − φ 1 / 2 (for 0<x<L1 and 1>Z>0). (43)
Z = 1 − φ 1/ 2 (1 + κ )1 / 2 with ϕ <1
κ = (44)
φ
Modélisation stochastique de l'intrusion saline en 2D plan 64
Fig. 2 Perspective view of ZSALT (x,y), H(x,y), and log K(x,y) for a gauss-
shaped isotropic covariance with σ = ln10 and L/λ = 30. Simulation grid:
300×300
Fig. 3 Perspective view of ZSALT (x,y), H(x,y), and log K(x,y) for a gauss-shaped
isotropic covariance with σ = ln10 and L/λ = 100. Simulation grid: 1000×1000
Modélisation stochastique de l'intrusion saline en 2D plan 65
30
homogenous field
Zsalt(x) Mean for N(1,1)
Zsalt(x) Mean for N(1,ln(10))
25
20
Zsalt (m)
15
10
0
0 100 200 300 400 500 600 700 800 900 1000
x (m)
We finally substitute Eq. (45) into Eq. (42) to calculate the standard deviation of Z. Neglecting
τ /64 and other “higher order terms” (“h.o.t.”), we obtain:
4
1 σφ
σZ ≈ + h.o.t. (46)
2 φ
This analytical expression can be used to predict σZ using either numerical estimates or theoretical
spectral estimates of φ-statistics : the two procedures yield similar results (see comments about
Fig. 6 further below).
We know need to determine the statistical moments of Φ, e.g. mean and variance. Two approaches
are possible concerning the transformed potential Φ:
a) Empirical evaluation of Φ-moments (sampling numerical simulation);
b) Theoretical evaluation of Φ-moments (analytical spectral perturbation).
Empirically, the first two lines in Table 2 show some of the numerically computed moments of
ΦSALT , assuming a linear trend <Φ>, and stationary fluctuations ϕ(x,y) around the linear trend:
Φ (x, y ) = Φ (x ) ≈ ax ϕ (x, y ) ≈ Φ (x, y ) − ax (47)
ϕ (x )2
1/ 2
σ Φ = ≈ constant (48)
Note: These relations hold only in a subdomain comprised between the sea boundary x = 0 (where
φ = 0) and the tip of the salt wedge x ≈ L1+O(σ) (where φ ≈ 1+O(σ)).
On the other hand, we demonstrate that the Φ-equation in the salt wedge zone is a stochastic PDE,
analogous to the Boussinesq equation for vertically averaged groundwater flow with random K(x,y).
Indeed, from eqs.(30), (31), (32), we have:
∂ ⎛ ∂H ⎞
⎜ − K ( x1 , x2 ) ( H − Z SALT ) ⎟ (i=1,2) (49)
∂xi ⎜ ∂xi ⎟
⎝ ⎠
The freshwater head H is given by the Ghyben-Herzberg relation Eq.(2.68):
H = (1 + ε ) Z SEA − ε ΔZ − ε Z SALT (50)
Substituting H in Eq.(49), and using the Φ-transform, we obtain:
∂ ⎛ ∂φ ⎞
⎜ − K ( x1, x2 ) ⎟ = 0 (i=1,2) (51)
⎜
∂xi ⎝ ⎟ ∂xi ⎠
We observe that this φ-equation is equivalent to a stochastic groundwater flow equation with 2D
random field transmissivity in a confined aquifer (cf. “infinite domain” spectral perturbation
solutions by [Mizell et al. 1982]). Thus, σΦ can be evaluated from the spectral solution of Eq.(51), at
least far enough from the sea and the saltwedge tip. The “theoretical” standard deviation of φ is
deduced from the Mizell et al (1982) solution, for a “modified Wittle” correlation structure:
(σ Φ )THEORY ≈ c σ ln K λln K J x ≈ c σ ln K λln K a (52)
where Jx is the mean Φ-gradient denoted “a” in this paper. The coefficient “c” is a dimensionless
constant of order 0(1) [Mizell et al 1982]. For the problem at hand, the value of “c” can be obtained
by matching numerical and theoretical “σΦ” at low levels of heterogeneity (σlnK ≤ 1). This procedure
gives:
Modélisation stochastique de l'intrusion saline en 2D plan 67
c ≈ 1.10 (53)
Similarly, the relevant value of the mean Φ-gradient, a = <dΦ/dx>, can be obtained from the exact
analytical solution Φ(x) in a homogeneous aquifer, which corresponds to the asymptotic case
σlnK → 0. Thus, asymptotically:
h 02
σ → 0: (a )THEORY = a 0 + O (σ ) = + O (σ ) (54)
ε Lx
σ ln K 0 1 1.60 2.30
σ̂ Φ NUM σΦ = 0 σˆ Φ ≈ 17 σˆ Φ ≈ 27 σˆ Φ ≈ 42
To check whether “a” is nearly constant and close to its predicted value “a0”, consider the results
summarized in Table 2. We conclude that the theoretical prediction of σΦ given by equation (52)
with a ≈ a0 is robust.
Finally - after some manipulations involving statistics from the Z-Φ transform (Eqs.(45),(46)) and
the spectral solution for σΦ - one obtains, to first order:
(a): σZSALT (x) ≈ c σlnK λlnK a or (b): σ ZSALT (x) ≈ c σ ln K λln K ∂Z (55)
2 x 2 ∂x
Both versions of this equation require mean gradient information: the first equation (a) requires
knowledge of the (stationary) mean φ-gradient “a”, while the second version (b) requires knowledge
of the (non-stationary) mean interface elevation gradient.
Fig. 6 shows 100 superimposed transects of the “potential” ΦSALT(x,yn), sampled at equally spaced
shorewise locations “yn”, and plotted versus (x), for σ = ln10. The figure also shows the analytical
profile ΦSALT(x) for a homogeneous aquifer (σ = 0), as well as the numerical average of ΦSALT(x,y).
The fluctuations of ΦSALT(x,y) around its mean trend were also plotted as transects (not shown here).
These numerical plots indicate the level of fluctuation of the salt interface in terms of the
transformed field ΦSALT. They also confirm the quasi-linear trend of ΦSALT.
We computed the fluctuations of ZSALT around its nonlinear mean trend, and we estimated σZ by
sampling ZSALT parallel to the seashore and plotting the resulting moment σZ as a function of
distance (x) from the sea. One result is shown in Fig. 7 for large heterogeneity (σ=ln10). The
standard deviation of ZSALT seems approximately stationary far enough from the seashore (x = 0) and
far enough from the salt wedge tip (x ≈ 700 m). In the stationary region of Fig. 7, we find
σZSALT ≈ 1.3 m. The 95% confidence band of the salt interface is several meters, which represents a
rather significant fraction of the total aquifer thickness of 30 m.
Modélisation stochastique de l'intrusion saline en 2D plan 68
Fig. 6 One hundred transects of ΦSALT (transformed from ZSALT); “analytical mean” curve ΦSALT (homogeneous
aquifer); and "numerical mean" curve ΦSALT (mean of ΦSALT sampled shorewise along “y”). The sea shore is at left.
Grid: 1000x1000. Heterogeneity: σ = ln(10)
The results appear different for lesser heterogeneity: see Fig. 8 for σ = 1.60, and note that similar
results were obtained for σ in the range 0 ≤ σ ≤ 2.0. In all these cases, σZSALT(x) is non-stationary
with respect to (x) and decreases with (x), as predicted by the theoretical Φ-transform analysis. This
can be seen by comparing the “numerical” and “analytical” (Eq.(55).a) σZSALT(x) curves in Fig. 8.
Fig. 7 Standard deviation of ZSALT vs. distance (x) from seashore (sea located at left), obtained by sampling ZSALT
fluctuations in the shorewise direction (y). The global value of σZ appears to be about σZ ≈ 1.3 m. Grid: 1000x1000 cells.
Heterogeneity: σlnK = ln10 = 2.30
Modélisation stochastique de l'intrusion saline en 2D plan 69
Fig. 8 Numerical and theoretical σZSALT vs distance from sea (x) for σlnK = 1.60
8 Acknowledgment
This study is part of the european project SWIMED on coastal aquifer management, funded by the
European Commission (Sustainable Water Management In MEDiterranean coastal aquifers):
[Link]
9 References
Ababou R (1996) Random Porous Media Flow on Large 3D Grids: Numerics, Performance, and
Application to Homogenization, Chap.1, pp.1-25. In: IMA Vol 79 Mathematics and its
Applications: Environmental Studies (Math. Comput. Statist. Anal.). Wheeler MF (ed.),
Springer, NY, 410 pp.
Ababou R, Bagtzoglou AC (1993) BIGFLOW: a Numerical Code for Simulating Flow in Variably
Saturated, Heterogeneous Geologic MediaTheory and User’s Manual Ver.1.1. NUREG/CR-
6028, US NRC Report, Washington DC.
Modélisation stochastique de l'intrusion saline en 2D plan 70
Ababou R, Bagtzoglou AC, Wood EF (1994) On the Condition Number of Covariance Matrices
Arising in Kriging, Estimation & Simulation of Random Fields. [Link].26(1), pp. 99-133,
1994.
Ababou R, Trégarot G (2002) Coupled Modeling of Partially Saturated Flows : Macro-Porous
Media, Interfaces, and Variability. Proc. CMWR 02, Comput Meth Water Resour, 23-28 June
2002, Delft, The Netherlands, Elsevier, 8 pp.
Ababou R, Sagar B, Wittmeyer G (1992) Testing Procedures for Spatially Distributed Flow Models.
Advances in Water Resources, Vol.15, pp. 181-198, 1992.
Matheron G. (1973) The Intrinsinc Random Functions & Applications. [Link]., 5, 439-468
Mizell SA, Gutjahr AL, Gelhar LW (1982) Stochastic Analysis of Spatial Variability in Two-
Dimensional Steady Groundwater Flow Assuming Stationary and Nonstationary Heads. Water
Resour Res 18(4) 1053-1067.
Spiller M (2004) Physical and Numerical Experiments of Flow and Transport in Heterogeneous
Fractured Media : Single Fracture Flow at High Reynolds and Reactive Particle Transport. PhD
thesis, Aachen Univ. (Germany) & Institut Nat. Polytech. Toulouse (France), October 2004.
Tompson AFB, Ababou R, Gelhar LW (1989) Implementation of the Three-Dimensional Turning
Bands Random Field Generator. Water Resour. Res., 25(10), 2227-2243.
Trégarot G (2000) Modélisation Couplée des Ecoulements à Saturation Variable avec
Hétérogénéités, Forçages, et Interfaces Hydrologiques. PhD thesis, Institut Nat. Polytech.
Toulouse, May 2000.
Modélisation stochastique de l'intrusion saline en 2D plan 71
1 INTRODUCTION
Seawater intrusion is a well known problem occurring in coastal aquifers. This phenomenon
can take the form of a seawater wedge extending inland below freshwater. Here we focus on
the effects of permeability heterogeneity on the extent of the seawater wedge in a phreatic
coastal aquifer. Seawater intrusion is analysed using 2D plane flow Dupuit-Boussinesq
approximation, combined to a sharp interface approach (Ghyben-Herzberg). Aquifer
heterogeneity is represented via a 2D random field permeability (vertically integrated). The
simulations are conducted for mildly to highly heterogeneous random fields, with natural
log-permeability standard deviation (σlnK) ranging from 1.0 to 4.0. This paper focuses on the
need for numerical accuracy in the case of highly heterogeneous aquifers. An optimized
continuation method for modelling flow in highly heterogeneous domains is introduced and
implemented in a Graphical User Interface (GUI) environment. While the domain geometry
is simple, the main computational challenge is to solve the seawater intrusion problem on
large numerical grids, with both nonlinear and highly variable coefficients. The numerical
results are then analysed in terms of spatial statistics, such as the standard deviation of the
seawater interface elevation versus distance from the seashore (σZSALT(x)).
Seawater intrusion is modeled in the framework of the sharp interface approach, while still
retaining the essential feature of aquifer variability in the 2D plane. Our approach assumes
that seawater and freshwater are immiscible fluids, and relies on pressure equilibrium
relations (Ghyben-Herzberg) to close the resulting system of equations. In this 2D
Modélisation stochastique de l'intrusion saline en 2D plan 72
where H is the total freshwater hydraulic head (m), ZSALT is the salt/freshwater interface
level (m), g is the gravitational acceleration (m s-2), ρF is the freshwater density (kg m-3), ρS
is the saltwater density (kg m-3), and ΔZ is the vertical depth of the freshwater outflow face
at the shoreline (m). The Ghyben-Herzberg closure relation (above) is inserted in the
vertically-averaged Boussinesq equations for freshwater, leading to a nonlinear system of
equations, namely : steady-state mass conservation for freshwater (2); vertically integrated
Darcy’s law (3); & freshwater transmissivity (4) :
∂Θ
= −div(Q) (2)
∂t
where K(x,y) is hydraulic conductivity (m s-1), ZINF is the impervious substratum level (m),
ZSALT is the salt/Freshwater interface level (m), Q is the Darcy specific discharge rate (m² s-
1
), Θ is the water content (m3 m-3), and T the transmissivity (m2 s-1). Here we focus on
steady state problems. For more details, see Ababou & Al-Bitar (2005).
Numerical simulations of seawater intrusion are carried out using the BigFlow code
BF 2000 (Ababou & Trégarot, 2002). The model is based on a single generalized flux-divergence
equation (conservative, mixed form) for either 3D or 2D plane flow. Space-time discretization is
based on implicit finite volumes, leading to sparse nonlinear systems. These are solved using two
interspersed loops: the outer loop implements fixed point (modified Picard) iterations; the inner loop
solves the linearized system by preconditioned Conjugate Gradients (DSCG). The 2D plane flow
option was specialized for sharp interface seawater intrusion problems (BF-SWIM2D).
Some special considerations were given to the non-linear aspects of the SWIM problem.
Indeed, Seawater Intrusion Modeling is “doubly non-linear” because of the non-linearity of
the Dupuit Boussinesq equation (head-dependent transmissivity) and because of the
nonlinearity due to the interception of the salt interface by the aquifer’s substratum. This
non-linearity is even stronger when the domain is highly heterogeneous (K(x,y). For these
cases, we used a special iterative “re-start method”, also known as “continuation” or
“homotopy” method. The continuation parameter is chosen to be σlnK, the standard
deviation of lnK.
Modélisation stochastique de l'intrusion saline en 2D plan 73
σi ln(i + 1)
= (5)
σ N ln( N + 1)
where σi is the standard deviation of the ith lnK field, N is the number of increments, σN is
the standard deviation of the objective (final) standard deviation. Other incremental
functions have been tested; the logarithmic increment proved to be the most efficient.
Fig. 1(a) shows the values of σlnK for 10 continuation steps (final value σlnK = 4.0). The
results are analysed in Fig. 1(b), in terms of error norm vs iteration number. The root-mean-
square norm of δH is used, where δH is the head variation over two successive iterations.
This δH is to be compared to the total head variation of 1 meter along the mean flow path.
The 2D random seawater intrusion problem was solved on a 1000x1000 node grid. The total
CPU time for the 10 continuation steps was on the order of one day on a PC Intel Pentium-4
processor, in double precision.
Inside each continuation step, there are two interlooped solvers: (1) the inner loop is the
iterative DSCG solver (stopping criterion δH < ε =10-9 (m), and ITER ≤ 1000); (2) the outer
loop is the fixed point (Picard) iteration solver (stopping criterion δH < ε = 10-6 and
ITER ≤ 20). In Fig. 1(b), the average number of iterations within each continuation step is
18596. However the actual number of iterations per step increases with successive steps.
The problem is more difficult as σlnK increases.
The continuation algorithm, as well as file management and graphic tasks, have been
implemented as an option in the GUI application (BF-Py) under development for the
BigFlow code. The BF-Py application has been implemented under Python® 3 and using the
wxPython® library.
Fig. 1 a Evolution of σlnK in the homotopy method with N = 10 and σN = 4.0 (left)
Fig. 1 b Error norm for ten continuation steps (tens of thousands of iterations).
One issue in simulating highly heterogeneous flows is to ensure numerical accuracy. For
this purpose, double precision computations need to be implemented, rather than single
precision. To inspect the solver’s robustness under double precision computations, a simple
configuration with high contrast field properties was used, as shown in Fig. 2(a). The
computational domain consists of a square domain with constant permeabilities in which a
bar shape with a different permeability is inserted perpendicularly to the flow direction and
over the entire width. Fixed heads are applied in the direction of flow, and no flux Neumann
conditions perpendicular to the flow.
This heterogeneous “bar shape” configuration has a lnK contrast that can be expressed
in terms of a geometric mean and a standard deviation (Ababou, 1988):
⎡⎛ ⎧ ( K1 )α1 ⎫ ⎞ 2 ⎛ ⎧ ( K 2 )α 2 ⎫ ⎞ 2 ⎤
K G = ( K1 ) ( K 2 ) σ ln K = ⎢⎜⎜ ln ⎨
α1 α2
⎬ ⎟⎟ + ⎜⎜ ln ⎨ ⎬ ⎟⎟ ⎥ (6)
⎣⎢⎝ ⎩ G ⎭ ⎠ ⎝ ⎩
K K G ⎭ ⎠ ⎦⎥
where αi is the area fraction of the domain associated with permeability Ki (α1+α2=1).
We conducted seawater intrusion simulations for several values of the permeability
contrast (K1/K2 or σlnK) and we analysed mass balance (not shown here for lack of space).
Mass balance remains quite good even for extremely high contrasts: the relative error on the
net computed flux is 0.5 % for K1/K2 = 108 or σlnK = 5.58.
Fig. 2(b) shows the analytical solution and the numerical hydraulic head transect H(x) in the
middle of the domain, for K1/K2 = 10+8; the local relative error is less than 10-7.
To complete the non-linear study, we examine in more detail the robustness of the DSCG
solver interlooped with the nonlinear Picard solver. First, note that a steady simulation with
BigFlow consists of one infinite time step. This is obtained by setting the mass storage term
in the linearized equation to zero. The Picard iterations are used to linearize the resulting
system.
Modélisation stochastique de l'intrusion saline en 2D plan 75
Let us consider the 1D homogeneous case. In this case the linearized system is tri-diagonal.
The CG-based solver should theoretically converge with a number of iterations smaller than or equal
to the matrix dimension (Golub & Van Loan, 1989).. This theoretical result was observed
experimentally with our double precision simulations using the DSCG solver.
Fig. 3(a) illustrates the results of the CG solver for one outer loop of the Picard iteration.
The grid has 1000×1×1 internal nodes. The results show that convergence is slow until the 999th
iteration is reached; then convergence occurs abruptly at the 999th iteration, within machine
precision. This behaviour was also observed by Ababou (1996) for linear saturated flow problems.
Fig. 3(b) shows the numerical results for all interlooped iterations (Picard and DSCG). There
is an increase in error at the beginning of each new Picard iteration, due to updates of the matrix
coefficients. However, this error decreases with iterations, i.e., there is a global convergence of the
nonlinear Picard iterations.
Fig. 3(a) Error norm of the DSCG matrix solver within a single Picard step; notice the abrupt
convergence occuring at the 999th iteration (left).
Fig. 3(b) Error norm of the DSCG matrix solver for all Picard iterations (right).
We now analyse the spatial statistics of the seawater wedge (interface elevation ZSALT(x,y))
for the 10 different levels of aquifer variability, as specified in the numerical continuation
method. This is shown in Fig. 4(a) and Fig. 4(b).
Modélisation stochastique de l'intrusion saline en 2D plan 76
Due to the randomness of the vertically integrated permeability K(x,y), the resulting
seawater wedge is very heterogeneous. This can be seen from Fig. 4(a), where ZSALT(x,y) is
shown for the case of largest variability (σlnK = 4.0). In addition, Fig. 4(b) shows the spatial
moment σZSALT(x) plotted as a function of distance from the seashore, for σlnK = 1.0 up to
σlnK = 4.0. This extends results by Ababou & Al-Bitar (2005); they analysed σZSALT(x) using
a perturbation approximation, valid up to σlnK = 1.6.
6 Acknowledgements
This study is part of the European project SWIMED on coastal aquifer management, funded
by the European Commission (Sustainable Water Management In MEDiterranean coastal
aquifers): [Link] .
7 REFERENCES
Ababou R (1988) Three-dimensional flow in random porous media. Ph.D. thesis. M.I.T., Cambridge MA, USA.
Ababou R. (1996) Random Porous Media Flow on Large 3-D Grids: Numerics, Performance, and Application to Homogenization.
Chap.1, pp.1-25, in IMA Vol. N .79 Mathematics and its Applications: Environmental studies - mathematical, computational
and statistical analysis. [Link] ed., Springer, New-York, 410 pp.
Ababou R. & Al-Bitar A. (2005) Random field approach to seawater intrusion in heterogeneous coastal aquifers: unconditional
simulations and statistical analysis. In: Geostatistics for Environmental Applications, Renard P., Demougeot-Renard H.,
Froidevaux R. (eds.), Springer Verlag.
Ababou R & Trégarot G (2002) Coupled Modeling of Partially Saturated Flows : Macro-Porous Media, Interfaces, and
Variability. In: Proc. CMWR 02, Comput. Meth. Water Resour., June 2002, Delft, The Netherlands, Elsevier, 8 pp.
Ababou R & Bagtzoglou A.C. (1993) BIGFLOW: a Numerical Code for Simulating Flow in Variably Saturated, Heterogeneous
Geologic Media - Theory and User’s Manual Ver.1.1. NUREG/CR-6028, Washington DC., USA.
Tompson AFB, Ababou R, Gelhar LW (1989) Implementation of the Three-Dimensional Turning Bands Random Field Generator,
2227-2243. Water Resour. Res. 25(10).
Golub G.H. & Van Loan C.F. (1989) Matrix computations. 2nd Ed., J. Hopkins Univ. Press,
Baltimore, Maryland, 642 pp.
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 78
Chapitre V MODELISATION 2D
COUPLEE SURFACE / SOUTERRAIN
AVEC OU SANS INTRUSION SALINE
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 79
V - 1 INTRODUCTION
V - 1.1 Introduction
Dans plusieurs configurations, l’écoulement dans un système hydrologique est régi par les
interactions entre l’écoulement en surface et l’écoulement souterrain. Dans ce type d’écoulement il
est important de prendre en compte le couplage entre le milieu souterrain et les volumes d’eau en
surface (rivière, lacs, ruissellement…).
Le couplage peut se faire implicitement ou itérativement.
Le couplage itératif consiste à faire des échanges entre les processus physiques à des instants
déterminés. Le temps d’échange dépend du temps caractéristique de chaque processus et de la
difficulté à résoudre le problème pour les codes à pas de temps adaptatif. Cette approche permet de
faire du calcul séquentiel ou parallèle indépendant durant un ou plusieurs pas de temps. Une telle
approche permet de réduire les demandes en mémoire de stockage, mais elle présente un algorithme
moins stable car la solution de chacun des systèmes non-linéaires est obtenue en considérant une
solution constante de l’autre processus durant un pas de temps donné. Ainsi les pas de temps doivent
être petits, et la demande en temps calcul peut s’avérer coûteuse. A cet inconvénient vient s’ajouter
le temps d’échange des données.
La méthode de couplage implicite consiste à écrire les équations de chacun des processus en
considérant des termes d’échange (sources/puits) entre les différents processus et à les discrétiser.
Ensuite, le système d’équations doit être assemblé dans un système algébrique global. L’inversion
de la matrice globale de ce système, correspondant aux deux processus simultanément, permet
d’obtenir une solution plus stable que dans le couplage itératif. Par contre, cette méthode est
coûteuse en mémoire et elle implique un temps calcul commun aux deux processus. Des méthodes
de décomposition de domaine et de calcul parallèle permettent de pallier à ces problèmes.
Une autre façon de faire du couplage implicite est de considérer une équation unique dont les
coefficients sont modifiés dans plusieurs parties du domaine de calcul, associés à des processus
physiques différents. Cette approche est inspirée principalement de l’analogie entre les équations
décrivant l’écoulement en milieu souterrain et en surface. Elle permet d’assurer une certaine
continuité des variables telles que la pression aux interfaces des différentes zones d’écoulement.
C’est cette approche qu’on utilise dans ce chapitre.
Le problème du couplage est encore plus complexe lorsque le système rivière/souterrain est sujet
à l’intrusion saline, comme c’est le cas dans les aquifères côtiers. La méthode la plus utilisée pour
résoudre le problème de l’intrusion saline est l’approche de zone de diffusion, où l’on couple les
équations d’écoulement à densité variable et l’équation de transport de sels (solutés). Cette dernière
approche est passée en revue dans le chapitre 2. Or, en présence d’un écoulement de surface à
cinétique rapide, ces modèles couplés écoulement-transport convergent plus difficilement.
L’approche interface abrupte à deux fluides (eau salée / eau douce) utilisée ici est plus adéquate.
Cette approche conduit au couplage de deux équations d’écoulement, mais si l’on considère de plus
que la zone salée est immobile, le système se réduit à une seule équation avec des coefficients non-
linéaires et variables suivant la configuration du système.
∂( ηU y ) ∂( ηU x U y ) ∂ ( ηU y )
2
⎡ ∂Z ⎤
+ + = − gη ⎢ s + Sfy ⎥ (5 58)
∂t ∂x ∂y ⎣ ∂y ⎦
avec :
Zs (x,y,t) , Zinf (x,y) : côtes de la surface libre Zs (m) et du fond Zinf (m), par rapport à un
repère fixe tel que le niveau de la mer (axe Oz vertical, vers le haut);
η(x,y,t) : tirant d'eau (profondeur, water depth) :
η(x,y,t) ≈ Zs (x,y,t) – Zinf (x,y)
Ux (x,y,t), Uy (x,y,t) : composantes horizontales suivant Ox et Oy de la vitesse U (m/s),
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 82
variation lente, en gardant tout de même à l’esprit que le passage Saint Venant à l’onde diffusive fait
intervenir des hypothèses plus fortes en 2D (x,y,t) qu’en 1D (x,t) 4.
∂Z s U ∂Z s U
= −c xx (η ) U x et = −c yy (η ) U y (5 65) et (5 66)
∂x η ∂y η
Ces expressions, une fois mises au carré, et additionnées, nous permettent d'exprimer le module
de la vitesse moyennée verticalement :
1/ 4
⎡ 1 ⎛ ∂Z ⎞ 2 1 ⎛ ∂Z s ⎞ ⎤
2
U =η ⎢ 2 ⎜ ⎟ + 2 ⎜ ⎟ ⎥
1/ 2 s
c yy ( η) ⎜⎝ ∂y ⎟⎠ ⎥⎦
(5 67)
⎢⎣ c xx ( η) ⎝ ∂x ⎠
Ceci, avec les équations (5.64) et (5 65), nous permet de ré-exprimer les vitesses moyennées Ux et
Uy en fonction du tirant d'eau η et des pentes locales de la surface libre :
−1 ∂Z s
Ux = η1/ 2
⎡ 1 ⎛ ∂Z ⎞ 2
1/ 4
∂x
1 ⎛ ∂Zs ⎞ ⎤
2
(5 68)
c xx ( η) ⎢ 2 ⎜ s
⎟ + 2 ⎜⎜ ⎟⎟ ⎥
⎢⎣ c xx ( η) ⎝ ∂x ⎠ c yy ( η) ⎝ ∂y ⎠ ⎥⎦
−1 ∂Z s
Uy = η1 / 2
⎡ 1 ⎛ ∂Z ⎞ 2
1/ 4
∂y
1 ⎛ ∂Z s ⎞ ⎤
2
(5 69)
c yy ( η) ⎢ 2 ⎜ s
⎟ + 2 ⎜ ⎟ ⎥
η ∂
⎣⎢ c xx ( ) ⎝ x ⎠ c yy ( η) ⎜⎝ ∂y ⎟⎠ ⎥
⎦
avec (5.14) et (5.15), l'équation de conservation de masse devient :
⎛ 2 − 1/ 4 ⎞
∂Z s ∂ ⎜ 1/ 2 ⎡ 1 ⎛ ∂Z ⎞ 2 1 ⎛ ∂Z s ⎞ ⎤ ∂Z s ⎟
= ⎜η c xx (η ) ⎢ 2
−1
⎜ s
⎟ + ⎜ ⎟ ⎥ η °
∂t ∂x ⎜ ⎢⎣ c xx (η ) ⎝ ∂x ⎠
2
c yy (η ) ⎜⎝ ∂y ⎟⎠ ⎥⎦ ∂x ⎟⎟
⎝ ⎠
(5 70)
⎛ − 1 / 4
⎞
⎡ 1 ⎛ ∂Z ⎞ 2 1 ⎛ ∂Z s ⎞ ⎤
2
∂ ⎜ 1/ 2 ∂Z s ⎟
+ ⎜η c yy (η ) ⎢ 2
−1
⎜ s
⎟ + ⎜ ⎟ ⎥ η °
∂y ⎜ ⎢⎣ c xx (η ) ⎝ ∂x ⎠
2
c yy (η ) ⎜⎝ ∂y ⎟⎠ ⎥⎦ ∂y ⎟⎟
⎝ ⎠
Cette dernière équation de conservation de masse s'écrit aussi :
4 R. Ababou (circa 2000) montre que l’équation obtenue en négligeant les termes d’accélération dans
Saint Venant 2D, est différente de la version 2D de l’équation d’onde diffusive, même si elle en est
qualitativement proche. Voir par exemple le rapport du projet « EiCRIN » (2000) sur la propagation des crues
rapides et les risques d’inondation.
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 84
∂Z s ∂ ⎛ ∂Z ⎞ ∂ ⎛ ∂Z ⎞
= ⎜ Tˆxx s ⎟ + ⎜ Tˆyy s ⎟ (5 71)
∂t ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠
η1/ 2
Tˆxx = 1/ 4
η
⎡ 1 ⎛ ∂Z ⎞2 1 ⎛ ∂Zs ⎞ ⎤
2
(5 72)
cxx (η)⎢ 2 ⎜ s
⎟ + 2 ⎜⎜ ⎟⎟ ⎥
⎣⎢ cxx (η) ⎝ ∂x ⎠ cyy (η) ⎝ ∂y ⎠ ⎦⎥
η 1/ 2
Tˆyy = 1/ 4
η
⎡ 1 ⎛ ∂Z ⎞2 1 ⎛ ∂Z s ⎞ ⎤
2
(5 73)
c yy (η )⎢ 2 ⎜ s
⎟ + 2 ⎜⎜ ⎟⎟ ⎥
⎣⎢ cxx (η ) ⎝ ∂x ⎠ c yy (η ) ⎝ ∂y ⎠ ⎦⎥
L’équation finalement obtenue est l'équation d’onde diffusive 2D dans le cas le plus général
(anisotrope), exprimée dans un repère horizontal. Cette forme de l'équation d’onde diffusive,
associée à la formule de Manning comme loi de perte de charge, a été notamment proposée et
résolue numériquement par Hromadka et al. (1985) et Di Giammarco et al. (1996).
η 1/ 2
Tˆ = 1/ 4
η
⎡ 2 ⎛ ∂Z ⎞ 2
⎛ ∂Z ⎞
2
⎤ (5 74)
⎢c (η )⎜ s ⎟ + c (η )⎜⎜ s ⎟⎟ ⎥
2
⎢⎣ ⎝ ∂x ⎠ ⎝ ∂y ⎠ ⎥⎦
Ceci peut encore s’écrire sous la forme Tˆ = Kˆ η avec :
η1/ 2 η1/ 2
Kˆ = =
c(η )1 / 2 ∇Z S
1/ 4 1/ 2
⎡⎛ ∂Z ⎞ 2 ⎛ ∂Z ⎞ 2 ⎤ (5 75)
c(η ) 1/ 2
⎢⎜ s ⎟ + ⎜⎜ s ⎟⎟ ⎥
⎢⎣⎝ ∂x ⎠ ⎝ ∂y ⎠ ⎥⎦
Dans le cas d’une loi de frottement de type Manning, avec cii (η)=(Maii )2 η -1/3, on obtient :
η 5/3 η2/3
Tˆ = Kˆ η = et Kˆ = (5 76)
Ma ∇Z S
1/ 2
Ma ∇Z S
1/ 2
~ 2K S
Loi de perte de charge : K ( h, ∇H ) = K R ( h)
δ + (δ 2 + 4γK ( h )3 / 2 ∇H ) 1/ 2 (5 78)
avec:
δ = 1 : pertes de charge linéaires (Darcy classique) ;
γ = C/(gν)1/2 ; C = 0.55 (constante adimensionnelle d’Ergun) ;
K (h ) = K S K R (h ) .
Cette formulation est, essentiellement, la formulation générique du code BigFlow. Elle provient
d’une reformulation de la loi généralisée linéaire/quadratique combinant additivement Darcy et
Forchheimer. Pour obtenir l’équation d’onde diffusante 2D / Manning, il faut utiliser la table de
correspondance suivante :
K ( h, ∇H ) = Tˆ (η , ∇Z s )
~ conductance non linéaire transmittance non linéaire
5 Ceci est nécessaire pour implémenter les formules classiques de type Chézy, Manning, etc. Cependant, en
réalité, pour les écoulements très faibles, rampants et/ou à très faible tirant d’eau, il y a un régime de pertes de
charges linéaires de type Darcy-Stokes (Ababou, 1998).
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 86
1/ 4
KS
Tˆ (η , ∇Z S ) = K R ( h) (5 79)
∇Z S
1/ 2
η 5/3
Tˆ (η , ∇Z S ) = (5 80)
Ma ∇Z S
1/ 2
Pour obtenir identité, on voit qu’il suffit de prendre dans les coefficients de BIGFLOW :
K S = Ma −4 et K R (η ) = η 5 / 3 .
Par ailleurs, une autre formulation de cette identification est possible en ré-écrivant la loi de
Manning sous la forme :
c(η ) = C η α avec C = Μa 2 et α = −1 / 3 .
On voit alors qu’on peut identifier les coefficients génériques de BIGFLOW comme suit en
fonction des coefficients de perte de charge C et α :
K S = 1 / C 2 et K R (η ) = η (3−α ) / 2 (5 81)
Plus généralement, la relation d’équivalence ci-dessus devrait être valable pour toute loi de perte
de charge :
Manning : Cii = Ma-2 et α = -1/3 ;
Chézy : Cii = Cz2 et α=0;
Darcy-Weisbach : Cii = λ/8g et α = 0.
• ZS (x,y): c’est la surface libre de l’eau ; elle représente la charge hydraulique et donc
l’inconnue du problème ; dans le sous sol elle est associée au niveau de la nappe et en
surface elle donne le niveau de la rivière ; elle est continue à l’interface des deux régions.
• η(x,y): on peut aussi définir le tirant d’eau (ηr) dans la rivière et dans l’aquifère (ηa).
Le modèle à bi-couche est utilisé pour coupler les écoulements 2D au sein de deux couches
superposées : une couche dite inférieure, comprise entre deux surfaces de cotes respectives
Zinf (x,y) et Zsup (x,y), et une couche dite supérieure reposant sur la surface de cote Zsup (x,y). La
nappe est comprise entre la surface inférieure de cote Zinf (x,y) sur laquelle elle repose, et l'unique
surface libre de cote Zs (x,y,t), située soit dans la couche inférieure, soit dans la couche supérieure.
La zone saturée délimitée par le substratum Zinf et la surface libre Zs, est donc continue, comme
pour le modèle à une couche.
Les écoulements au sein de chaque couche peuvent être de même nature (e.g., régis dans chaque
couche par l'équation de Boussinesq), ou de nature différente (e.g., équation de Boussinesq dans une
couche, et onde diffusive dans l’autre). Ainsi, plusieurs configurations d'écoulements sont possibles :
• écoulements dans un aquifère stratifié horizontalement (bi-couche);
• écoulements dans une plaine alluviale constituée d'un cours d'eau connecté à sa nappe
d'accompagnement (figure 5.1) ;
• affleurement de nappe en un point bas d'un versant (dans un modèle comme MODFLOW,
ces affleurements sont modélisés par une limite à potentiel imposé, là où le taux
d'infiltration de la pluie est supérieur aux possibilités de stockage de la nappe) ;
• écoulements dans un réseau de rivières souterraines (couche inférieure), alimentées par le
drainage d'un massif poreux sus-jacent (couche supérieure).
Zs (x,y,t)
ηr
η = ηa
η ηa
x
x=L Zone B
Zone B Zone A
où η est le tirant d’eau total. En remplaçant ηa et ηr dans l’eq (5 84), on obtient donc :
∂
[θ r η r + θ a η a ] = div(Tr ∇η + Tr ∇Z inf + Ta ∇Z sup ) (5 86)
∂t
∂θ a η
= div(Ta ∇η + Ta ∇Z inf ) (5 89)
∂t
Trois termes de gradients interviennent dans cette équation : un terme de gradient de charge
hydraulique, un terme de gradient gravitaire qui dépend de la pente du substratum, et un terme de
gradient gravitaire qui dépend de la pente du fond de rivière. Les coefficients T’ et T” traduisent les
transmissivités effectives de l’aquifère, de la rivière, et/ou de l’ensemble du système nappe-rivière :
Tˆ ′ = ωTa + λTr
(5 91)
Tˆ ′′ = λ (T + ωT )
a r
où pour ηr = 0 Î {ω = 1 et λ = 0}
pour ηr ≠ 0 Î {ω = 0 et λ = 1}
On peut enfin vérifier la continuité de pression ou de charge à l’interface des zones A et B. On se
place au niveau de l’abscisse x = L entre la zone A et la zone B (voir figure 5.1) :
⎧η → η a
⎪
pour x → L : ⎨η a = Z sup − Z inf
⎪
⎩η r → 0
En insérant ceci dans l’équation (5 90) pour x → L − on obtient :
∂θ a η
= div(Ta ∇η + Ta ∇Z inf ) (5 92)
∂t
Et en insérant ces mêmes résultats dans l’équation (5 90) pour x → L + on obtient :
∂θ a η
= div(Ta ∇η + Ta ∇Z inf ) (5 93)
∂t
Ainsi la continuité de la charge hydraulique est bien assurée entre les zones A et B.
Z SALT = Z SEA −
(H − Z SEA )
− ΔZ (5 94)
ε
où : Zsea est le niveau de la mer ;
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 90
Cette modification revient à remplacer le paramètre Zinf(x,y) par l’inconnue Zsalt(x,y,t) dans la
zone du biseau salé.
Nous simulons avec BIGFLOW l'écoulement permanent 1D dans un canal horizontal, de section
uniforme, de coefficient de rugosité de Manning Ma uniforme, d'extension longitudinale L et limité
en amont et en aval par deux réservoirs de cotes respectives Zs0 et ZsL . Les valeurs des paramétres
utilisés sont :
Ma = 0.015 m -1/3 .s
L = 100 m
Zs0 = 1.5 m
ZsL = 0.5 m
Pour ce problème d'écoulement permanent 1D, l'équation qui régit l’écoulement dans le canal
s'écrit :
⎡ ⎤
⎢ ⎥
∂ ⎢ 1 / Ma 2
5 / 3 ∂Z s
⎥
⎢ η ⎥=0 (5 96)
∂x ⎢ ⎡ 1 ⎛ ∂Z ⎞ 2 ⎤ 1/ 4 ∂x ⎥
⎢⎢ 4 ⎜ ⎟ ⎥ ⎥
s
∂
⎢⎣ ⎢⎣ Ma ⎝ x ⎠ ⎥⎦ ⎥⎦
Soit en supposant le gradient (5.97) de même signe sur toute la longueur du canal.
∂Z s ∂η ∂Z inf
= + (5 97)
∂x ∂x ∂x
et
1
∂Z s ∂x η 5 / 3 = C (5 98)
Ma
d'où :
∂η ∂Z
= Cη −10 / 3 − inf (5 99)
∂x ∂x
Cette équation différentielle admet une solution simple pour un canal horizontal,
(∂Zinf /∂x = 0 → η = Zs ), et avec les conditions limites Zs (0) = Zs0 et Zs (L) = ZsL . Nous obtenons,
avec ZsL > Zs0 , la solution analytique suivante pour la cote de la surface libre :
3 / 13
⎡ 13 / 3 Z s130 / 3 − Z sL
13 / 3
⎤
Z s ( x) = ⎢Z s0 − x⎥ (5 100)
⎣ L ⎦
Ce qui donne, de plus, le débit spécifique (m2 /s) constant :
1 1 3 Z s130 / 3 − Z sL
13 / 3
Q sx = ∂Z s ∂x η =
5/ 3
(5 101)
Ma Ma 13 L
La figure 5.3 montre la comparaison du profile Zs(x) analytique et du profile Zs(x) numérique
(BigFlow). Le profile est correctement simulé. Le débit spécifique est aussi correctement calculé :
Calcul analytique : Qs-ana ≈ 7.676 m2 s-1 ;
Calcul numérique : Qs-num ≈ 7.678 m2 s-1 ;
L'erreur relative sur le débit est de 0.03 %.
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 92
2.0
1.8
1.6
1.4
1.2
1.0
Zs (m)
0.8
0.6
0.4
0.2
0.0
0 10 20 30 40 50 60 70 80 90 100
x (m)
Fig. 5.3 Position de la surface libre Zs (x) dans le cas d'un canal uniforme sur un fond
horizontal. Comparaison entre solutions numérique ( ο ) et analytique ( - ) (d’après Trégarot
2000).
Cette dernière hypothèse implique un écoulement quasi-horizontal, ce qui peut se justifier pour des
pentes du fond (z = Zsup) et de la surface libre (z = Zs) relativement faibles. On obtient les équations
de Saint-Venant (voir plus haut), avec :
En ignorant les termes d'accélération locale et convective, on obtient les équations d’ondes
diffusantes (voir plus haut et ci-dessous) :
∂( ηU ) ∂ ( ηV ) ∂Z s
+ + =0
∂x ∂y ∂t
∂Z s
+ Sfx = 0
∂x
∂Z s
+ Sfy = 0
∂y
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 93
n 2x (U 2 + V 2 ) n 2y (U 2 + V 2 )
1/ 2 1/ 2
avec Sfx = 1 / 3 U et S fy = V
η η η 1/ 3 η
où nx et ny sont les coefficients de frottement de Manning (en m-1/3.s) dans les directions x et y,
supposées directions principales d'anisotropie pour ces coefficients. Ce sont ces dernières équations
qui sont résolues par BIGFLOW.
Application
Nous simulons avec BIGFLOW l'écoulement 1D dans un canal de pente I = -0.001, de section et
de coefficient de rugosité n = 0.001 uniformes, d'extension longitudinale L = 2000 m et limité en
amont et en aval par deux réservoirs de hauteurs respectives η0 = 10 m et ηL = 5 m.
Initialement, un "barrage" partage pratiquement le domaine en deux (situé en X = 1025 m): le niveau
est à 10 m en amont, et de 5 m en aval. Pour t > 0, le "barrage" est "rompu". Nous avons réalisé la
simulation jusqu'à t = 3 s. Les temps de sortie sont 0.01, 0.05, 0.1, 0.5, 1 et 3 unités de temps (s).
La taille des mailles choisie est de 12.5 m (domaine de 161 noeuds selon la direction X).
Voici pour exemple, la structure du fichier d'entrée principal utilisé par BIGFLOW 2000 :
FICHIER INPUT1 DU TEST « RUPTURE DIGUE 1D » (BIGFLOW 2000)
IDATE (DAY-MONTH-YEAR):
091098
IDRUN (SIMULATION ID NUMBER):
109801
LRUN (0:TEST RUN/1:FULL RUN):
1
A A’
B’
Vue plan du domaine de calcul Section BB’ – coupe longitudinale du domaine de calcul
Fig. 5.7 Simulation transitoire simplifiée rivière/aquifère avec le modèle couplé 2D plan de
BIGFLOW (méandre rectangulaire sans pente) (t=45s).
Fig. 5.8 Simulation transitoire simplifiée rivière/aquifère avec le modéle couplé 2D plan de
BIGFLOW (méandre rectangulaire sans pente) (t=100s).
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 99
Fig. 5.11 Représentation des données disponibles dans la base de données SIG
de la Garonne au niveau de Montbéqui.
V - 5.2.1 Krigeage
Le krigeage consiste à calculer les coefficients Wi de la fonction F à l’aide des valeurs du
variogramme. L’hypothèse principale indispensable pour utiliser le krigeage est que la moyenne μ et
la variance Var de la fonction F soient stationnaires, donc qu’elles ne dépendent pas de la position
des points, mais seulement de la distance entre les points.
μ ( x) = E{Z ( x)} (5 103)
Var{Z ( x )} = E{ [ Z ( x ) − μ( x )] 2 } (5 104)
Cette hypothèse est applicable dans le cas où les structures qu’on mesure présentent une forme
assez régulière.
1 n(h)
γ ( h) = ∑
2n(h) i =1
( xi − y i ) 2 (5 105)
A l’aide de la méthode des moindres carrés, une fonction est ajustée sur les valeurs
expérimentales du semi-variogramme. On obtient ainsi la fonction du semi-variogramme. Le choix
et l’ajustement de la fonction au semi variogramme constituent la partie délicate du krigeage
(Gratton, 2002).
comme si c’étaient de nouvelles données, alors les krigeages subséquents ne s’en trouvent
pas modifiés (sauf pour la variance de krigeage).
L’orientation de chaque profile par rapport à l’horizontal est obtenue en interpolant linéairement
l’angle β suivant l’axe τ :
di
β i = β1 + (β 2 − β 1 ) (5 110)
d
où β 1,2 est l’angle entre l’horizontale et le profile mesuré P1,P2 définis.
par tan( β 1, 2 ) =
(y B 1, 2 − y A1,2 )
(x B 1, 2 − x A1,2 )
avec A et B sont les points extrêmes de chaque profile. La position
sur l’axe n de chaque point des profiles interpolés est donnée par :
−1
⎛π ⎞
τi j = τc i = d i cos( α i ) et ni j = Pi j cos ⎜ − α i − β i ⎟ (5 111) et (5 112)
⎝2 ⎠
avec lij = ( xij − xci ) 2 + ( yij − yci )2 la distance entre le centre et un point du profile
L’interpolation entre deux lignes sera effectuée de façons linéaires entre deux profiles. La distance
∑
N
entre les deux profiles est donné par d ≈ 1
d i . Plusieurs pairs de points sont considérés entre les
deux profiles. Les paires de points suivant sont relié automatiquement entre deux profiles :
• les points de centre ;
• la rive gauche du lit majeur ;
• la rive gauche du lit mineur ;
• la rive droite du lit majeur ;
• la rive droite du lit mineur.
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 104
L’élévation du fond de rivière sur tous les profiles intermédiaires et pour chaque paire de points
est calculée avec la formule suivante :
Z i j = Z1 j +
di
d
(
Z 2 j − Z1 j ) (5 113)
Le calcul des coordonnées (x,y) à partir des coordonnées (t,n) pour les points interpolés se fait à
partir des formules suivantes :
X i j = Xc i + lij cos( β i ) = X i −1 + d i cos( α i ) + l ij cos( β i ) (5 114)
A2
P2
Xi
Ci
di B
τ P
γ n
α1
y A
β1
P1
B
x
Fig. 5.12 Transformation des coordonnées cartésiennes en coordonnées curvilignes
6. Interpolation intégrée sur la zone complète : finalement les données regroupées dans
l’étape (1) sont associées aux données obtenues par l’interpolation à l’étape 5. Une
dernière interpolation permet d’obtenir le MINT (voir figure 5.13).
Transformation Assemblage
curviligne Transformation des données Triangulation
P3 P3 P3
P3
P2
P2 P2 P2
τ
n
P1 P1
P1
y P1 τ y y
n
x x x
Le lecteur est aiguillé vers Weng et al. (2003) concernant les données disponibles sur le site, et
pour les résultats de modélisation de l’eau souterraine (non couplé) et pour la calibration effectuée.
La calibration de la distribution de la conductivité hydraulique effectuée dans Weng et al. (2003) est
utilisée dans cette étude.
On remarque aussi que la taille du domaine d’écoulement de l’eau dans la rivière évolue de
façonsdynamique dans le temps. Cela est perceptible en comparant l’eau au niveau du bancs de
galets entre les figures 5.16 à 5.17. Cette caractéristique rend le modèle spécialement adapté pour la
modélisation des processus d’inondations lentes.
Dans la figure 5.20 les niveaux d’eau dans les piézomètres en régime permanent sont comparée
aux niveau d’eau moyen annuel mesuré dans les piézomètres (1 à 27). On remarque une bonne
correspondance entre les résultats sauf pour les piézomètres (23 et 26). Il est difficile d’expliquer la
différence au niveau des piézomètres 23 et 26, mais on suppose que les paramètres de conductivité
ne sont pas bien estimés.
V - 7 CONCLUSIONS
Dans ce chapitre, on s’est intéressé aux interactions « surface/souterrain », en présence ou non de
l’intrusion saline. Les équations d’écoulements plans (x,y) pour les eaux de surface sont développées
avec les différentes simplifications possibles.
Ensuite une analogie est montrée entre les équations d’écoulement d’eau de surface et les
équations d’écoulement souterrain intégré verticalement (Dupuit-Boussinesq). On montre l’analogie
avec la formulation plus générale des écoulements en milieux macroporeux à perte de charge
quadratique (Ward-Forchheimer).
Puis une procédure de couplage surface/souterrain est développée. La procédure consiste à
superposer les écoulements dans la rivière et l’aquifère, en additionnant les équations de
conservation de masse dans chaque couche.
Cette procédure de couplage est généralisée au cas de l’intrusion saline en utilisant le modèle
d’intrusion saline 2D plan proposé dans le chapitre 4 pour remplacer le module Dupuit-Boussinesq
de la couche sol.
Le modèle d’écoulement de surface est validé avec la solution analytique de l’équation d’onde
diffusive en 1D en permanent et en transitoire. Les tests montrent une bonne correspondance entre
les profiles de la ligne d’eau et entre les flux calculées par le modèle numérique et la solution
analytique.
Il n’existe pas de solution analytique pour le modèle couplé surface/souterrain. Celui là a été testé
qualitativement par deux cas d’écoulement dans des géométries simplifiées. Le premier test est un
domaine avec un canal à pente uniforme et un substratum à pente uniforme aussi. La dynamique du
test montre le couplage entre la nappe et la rivière correctement. Le deuxième teste est un méandre
Modélisation 2D couplée surface / souterrain avec ou sans intrusion saline 111
en U. Dans ce cas on considère le modèle couplé rivière aquifère avec intrusion saline. Dans ce cas
aussi la dynamique semble être valide.
Enfin la méthode est appliquée à la modélisation de l’écoulement dans une partie de la vallée
alluviale de la Garonne au niveau de Monbéqui, entre les villes de Toulouse et Moissac. Pour
effectuer cette simulation, un modèle numérique intégrant le fond de la rivière et la topographie était
indispensable. Une méthode pour l’obtention d’un modèle intégré numérique de terrain (MINT) à
partir d’un nombre limité de sections de rivières est également élaborée et appliquée au site
« Garonne ». Les résultats indiquent une bonne correspondance entre les résultats du modèle en
permanent et les moyennes annuelles des mesures piézoélectriques.
Conclusions 112
Chapitre VI CONCLUSIONS
Conclusions 113
dynamique entre l’eau douce et l’eau salée en chaque point du domaine. On peut ajouter que le
modèle donne une estimation exacte de l’intersection de l’interface avec le substratum.
L’implémentation du modèle à interface abrupte en 3D avec deux fluides en mouvement sera un
atout majeur vers la modélisation des aquifères côtiers avec de forts taux de pompages. Ce modèle
aura l’avantage de ne pas ajouter de la diffusion numérique. En pratique les modélisateurs
surévaluent la diffusion pour des raisons de stabilité numérique lors de l’utilisation des modèles avec
zone de mélange, dans les cas d’études à grandes échelles. Par contre il sera difficile d’utiliser le
modèle a interface abrupte dans le cas des aquifères qui sont sujets à la dissolution ou à la
précipitation de sels. Enfin le modèle à interface abrupte n’est pas adapté au problème à fort
contraste de densité ou des problèmes à double convection qui font intervenir de la
thermoconvection.
♦ Modélisation stochastique analytique et numérique
Un autre aspect étudié au cours de cette thèse est l’effet de l’hétérogénéité du milieu poreux sur la
position de l’interface entre l’eau salée et l’eau douce.
Plusieurs simulations numériques sont effectuées sur des milieux moyennement jusqu'à fortement
hétérogènes en perméabilité. Ces résultats sont confrontés à une solution analytique stochastique
obtenue par analyse spectrale Fourier. La position de l’interface, qui est la variable d’intérêt, est
statistiquement non-stationnaire. Pour appliquer la méthode spectrale/Fourier, on applique une
transformation quadratique de la position de l’interface. L’étude spectrale Fourier est ainsi faire sur
la transformée (Φ) statistiquement stationnaire obtenue. Enfin, on récupère les moments de la
variable d’origine (la position de l’interface) à partir de ceux de la variable transformée (Φ). La
comparaison entre les résultats de la solution analytique stochastique et les simulations numériques
montre une bonne concordance loin de la mer.
Les résultats montrent aussi que l’augmentation de l’hétérogénéité en 2D plan de la perméabilité
accentue le phénomène de l’intrusion saline en moyenne et l’incertitude de la position. L’analyse de
l’incertitude de la position de l’interface est de l’ordre de plusieurs mètres pour un aquifère, de
profondeur 30 m et de longueur 1000 m.
La solution analytique stochastique obtenue peut être utilisée pour une première estimation de
l’incertitude de la position de l’interface eau douce / eau salée dans le cas d’une étude
hydrogéologique préliminaire sur l’intrusion saline. Les données nécessaires pour la formule sont :
la description de la structure d’hétérogénéité du milieu poreux à travers la moyenne, l’écart type et
les longueurs de corrélation. Celles-ci peuvent être associées à des longueurs caractéristiques des
structures d’hétérogénéité observées sur site.
Une future étude stochastique numérique et perturbative avec des modèles numériques et
mathématiques 3D permettrait d’obtenir des résultats plus réalistes. On peut s’attendre à ce que ces
modèles 3D donnent une incertitude plus faible de l’interface, vu que les chemins préférentiels
seront plus nombreux que dans le cas d’une configuration 2D.
Une prise en compte des conditions limites du biseau salé à travers une approche stochastique
prenant en compte la non-stationnarité due aux conditions limites, permettrait de se rapprocher des
cas réels.
♦ Couplage surface / souterrain en écoulements plans
Une méthode innovante est utilisée pour la modélisation couplée d’un système rivière/aquifère en
écoulement plan.
Le modèle développé considère une seule équation générique pour les deux domaines de calcul.
Cette équation exprime un écoulement de surface régi par les équations d’onde diffusive, et un
écoulement en milieu poreux régi par les équations de Boussinesq. Les résultats obtenus montrent
les échanges entre les deux systèmes, surtout au niveau des structures morphologiques complexes
des rivières, telles que les méandres.
Conclusions 115
RÉFÉRENCES
Ababou R. Random Porous Media Flow on Large 3D Grids: Numerics, Performance, and
Application to Homogenization, Chap.1, pp.1-25. In: IMA Vol 79 Mathematics and its Applications:
Environmental Studies (Math. Comput. Statist. Anal.). Wheeler MF (ed.), Springer, NY, 410 pp.
1996.
Ababou R., Al-Bitar A., Peyrard D., Quintard M., Sanchez-Perez J.M., Sauvage S., Vervier P.,
Weng P. Modeling coupled surface / subsurface flow interactions : implementation and comparison
of three models based on Darcy, Boussinesq / Saint Venant, and Boussinesq / diffusive wave, with
application to the Garonne floodplain, Midi-Pyrénées, France. Proceedings IAHR-GW 2006 on
“Groundwater Hydraulics in Complex Environments”, Toulouse, France, 12-14 June 2006.
Ababou R., Trégarot G. Coupled Modeling of Partially Saturated Flows : Macro-Porous Media,
Interfaces, and Variability. Proceedings CMWR 2002, Computational Methods in Water Resources,
23-28 June 2002, Delft, The Netherlands, 8pp.
Ababou R., Trégarot G., Larabi A. Partially Saturated Hydrological Flows : Numerical Experiments
and Analyses. XII Internat. Conf. on Computational Methods in Water Resources (CMWR)
Proceedings, The Chersonese, Crete, Greece, 8 pp., 15-19 June 1998.
Ababou R., Trégarot G., Bouzelboudjen M. Variably Saturated Subsurface Flow with Layers and
Interfaces : Perched Water Tables and Stream-Aquifer Connection. ModelCare'96 Proc.,
International Groundwater Modeling Center (IGWC), GWMI Series No.96-OX, Colorado School of
Mines, 10 pp., 25-27 Sept. 1996.
Ababou R. et Bagtzoglou A.C. BIGFLOW: A Numerical Code for Simulating Flow in Variably
Saturated, Heterogeneous Geologic Media (Theory & User's Manual, Version 1.1), Report
NUREG/CR-6028, U.S. Nuclear Regulatory Commission, Government Printing Office, Washington
D.C., U.S.A., 139 pp. 1993.
Ababou R., Sagar B., and Wittmeyer G. Testing Procedures for Spatially Distributed Flow Models.
Advances in Water Resources, Vol.15, pp. 181-198 1992.
REFERENCES 118
Ababou R., Gelhar L.W., et McLaughlin D., 1988. Three-Dimensional Flow in Random Porous
Media, Report No. 318, Ralph Parsons Laboratory for Water Resources & Hydrodynamics,
Massachusetts Institute of Tech., Cambridge, MA 02139, March 1988 (2 vols., 833 pp.).
Abarca E., Carrera J., Sanchez-Vila. Anisotropic Dispersive Henry Problem., Advances in Water
Ressources, doi :10.1016/[Link].2006.08.005, 2006.
Ackerer Ph., Younes A., Mosé R. Modelling variable density flow and solute transport in porous
medium : 1. Numerical model and verification, Transport in Porous Media, Vol. 35, pp.345-373,
1999.
Aharmouch A. Intrusion marine dans les aquifères côtiers : développement et tests de modèles
numériques à éléments finis 3D, Thèse de doctorat Ecole Mohammadia d'Ingénieurs, Rabat, Maroc.
2003.
Albitar A. & Ababou R. Random Field Approach to Seawater Intrusion in Heterogeneous Coastal
Aquifers:Unconditional Simulations and Statistical Analysis, Chapter in GeoENV: Geostatistics for
Environmental Applications, Renard P., Demougeot-Renard H., Froidevaux R. (Eds.), ISBN: 3-540-
26533-3, Springer, 2005.
Al-Bitar A. & Ababou R. Modeling of Salt Water Intrusion In Coastal Aquifers: Studying the effect
of Heterogeneity and Uncertainties towards seawater intrusion remediation, Chapter in “Monitoring
Modeling and Management of Coastal aquifers”. Benavente J., Larabi A., El Mabrouki K. (eds.),
pp.391-409, Water Research Institute, University of Granada, Spain, 2004. (ISBN: 84-699-7903-5).
Al-Bitar A., Ababou R. Studying the effect of Heterogeneity and Uncertainties of Salt Water
Intrusion In Coastal Aquifers. Proceedings of the 1st International Conference on Thermal
Engineering Theory and Applications (ICTEA), Beyrouth, Lebanon, 31 May–4 June 2004. Proc.
paper ICTEA-ES1-05 (29 pp.), 2004.
Al-Bitar A. & Ababou R., Coupled surface / groundwater flow via a Diffusive wave / Boussinesq
model and undersharp interface saltwater intrusion conditions. Extended abstract in Pre-Proceedings
IAHR-GW2006 International Conference on “Groundwater Hydraulics in Complex Environments”,
Toulouse, France, June 12-14, 2006.
Aharmouch, A., Larabi A., et Hilali M. A 3D model for groundwater flow and seawater intrusion
interface: Application to the Martil coastal aquifer system. In Proc XIV Int. Conf. on Computational
Methods in Water Resources, Vol, 1, Density Dependent Flow and Transport, ed. S.M.
Hassanizadeh, R.J. Schotting, W.G. Gray and G.F. Pinder. Elsevier, pp. 539-546, Delft, The
Netherlands, 2002.
Badon-Ghyben W. Nota in verband met de voorgenomen putboring nabij Amsterdam [Notes on the
probable results of well drilling near Amsterdam] In: Tijdschrift van het Kononklijk Instituut van
Ingenieurs The Hague, vol. 9. 1888. p. 8–22.
Bakker, M., Oude Essink, G.H.P., Langevin C.D. The rotating movement of three immiscible
fluids—a benchmark problem. Journal of Hydrology 287, 270–278, 2004.
REFERENCES 119
Bear J., Cheng A.H.-D, Sorek S., Quazar D., and Herera I., Editors, Seawater Intrusion into Coastal
Aquifers-Concepts, Methods, and Practices, Kluwer Academic Publications, 625pp., 1999.
Bedient P. B., Huber C. W. Hydrology and Floodplain Analysis (3rd ed.), 2002.
Bouzouf B., Gloth O., Hanel D., Vilsmeier R., Simulation of discontinuous flows in porous media
First international conference on Saltwater Intrusion and Costal Aquifers Monitoring, Modeling and
Management. Essaouira, Morocco, April 23-25, 2001
Boussineq J., Recherches théoriques sur l'écoulement des nappes d'eau infiltrées dans le sol et sur le
débit des sources, C.R.H Acad. J. Math. Pures Appliquées 10 5–78, 1904.
Brooks R.H., Corey A.T. Hydraulic Properties of Porous Media. Hydrology paper 3, Colorado State
University, Fort Collins, 1964.
Brutsaert W. Probability laws for pore-size distributions. Soil Sci. 101:85.92, 1966.
Buckingham E. Studies on the Movement of Soil Moisture. Bulletin No. 8, Washington, DC, US,
Department of Agriculture, Bureau of Soils, 1907.
Burdine N.T. Relative Permeability Calculations From Pore Size Distribution Data. Petroleum
Transactions, AIME, Vol. 198, pp. 71-78, 1953.
Campbell G.S. A Simple Method for Determining Unsaturated Conductivity from Moisture
Retention Data. Soil Science, Vol. 117, No. 6, pp. 311-314, 1974.
Celia M.A., Bouloutas E.T., Zarba E.L. A General Mass-Conservative Numerical Solution for the
Unsaturated Flow Equation, Water Resour. Res., 26, 1483-1496, 1990
Cheng, A.H.-D. and Ouazar, D. (eds.), Coastal Aquifer Management-Monitoring, Modeling, and
Case Studies, Lewis Publishers, 280 p., 2003.
Childs E.C., Collis-George N. The Permeability of Porous Materials. Proceedings of the Royal
Society of London, Series A : Mathematical and Physical Sciences, Vol. 201, pp. 392-405, 1950.
Chow V.T., Maidment D.R., Mays L.W. Applied Hydrology. Mc Graw-Hill International Editions,
Civil Engineering Series, 572 pp. 1988.
Croney D. The movement and distribution of water in soils. Géotechnique, 3 (1), p. 1-16 (1952).
REFERENCES 120
Darcy H.P.G. Les Fontaines Publiques de la Ville de Dijon, Exposition et Application des Principes
à Suivre et des Formules à Employer dans les Questions de Distribution d'Eau. Paris, France, Victor
Dalmont Ed., 1856.
Delhomme, J.P. Spatial variability and uncertainties in groundwater flow parameters : a statistical
approach. Water Resources Research, 152: 269-280, 1979
de Marsily G. Quelques réflexions sur l’utilisation des modèles en hydrologie Rev Sci Eau
7(3) :219-234, 1994.
Diersch H.J.G., Kolditz O. Variable-density flow and transport in porous media: approaches and
challenges Advances in Water Resources 25, 899–944, 2002.
Diersch H.J.G. Finite element modeling of recirculating density driven saltwater intrusion processes
in groundwater. Adv Water Resour, 11(1):25–43, 1988.
Diersch H.J.G., Nillert P. Saltwater intrusion processes in groundwater: novel computer simulations,
field studies and interception techniques. Int Symp Groundwater Monitoring and Management 1987,
Dresden, IAHS Publ no 173 1990:319–29.
Di Giammarco P., Todini E., Lamberti P. A Conservative Finite Elements Approach to Overland
Flow: the Control Volume Finite Element Formulation. Journal of Hydrology, Vol. 175, pp. 267-
291, 1996.
Dupuit J. Etudes Théoriques et Pratiques sur le Mouvement des Eaux dans les Canaux Découverts et
à Travers les Terrains Perméables. 2nde édition, Dunod, Paris, 1863.
Essaid H. I. A multilayered sharp interface model of coupled freshwater saltwater flow in coastal
systems : Model Development and application, Water Resour. Res., 26(7), 1431-1453, 1990a.
Essink G. O. Improving fresh groundwater supply - problems and solutions. Ocean and Coastal
Management 44:429-449, 2001.
Freeze, R. A. et J. A. Cherry Groundwater, Prentice-Hall, New Jersey, 604 pp.,1979.
Freeze R.A. and Witherspoon P.A. Theoretical analysis of regional groundwater flow: analytical and
numerical solutions to the mathematical model: Water Resources Research, v. 2, n. 4, p. 641-656,
1966.
REFERENCES 121
Freeze, R.A. and Cherry J.A. Groundwater. Prentice-Hall, Inc., Englewood Cliffs, NJ, 604 pp.,
1979.
Gambolati G., Pini G. & Zilli G. Numerical comparison of preconditionnings for large sparese finite
element problems, Numerical Methods for Partial Differential Equations, John Wiley & sons, Inc.,
139-157, (1988b).
Gambolati G., Putti M. & Paniconi C. Three-dimensional model of coupled density-dependent flow
and miscible transport in groundwater. In Bear et al. (eds.) Seawater intrusion in coastal aquifers:
Concepts, Methods, and Practices, chapter 10, pp. 315-362, Kluwer Academic, Dordrecht, The
Netherlands, 1999.
Golub G.H. & Van Loan C.F. Matrix computations. 2nd Ed., J. Hopkins Univ. Press, Baltimore,
Maryland, 642 pp, 1989.
Gunduz O. and Aral M.M. River networks and groundwater flow: a simultaneous solution of a
coupled system. J Hydrol ;301:216–34, 2005.
Haverkamp R., Vauclin M., Touma J., Wierenga P.J. A comparison of numerical simulation models
for one dimensional in_ltration, Soil Sci. Soc. Amer. J., 41: 285-294. 1977.
Henry H.R. Effects of dispersion on salt encroachment in coastal aquifers. In: Sea water in coastal
aquifers US Geol Surv Water Supply Pap, 1613-C. p. 70–84, 1964..
Henry,H. R. Interfaces between saltwater and freshwater in coastal aquifers, U.S. Geological Survey
Water-Supply Paper 1613-C, Sea Water in Coastal Aquifers: C35–C70, 1964.
Herbert A.W., Jackson C.P. and Lever D.A. Coupled groundwater flowand solute transport with
fluid density strongly dependent on concentration, Water Resour. Res. 24(10), 1781–1795, 1988.
Herzberg A. Die Wasserversorgung einiger Nordseebader (the water supply of parts of the North Sea
coast in Germany). Z Gasbeleucht Wasserversorg;44:815–9, and 1901;45:842–4., 1901.
Holzbecher E. Modeling of saltwater upconing. In: Wang S, editor. II Int Conf Hydro-Science and
Hydro-Engin Proc vol. 2, Part A, Beijing. p. 858–65, 1995.
Hromadka T.V., Berenbrock C.E., Freckleton J.R., Guymaon G.L. A Two-Dimensional Dam Break
Flood Plain Model. Advances in Water Resources, Vol. 8, No. 1, pp. 7-14, 1985.
REFERENCES 122
Huyakorn P.S., Anderson P.F., Mercer J.W., White Jr H.O. Saltwater intrusion in aquifers:
development and testing of a three-dimensional finite element model. Water Resour Res
1987;23:293–312.
Johannsen K., Kinzelbach W., Oswald S., Wittum G. The saltpool benchmark problem––numerical
simulation of saltwater upcoming in a porous medium. Adv Water Resour 25(3):335–48, 2002.
Knupp P.M., Lage J.L. Generalization of the Forchheimer-Extended Darcy Flow Model to the
Tensor Permeability Case via a Variational Principle. Journal of Fluid Mechanics, Vol. 299, pp. 97-
104, 1995.
Kolditz O, Ratke R, Diersch H.J.G., Zielke W. Coupled groundwater flow and transport: 1.
Verification of variable-density flow and transport models. Adv Water Resour 1998;21:27–46.
Larabi A., and De Smedt F., Numerical solution of a 3D groundwater flow involving free boundaries
by a fixed finite element method, J. of Hydrology, 201, 161-182, 1997.
Lecca, G. Implementation and testing of the CODESA-3D model for density-dependent flow and
transport problems in porous media. CRS4-TECH-REP-00/40, Cagliari, Italy, 2000.
Lecca, G., & P. Cau. Estimating the impact of a discontinuous confining layer on the seawater
intrusion using a stochastic approach: the Oristano coastal aquifer site (Sardinia, Italy). To be
published in the Proceedings of the Third International Conference on Saltwater Intrusion and
Coastal Aquifers—Monitoring, Modeling, and Management (SWICA M3). Merida, Mexico, March
30 – April 2., 2003.
Liggett, J. A., and P. L-F. Liu. The Boundary Integral Equation Method for Porous Media Flow,
George Allen and Unwin, Boston, 1983.
Matheron G. The Intrinsinc Random Functions & Applications. [Link]., 5, 439-468, 1973
Merwade V. M., Maidment D. R., and Goff J. A. Anisotropic considerations while interpolating
river channel bathymetry. Journal of Hydrology, Volume 331, Issues 3-4, Pages 731-741, 2006.
Mizell S.A., Gutjahr A.L., Gelhar L.W. Stochastic Analysis of Spatial Variability in Two-
Dimensional Steady Groundwater Flow Assuming Stationary and Nonstationary Heads. Water
Resour Res 18(4) 1053-1067, 1982.
Mualem Y. A New Model for Predicting the Hydraulic Conductivity of Unsaturated Porous Media.
Water Resources Research, Vol. 12, No. 3, pp. 513-522, 1976.
Mualem Y., Dagan G. Hydraulic Conductivity of Soils : Unified Approach to the Statistical Models.
Soil Sci. Soc. Am. J., Vol. 42, pp. 392-395, 1978.
Muskat, M. The Flow of Homogeneous Fluids through Porous Media, McGraw-Hill Inc., New
York. Reprinted by J. W. Edwards, Ann Arbor, 1937.
REFERENCES 123
Oberbeck A., Ueber die Wärmelaitung der Flüssigkeiten bei Berucksichtigung der Stromung infolge
von Temperaturdifferenzen, Ann. Phys. Chem. 7 271–292, 1879.
Panday S., Huyakorn P.S. A fully coupled physically-based spatially-distributed model for
evaluating surface/subsurface flow. Adv Water Resour 27:361–82, 2004.
Pinder, G. H., & Bredehoeft J. D. Application of the digital computer for aquifer evaluation, Water
Resour. Res.,4, 1069-1093, 1968.
Putti M. and Paniconi C. Picard and Newton linearization for the coupled model of saltwater
intrusion in aquifer. Adv Water Resour 18(3):159–70, 1995.
Putti M. and Paniconi C. Time step and stability control for a coupled model of surface and
subsurface flow. Proceedings of the XV International Conference on Computation Methods in Water
Resources (CMWR XV), vol. 2. Chapel Hill, NC, USA: Elsevier; p. 1391–402, 2004.
Reilly T.E., Goodman A.S. Analysis of saltwater upconing beneath a pumping well. J Hydrol ;
89:169–204, 1987.
Remson, I.,Appel C. A., & Webster R. A. Groundwater models solved by digital computer, [Link].
Div. ASCE, Vol. 91, no. HY3, pp 133-147, 1965.
Richards L.A. Capillary conduction of liquids through porous media, Physics 1, 318-333, 1931.
Sbai, M. A. Modeling three dimensional groundwater flow and transport by hexahedral finite
elements.,Ph.D. Thesis of Free University of Brussels (VUB), 1999.
Sbai. M. A., De Smedt F. & Larabi A. A Generalized Approach for Modeling 3D Transient Free and
Moving Boundaries in Coastal Aquifers. In Proc. First International Conference on Saltwater
Intrusion and Coastal Aquifers- Monitoring, Modeling and Management. Essaouira Morocco, April
23-25, 2001.
Segol G., Pinder G.F., Gray W.G. A Galerkin-finite element technique for calculating the transient
position of the saltwater front. Water Resour Res;11(2):343–7, 1975
Tang D.H. and Pinder G.F. Analysis of mass transport with uncertain physical parameters. Water
Resources Research, 15 (5), 1147-1155, 1979
Touma, J., Vachaud, G., and Parlange, J.-Y. Air and water flow in a sealed ponded vertical soil
column. Soil Sci. 137: 181-187, 1984.
Vanmarcke, E.H. Random Fields: Analysis and Synthesis, The MIT Press, Cambridge,
Massachusetts, 1983.
REFERENCES 124
Voss C. and Souza W. R. Variable density flow and solute transport simulation of regional aquifers
containing a narrow freshwater-saltwater transition zone, Water Resour. Res. 23(10), 1851–66,
1987.
Voss C.I. USGS SUTRA Code––History, practical use, and application in Hawaii. In: Bear J, Cheng
AHD, Sorek S, Quazar D, Herrera I, editors. Seawater intrusion in coastal aquifers. Dordrecht:
Kluwer Publ; p. 249–313, 1999.
Voss C.I., Souza WR. Variable density flow and solute transport simulation of regional aquifers
containing a narrow freshwatersaltwater transition zone. Water Resour Res;26:2097–106, 1987.
Vogel P., Zur Theorie binärer Fluidgemische in porösen Medien, Bundesanstalt für
Geowissenschaften und Rohstoffe, Hannover, 1995.
Volker R.E., Rushton K.R. An assessment of the importance of some parameters for seawater
intrusion and a comparison of dispersive and sharp-interface modeling approaches. J Hydrol
;56:239–50, 1982.
Ward J.C. Turbulent Flow in Porous Media. Journal of the Hydraulics Division, Proceedings of the
ASCE. Vol. HY 5, pp. 1-12, 1964.
Weng P., [Link]. J.M., Sauvage.S., Vervier.P., and Giraud F. Assessment of the quantitative
and qualitative buffer function of an alluvial wetland : hydrological modelling of a large floodplain
(Garonne River, France), Hydrological Processes, Vol.17 (2), pp. 2375-2392, 2003
Wösten J.H.M., Finke P.A., Jansen M.J.W. Comparison of Class and Continuous Pedotransfer
Functions to Generate Soil Hydraulic Characteristics. Geoderma Vol. 66, pp. 227-237, 1995.
Van Genuchten M.T. A Closed-form Equation for Predicting the Hydraulic Conductivity of
Unsaturated Soils. Soil Sci. Soc. Am. J., Vol. 44, pp. 892-898, 1980.
Wu Y.-S., Forsyth P.A., On the selection of primary variables in numerical formulation for modeling
multiphase flow in porous media Journal of Contaminant Hydrology 48,277–304, 2001
Younes A. On modelling the multidimensional coupled fluid flow and heat or mass transport in
porous media. International Journal of Heat and Mass Transfer 46 367–379, 2003
Younes A., Ackerer Ph., Mosé R., Modelling variable density flow and solute transport in porous
medium : 2. Re-evaluation of the salt dome flow problem, Transport in Porous Media, Vol. 35,
pp.375-394, 1999.
Zammit C. Analyse et Evaluation des Paramètres des Caractéristiques Hydrodynamiques des Sols.
Prediction par un Modèle Analytique à Base Physique à partir de Données Texturales. Thèse de
Docteur de l'Université Joseph Fourier - Grenoble 1, 200 pp., 1999.
Zhang D. Stochastic Methods for Flow in Porous Media Coping with Uncertainties Academic Press,
San Diego, CA. 350 pp. ISBN: 0-12-779621-5, 2002.
REFERENCES 125
ANNEXES
Based on the classification and criteria presented in chapter 2, we present here a representative
selection of computer codes geared towards GW flow modelling, with at least some coupling aspects
taken into consideration (first and foremost, the saltwater density coupling, and then, the coupling to
surface flow and other hydrogeological processes).
The list below is meant to complete the list given in Chapter 2 even though it isn’t exhaustive. To
give just one example, the FRAC3DVS code (Ed Sudicky, University Waterloo, Canada) has not
been documented here, although it includes several features of interest in the framework of this
review. Also, it should be noted that some of the codes are commercial, while others are non-
commercial and/or research codes (publicly available).
We present here a representative selection of computer codes geared towards GW flow Modeling,
with at least some coupling aspects taken into consideration: first and foremost saltwater density
coupling, and secondly, the coupling to surface flow and other hydrogeologic processes.
MT3DMS (Modular 3-D Transport model with Multi-Species add-on reaction packages)
Author: C. Zheng et al.
Code status: Publicly available (The Hydrology Group - University of Alabama)
Description: Three types of solutions methods have been implemented in MT3DMS : the standard
finite difference method; the particle-tracking-based Eulerian-Lagrangian methods; and the
higher-order finite-volume TVD method. MT3DMS uses the flow computed from a finite-
difference flow model to model the transport of species assuming that changes in the
concentration field will not affect the flow field significantly. This approximation is not suitable
a priori for saltwater intrusion modeling.
Web site: [Link]
128
A.3 - SEAWAT–2000
A.5 - MOCDENS3D
reference density in the mass conservation equation. The computed velocities are then used in
the (non buoyant) solute transport model, which is the original transport model of MOC3D.
seawater intrusion on a regional scale. Water density are considered to vary either uniformly or
linearly within in each layer. That is, in the case of horizontal layers, the vertical density profiles
are either linear or else constant within each layer. A single input file is needed to obtained the
required modifications of MODFLOW to treat the flow problem with variable density flow
(specification of a distribution of sources).
Web site: [Link]
A.10 - SHARP
with the bottom of the aquifer) are located by linearly projecting a line defined by the interface
elevations in adjacent blocks until it intersects the top and bottom of the aquifer.
Web site: [Link]
A.12 - SALTFLOW
A.13 - HBGC
The BIGFLOW code, named “BF” for short, is the research code used by the french partner
(Partner 3 of SWIMED project) to develop analyses on porous media hydrodynamics and
hydrogeologic processes, particularly in the presence of heterogeneities and couplings.
Historically, BIGFLOW was initiated as a numerical tool for modeling 3D flow systems in
randomly heterogeneous geologic formations, considering first only saturated GroundWater flow
(Ababou et al. 1985) then unsaturated flow in heterogeneous soils (Ababou et al. 1987 & 1988)
porous media with high resolution. The code was extensively tested between 1988 and 1992
resulting in a published manual of BGFLOW version 1.1 was published by the U.S. Nuclear
Regulatory Commission (Ababou & Bagtzoglou 1993). The BIGFLOW code was later reshaped and
enhanced, leading to BF-2000, a computer code for modeling diverse flow processes in
heterogeneous variably saturated hydrologic media, including some surface/subsurface coupling.
The BF2000 code is well documented in the Ph.D. thesis of Trégarot (2000).
The spatially distributed flow model BF2000 generalizes variably saturated flow processes for
heterogeneous, anisotropic, macroporous, and partially saturated media in 3D. The model can also
efficiently track the dynamics of multiple interacting free surfaces in 3D. The model can
accommodate nonlinear, anisotropic and heterogeneous material properties, in (x,y,z) for the 3D
case, or in (x,y) for the 2D vertically integrated case. A new 2D module BFSWIM was introduced
134
more recently in BF2000 to take into account saltwater intrusion modeling, and also, internal
sink/source terms .
The BF code has two rather different options, 2D/3D, each one leading to various possible flow
regimes as illustrated in the flow-charts below option (BF-3D) and 2D option (BF2D):
• 3D option (BF-3D): finite volume model for fully heterogeneous, three-dimensional,
variably saturated and nonlinear (porous or macroporous) soil-aquifer systems (figure 1).
• 2D option (BF-2D): finite volume model for vertically integrated quasi-plane flows,
including Boussinesq-Dupuit aquifer flows, free surface hydraulics based on kinematic
wave, Darcy-Forchheimer flow in rough fractures, and the seawater intrusion module.
Figure 1: 3D options with BIGFLOW (note the 3D Saltwater intrusion model is not listed here).
Remarks on 2D vs. 3D
The 3D option is fully three-dimensional, in the sense that all material properties being arbitrarily
distributed in 3D space. On the other hand, the 2D option is vertically integrated and assumes quasi-
plane flow (vertically hydrostatic). It should be emphasized, however, that the internal
computational structure of the code is 3D. Indeed, the 2D option was developed by specializing the
general algorithms of the code without modifying its internal 3D structure.
The 3D algorithms of the BF code have been tested quite extensively since the late 1980’s, and
many features of the more recent 2D plane flow algorithms have been tested since the late 1990’s
(see references: Ababou et al.; Trégarot et al.).
Remarks on Sea Water Intrusion Modeling (SWIM)
In this report, we choose to emphasize the general capabilities of the 2D/3D code as they existed
at the start of the project, rather than the new seawater intrusion modeling capabilities. Indeed, for
the present project on coastal aquifers, a new algorithm for seawater intrusion was developed into
the 2D version of BF, called BF-SWIM2D, to be described in deliverable “D7: Seawater
B.2 - BIGFLOW - Model equations and physical basis (2D & 3D)
BIGFLOW is an integrated model for groundwater and hydrologic flow processes based on a
unified, generalized, conservative flux-divergence equation obtained from Darcy-Buckingham and
mass conservation (see equations in Figure 3 below).
Table 0.1: Variables and parameters of the generic flow equation BF (3D option)
Table 0.2: Variables and parameters of the generic flow equation BF (2D option)
137
• General inputs
These are to define the main aspect of the simulation like type of simulation, transient/steady state
simulation, and linear/nonlinear prob. There are also the less important identification code and date
of simulation.
• Domain geometry
BF assumes that the 3D domain has the shape of a rectangular box (parallelepiped rectangle).
However the box can be slanted at any angle with respect to the vertical. The finite difference mesh
is a fixed regular network of orthogonal links and nodes. Never the less physical boundaries like
substratum bathymetry and surface topography, used in the coupled surface-subsurface flow, can be
distributed uniformly or spatially over the domain.
• Boundary conditions and initial conditions
BF enables the use of three types of boundary conditions (BC) Dirichlet (fixed head), Neumann
(fixed flux) and zero head gradient BC. All types can be fixed, uniformly or spatially distributed.
Time dependence is also available. Mixed type BC can be defined separately for each BC. Mixed
type means that the type of BC can vary within the boundary plane. Initial head pressures can be
defined as fixed, uniformly distributed or spatially distributed.
• Forcing Terms
Internal forcing terms are given in two different ways either locally for source/sink terms (e.i.
pumping/recharge wells), or uniformly distributed over the domain. In the latter case these represent
recharge or infiltration rates. Forcing terms are also time-varying.
• Physical Properties
• For all cases: hydraulic properties are specified cell-by-cell (pixels in 2D, voxels in 3D)
• For subsurface flow: hydraulics properties for saturated media are the hydraulic conductivity (Ks)
and the specific storativity (Ss) both are spatially variable if needed. For spatially saturated or
unsaturated media the hydraulic conductivity (K) is a function of pressure head h. The soil
moisture retention curve can be chosen among several functions with multiple parameters, two of
which can be spatially distributed.
• For surface flow: where the diffusive-kinematic wave equations are applied, the hydraulic
properties are the friction coefficients (C) derived from the Chezy, Manning or Darcy-Weisbach
formulas.
• For macroporous media flow: another coefficient (γ) of the quadratic term in the generalized
Ward’s law is defined.
• Numerical solver parameters
These parameters allow the user to specify the choice of matrix solver and preconditioner, as well
as its convergence criteria, like the maximum number of iterations, the minimum error criteria, and
the type of norm used to compute those. These parameters also include the nonlinear Picard iteration
convergence criteria. In addition, machine dependent criteria can be specified, such as machine
precision, smallest and largest real floating point numbers, etc.
• Limitations
The geometry of the 3D domain is rectangular
(however, in 2D, the bedrock is specified as ZINF(x,y))
The finite volume grid is rectangular, uniform, cartesian
(similarly in 2D and in 3D)
138
Some useful features such as leaky layers in the 2D case are not yet implemented
The time step adaptivity and the nonlinear solution process are not always succesfull
(a problem common to most computer models for highly nonlinear and transient problems)
(a) (b)
Figure 4 Bigflow stencil (a) for Darcy model (b) for Darcy Forcheimer quadratic loss model
(notice the used green cells in that case)
139
A graphical user interface under Python® 6 and using the wxPython® library has been developed
for pre and post processing. Since the Python language is cross-platform, the interface is available
on virtual all operating systems. Part of the interface has been programmed during the WP3.1 of
SWIMED project.
The objective of the GUI (Graphical User Interface) is not to supply a commercial software
package for BIGFLOW, but to provide BIGFLOW users with a tool to manage the model inputs and
outputs for more efficient usage. This is mostly important for problems with complex geometry, but
also, for optimisation problems and multiple scenario analyses.
The Python interface manages directories and files with minimal user intervention. It also enables
users to input simulation parameters like BC’s, grid size, time step, convergence criteria, initial head
distribution, etc. Simulations are launched and tracked in real-time by numerical and graphical
outputs : see for example Figure(5 & 6) below. Outputs are also managed by the Python
application; they can be plotted and exported in different graphical formats. Multiple simulations can
be launched and managed interactively, for example Monte-Carlo simulations for stochastic
uncertainty analyses, or iterative re-starts as occured in the sigma-homotopy method 7.
The Python interface is intended to be a platform for hydrological modelling. In the latest version
the interface now has incorporated
a 3D conditional random media generator (XiMUL) (Ababou et al.)
a discrete random object generator (Bailly et al.)
a particule tracking code LPT3D (Ababou et al.)
The interface enables 3D visualization using the VTK library and the MayaVi interface.
Additionally GIS importing and exporting capabilities has been implemented.
represents the degree of heterogeneity (standard deviation of lnK); the method was developed during the WP
3.1 and was used to model seawater intrusion in highly heterogeneous stochastic aquifers.
140
Source term modelling: contours of the head field Unsaturated aquifer modelling: Vector plot
of fluxes and streamlines
In BIGFLOW (ref, RA 1993) all arrays are assembled in a total array. The total array is
dynamically allocated in the code. Its size changes depending on domain size and simulation
options. If an array is not used in a simulation its size will be reduced to minimal value (e.g. (1,1,1)
since all data arrays are three dimensional in BIGFLOW), the size of total array needed is then
reduced. This feature has been conserved in later versions of BIGFLOW namely BF2000 and
BFSWIM.
An estimate of the size of total array is set by the user before each compilation. This estimate
should be greater or equal to the effective size of the total array, but not much greater. Since this will
reduce performance and increase memory requirements.
The array size set before compilation defines the memory requirements of the code. The reserved
memory is than the product of the integer size and the total array dimension.
For a simple precision calculation the integer size is 4 bytes and for double precision its 8 bytes.
The used Random Access Memory RAM is than:
This value is limited by two features: the system RAM capacity and the reserved STACK. The
system RAM is an intrinsic feature of the PC or the station. As for the STACK it’s the maximum
reserved memory in the compilation options. The reserved STACK shouldn’t exceed system RAM
and not less than the required RAM value.
Since edge nodes are only considered for initial mid nodal fluxes and do not intervene in
computation, the effective number of nodes is therefore obtained by:
Nt = ∏i =1 ( Ni − 2)
3
Figure 1 illustrates the fitted function and the corresponding correlation coefficient for 2D sharp
interface salt water intrusion modeling (BFSWIMv1.1). The total array size is :
ABIG = 54.05 Nt + 0.0793
with ABIG in million elements ans Nt in million nodes
the memory requirement using single precision is :
RAM (Kb) = 211.4 Nt + 5.648
500 2000
Required Memory in Mb
R =1
Totla Array in millions
Required memory Mb
100 400
0 0
0 1 2 3 4 5 6 7 8 9 10
Nt (N1xN2x1) million nodes
144
∂ 2ϕ ∂ 2ϕ ∂f
+ − J = 0 avec J = −∂H / ∂x (on néglige les termes de second ordre)
∂x 2 ∂y 2 ∂x
+∞ +∞
ϕ = ∫e i ( kx )
dZ ϕ (k ) ϕ = ∫ ∫e
i ( kx )
dZ ϕ (k )
RN − ∞− ∞
+∞ +∞
f = ∫e dZ f (k ) f = ∫ ∫e
i ( kx ) i ( kx )
dZ f (k )
RN − ∞− ∞
The Spectral/Fourier transform implies the unicity of the transform and the orthogonality
The unicity of the Fourier representation gives:
− ik1 J
dZ ϕ (k ) = dZ f (k )
(k1 + k 2 ) 2
The Weiner-khinchine theorem gives the relation between the spectral density function and the
correlation by:
Rϕϕ (τ )dτ
1
∫e
− ikτ
Sϕϕ (k ) =
2π N
RN
Replacing the correlation function considered her is:
145
⎡ πτ ⎛ πτ ⎞ 1 ⎛ πτ ⎞ 2 ⎛ πτ ⎞⎤
R (τ ) = σ ⎢ K1 ⎜ ⎟ − ⎜ ⎟ K 0 ⎜ ⎟⎥
2
⎣⎢ 4λ ⎝ 4λ ⎠ 2 ⎝ 4λ ⎠ ⎝ 4λ ⎠⎦⎥
where Jx is the mean Φ-gradient denoted “a” in this paper. The values of a are given in the
following table:
sigma(ln(k)) L a=d(phi)/dx
0 581 1.5461
1.15 640 1.4036
1.32 690 1.3019
2.312519 701 1.2815
2.684751 704 1.2760
2.988887 705 1.2742
3.246030 706 1.2724
3.46877 706 1.2724
3.665255 706 1.2724
3.841010 706 1.2724
4.00 706 1.2724
In this paragraph we use an empirical upscaling technique for the saltwater/freshwater interface.
The upscaling is considered with reference to the total water flux
Qˆ =
− Kˆ p
2L
2
(
H1 − Z SALT (ε + 1) + 2 ε ΔZ Z SALT
2
)
The power average model (Ababou, 1994) for the empirical upscaling technique is then used. This
is given by the following formula:
Kˆ p = E ( K p )1/ p
with p=directional averaging exponents
2 λh
p = 1−
D λi
D space dimentions
λi=directional fluctuation scales
λh=D-dimentional harmonic mean fluctuation scale
146
−1
⎡1 D ⎤
λh = ⎢ ∑ λi −1 ⎥
⎣ D i =1 ⎦
For a Gaussian field the equivalent permeability is given by:
⎛1 ⎞
Kˆ p = E ( K p )1/ p = K G exp⎜ pσ y2 ⎟
⎝2 ⎠
Results are shown in figure below. The numerical results (x) are plotted against the analytical
results ( ) for different values of p (p=-1,0, 0.08, 1). The best fitted value of P is 0.01. This shows a
small deviation from the geometric mean.
Other values of interest for the stochastic analytical model are shown in this table.
147
Axial transects of fluctuations around the mean of One hundred transects of ΦSALT (transformed from
ΦSALT (see previous figure); simulation on ZSALT); “analytical mean” curve ΦSALT
1000x1000 grid for σ = 2.0. (homogeneous aquifer); and "numerical mean" curve
ΦSALT (mean of ΦSALT sampled shorewise along
“y”). The sea shore is at left. Grid: 1000x1000.
Heterogeneity: σ = 2.0.
std of the transformed of the transform ΦSALT (see Numerical and theoretical σZsalt vs distance from
previous figure); simulation on 1000x1000 grid for sea (x) for σlnK = 2.0.
σ = 2.0.
150
std of the transformed of the transform ΦSALT (see Numerical and theoretical σZSALT vs distance from sea (x) for σlnK =
ln(10).
previous figure); simulation on 1000x1000 grid
for σ = ln(10).
151