Joseph Lmpa

Liouville

Christophe Bourel
email: christophe.bourel@univ-littoral.fr

Catherine Choquet

Carole Rosier

Munkhgerel Tsegmid

Dérivation du modèle

Keywords: Fluid flow modelling, Saturated and unsaturated porous media, Numerical simulations, Asymptotic analysis, Vertical Richards equations, Dupuit Hypothesis

Tout d'abord, j'exprime ma profonde

Résumé

Nous présentons une classe de nouveaux modèles pour décrire les écoulements d'eau dans des aquifères peu profonds non confinés. Cette classe de modèles offre une alternative au modèle Richards 3d plus classique mais moins maniable.

Leur dérivation est guidée par deux ambitions : le nouveau modèle doit d'une part être peu coûteux en temps de calcul et doit d'autre part donner des résultats pertinents à toute échelle de temps. Deux types d'écoulements dominants apparaissent dans ce contexte lorsque le rapport de l'épaisseur sur la longueur de l'aquifère est petit : le premier écoulement apparaît en temps court et est décrit par un problème vertical Richards 1d ; le second correspond aux grandes échelles de temps, la charge hydraulique est alors considérée comme indépendante de la variable verticale. Ces deux types d'écoulements sont donc modélisés de manière appropriée par le couplage d'une équation 1d pour la partie insaturée de l'aquifère et d'une équation 2d pour la partie saturée. Ces équations sont couplées au niveau d'une interface de profondeur h(t , x) en dessous de laquelle l'hypothèse de Dupuit est vérifiée. Le couplage est assuré de telle sorte que la masse globale du système soit conservée. Notons que la profondeur h(t , x) peut être une inconnue du problème ou être fixée artificiellement. Nous prouvons( dans le cas d'aquifères minces) en utilisant des développements asymptotiques que le problème Richards 3d se comporte de la même manière que les modèles de cette classe à toutes les échelles de temps considérées (courte, moyenne et grande).

Nous décrivons un schéma numérique pour approcher le modèle couplé non linéaire. Une approximation par éléments finis est combinée à une méthode d'Euler implicite en temps. Ensuite, nous utilisons une reformulation de l'équation discrète en introduisant un opérateur de Dirichlet-to-Neumann pour gérer le couplage non linéaire en temps. Une méthode de point fixe est appliquée pour résoudre l'équation discrète reformulée. Le modèle couplé est testé numériquement dans différentes situations et pour différents types d'aquifère. Pour chacune des simulations, les résultats numériques obtenus sont en accord avec ceux obtenus à partir du problème de Richards original.

Nous concluons notre travail par l'analyse mathématique d'un modèle couplant le modèle Richards 3d à celui de Dupuit. Il diffère du premier parce que nous ne supposons plus un écoulement purement vertical dans la frange capillaire supérieure. Ce modèle consiste donc en un système couplé non linéaire d'équation Richards 3d avec une équation parabolique non linéaire décrivant l'évolution de l'interface h(t , x) entre les zones saturée et non saturée de l'aquifère. Les principales difficultés à résoudre sont celles inhérentes à l'équation 3D-Richards, la prise en compte de la frontière libre h(t , x) et la présence de termes dégénérés apparaissant dans les termes diffusifs et dans les dérivées temporelles.

Abstract

We present a class of new efficient models for water flow in shallow unconfined aquifers, giving an alternative to the classical but less tractable 3D-Richards model. Its derivation is guided by two ambitions : any new model should be low cost in computational time and should still give relevant results at every time scale. We thus keep track of two types of flow occurring in such a context and which are dominant when the ratio thickness over longitudinal length is small : the first one is dominant in a small time scale and is described by a vertical 1D-Richards problem ; the second one corresponds to a large time scale, when the evolution of the hydraulic head turns to become independent of the vertical variable. These two types of flow are appropriately modelled by, respectively, a one-dimensional and a two-dimensional system of PDEs boundary value problems. They are coupled along an artificial level below which the Dupuit hypothesis holds true (i.e. the vertical flow is instantaneous below the function h(t , x)) in a way ensuring that the global model is mass conservative. Tuning the artificial level, which even can depend on an unknown of the problem, we browse the new class of models. We prove using asymptotic expansions that the 3D-Richards problem and each model of the class behaves the same at every considered time scale (short, intermediate and large) in thin aquifers.

We describe a numerical scheme to approximate the non-linear coupled model. The standard Galerkin's finite element approximation in space and Backward Euler method in time are used for discretization. Then we reformulate the discrete equation by introducing the Dirichlet to Neumann operator to handle the nonlinear coupling in time. The fixed point iterative method is applied to solve the reformulated discrete equation. We have examined the coupled model in different boundary conditions and different aquifers. In the every situations, the numerical results of the coupled models fit well with the original Richards problem.

We conclude our work by the mathematical analysis of a model coupling 3D-Richards flow and Dupuit horizontal flow. It differs from the first one because we no longer assume a purely vertical flow in the upper capillary fringe. This model thus consists in a nonlinear coupled system of 3D-Richards equation with a nonlinear parabolic equation describing the evolution of the interface h(t , x) between the saturated and unsaturated zones of the aquifer. The main difficulties to be solved are those inherent to the 3D-Richards equation, the consideration of the free boundary h(t , x) and the presence of degenerate terms appearing in the diffusive terms and in the time derivatives.

. CONTEXTE GÉNÉRAL

L'eau est l'une des exigences essentielles de notre vie. Dans le monde, la consommation d'eau utilisée pour les activités humaines quotidiennes provient principalement de l'eau d'aquifères souterrains. Dans certains pays, la consommation totale d'eau provient des ressources souterraines. Par exemple, en Mongolie, 80% de la consommation d'eau provient des ressources souterraines (voir [26]), le problème étant que ces eaux ne sont souvent pas réutilisées. De plus l'urbanisation et l'industrialisation récentes sont à l'origine de la contamination des sols et de la détérioration des aquifères d'eau douce. Ainsi, l'étude des écoulements d'eau dans les milieux poreux (saturés et non saturés) est un enjeu important pour la consommation d'eau dans de nombreux domaines entre autre tels que l'agriculture, l'environnement, la gestion des déchets, le développement urbain et les processus industriels. Par ailleurs, les écoulements en milieux poreux sont de nature complexe. La nature de la formation des milieux poreux dans les aquifères souterrains sont souvent difficiles à décrire du fait de leur spécificité (comme par exemple la géométrie de l'aquifère, l'hétérogénéité du sol, l'anisotropie des sédiments ou du substrat rocheux dans l'aquifère, les mécanismes de transport des contaminants et des réactions chimiques). La surveillance de tels écoulements est donc prohibitive. Ceci oblige les scientifiques et les ingénieurs à s'appuyer fortement sur des modèles mathématiques pour comprendre et prévoir le comportement du débit d'eau souterraine.

En 1856, Henry Darcy présenta une formule décrivant l'écoulement à travers un matériau poreux. Il déduisit cette formule d'une manière phénoménologique en consignant ses expériences dans des notes pour un système d'alimentation en eau de la ville française de Dijon. Lorenzo Adolph Richards a ensuite publié un modèle décrivant un écoulement d'eau dans un domaine non saturé. Ce modèle est basé sur une extension de la loi Darcy à un écoulement polyphasique réalisé par Edgar Buckingham. La modélisation des écoulements souterrains a connu un grand essor au cours de ces sept ou huit dernières décennies. Ces modèles mathématiques sont basés sur 9 des hypothèses simplificatrices des milieux poreux. Ils sont en effet une approximation et non une duplication exacte des conditions sur le terrain. Toutefois, même approximatifs, ils constituent un outil d'investigation utile que les hydrologues peuvent utiliser pour diverses applications, telles que la simulation du débit d'eau, la migration chimique dans les zones saturées, les échanges entre les rivières et les eaux souterraines. Ces simulations peuvent aussi servir de base pour la création de zones de protection des eaux souterraines.

Le présent travail traite des phénomènes naturels liées à ces questions hydrogéologiques. En général, le mouvement des eaux souterraines est considéré comme un problème de fluide polyphasique qui est décrit par les équations de Richards [38]. Il s'agit d'un système tridimensionnel d'équations non linéaires dégénérées de type parabolique. Notre travail repose principalement sur ces équations. Nous nous concentrons sur un écoulement fluide monophasique incompressible (l'eau) à densité constante dans des aquifères très minces et grands. Nous observons alors qu'il existe deux types d'écoulements dominants dans ce type d'aquifère. Ils se produisent dans différentes régions de l'aquifère et à des échelles de temps différentes. Le premier flux dominant apparaît dans la région supérieure et il est globalement vertical. Ce flux est très rapide comparé à l'autre et il se produit sur des temps courts. Nous l'avons donc appelé composante rapide de l'écoulement. Le second apparaît dans la région inférieure de l'aquifère et il se fait globalement dans la direction horizontale. Il a une vitesse très lente comparée à la vitesse verticale dans la frange capillaire et il se produit sur une période de temps très longue. Pour cela, nous l'avons appelé composante lente de l'écoulement. Ces deux composantes de l'écoulement sont séparées par une interface. Cette interface est donc l'intersection de ces deux régions. La région inférieure est complètement saturée tandis que la partie supérieure est partiellement saturée (pouvant même, éventuellement être complètement sèche). En particulier au-dessus de l'interface, le flux vertical est dominant alors qu'en dessous, il est quasiment instantané. Ainsi, l'hypothèse de Dupuit est satisfaite dans la région inférieure du réservoir. Cela permet l'intégration verticale des équations de Richards dans cette partie et conduit à l'utilisation d'une famille de modèles 2D fortement développés depuis les années 60 (voir, par exemple, les travaux de Jacob Bear, [7,8]).

Néanmoins, l'approximation de Dupuit est surtout utile pour des phénomènes à long terme, la composante rapide de l'écoulement n'étant pas prise en compte dans ces modèles. Il est pourtant important de bien décrire l'écoulement dans la zone insaturée par exemple à cause des réactions chimiques qui s'y produisent principalement. Par ailleurs, la partie supérieure constitue souvent la recharge principale en eau de la zone basse de l'aquifère. Cela tend à encourager les modèles couplant les équations de Richards à l'approximation de Dupuit. De nombreux travaux ont été proposés dans ce sens. Par exemple, dans [37], les auteurs considèrent l'équation Richards 1D couplée à un modèle simplifié dans la partie saturée. Cette étude est purement numérique et le modèle n'est pas mathématiquement justifié. Dans [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF], ce type de modèle est intégré à un code de calcul appelé "SHE" (pour "Système hydrologique européen", devenu plus tard SHETRAN) dans le cas où la nappe reste éloignée du sol. Signalons aussi les travaux de [START_REF] Yakirevich | A quasi three-dimensional model for flow and transport in unsaturated and saturated zones : 1. implementation of the quasi two-dimensional case[END_REF] et [35] pour des modèles similaires et le modèle bidimensionnel directement couplé à un modèle de surface proposé dans [29]. Nous soulignons enfin que la vraie difficulté réside dans la manière dont on couple les deux modèles. Dans le second chapitre de ce travail, nous proposons un nouveau modèle pour lequel le couplage se biais par le frais d'un terme assurant la conservation de la masse totale du système. Nous justifions ensuite ce modèle en montrant formellement qu'il a les mêmes comportements asymptotiques à différentes échelles de temps que ceux obtenus pour Richards 3d .

Le deuxième objectif est d'illustrer numériquement l'efficacité du modèle couplé. C'est l'objet du troisième chapitre. Nous préconisons une discrétisation éléments finis combinée avec un schéma temporel de type Euler implicite (ce qui n'est pas sans difficulté en raison de la non linéarité du problème). Les simulations numériques illustrent parfaitement les performances du nouveau modèle testé dans plusieurs situations et donnant des solutions très proches de celles obtenues avec les équations Richards 3d dans tout le domaine d'étude.

Enfin nous concluons ce travail par une première étape dans l'étude théorique de ce modèle (ce qui constitue le dernier chapitre). A cette fin, nous présentons un modèle appartenant à la classe des modèles "Dupuit-Richards". Ce modèle diffère du précédent, déjà parce que nous considérons les équations complètes de Richards 3d dans la partie non saturée. Mais la principale différence réside dans le couplage entre les deux zones. En effet, pour le premier modèle, le couplage se fait par des termes de flux garantissant la conservation de la masse, alors que dans le second, nous imposons des propriétés de transmission pour la pression et pour les flux normaux à l'interface de saturation.

Les résultats de ce travail sont donnés sous forme de trois articles (correspondants aux trois derniers chapitres) qui sont rédigés en anglais. L'introduction rédigée en français donne un aperçu détaillé de ces résultats.

ÉQUATION DE RICHARDS -APPROXIMATION DE DUPUIT

Les aquifères sont souvent caractérisés par une forme de stratification des écoulements qui permet la définition d'interfaces, la lenteur de la dynamique naturelle assure que ces interfaces sont régulières et ont un comportement stable. Ce point ajouté au fait que les écoulements sont essentiellement orthogonaux aux équipotentielles (hypothèse de Dupuit) permet l'intégration verticale de l'équation de Richards dans la zone saturée. Dans cet esprit, de nombreux modèles 2D ont été développés et sont souvent utilisés depuis les années 60 (voir par exemple les travaux de Jacob Bear, [7,8]). Nous allons rappeler dans cette section l'obtention des équations de Richards ainsi que le principe de l'approximation de Dupuit.

Lois de conservations

Dans cette section, nous allons introduire les équations fondamentales qui sont communément utilisées en hydrogéologie ainsi que les paramètres physiques impliqués dans ces équations. Les équation de Richards constitue un modèle classique pour décrire le mouvement d'un fluide en milieu poreux soumis à l'action de la capillarité et de la gravité. Il est attribué à Lorenzo A. Richards qui a publié l'équation en 1931. Mais, en fait, cette équation a été publiée pour la première fois par le mathématicien et physicien anglais Lewis Fry Richardson dans son livre "Prévision météorologique par processus numérique" publié en 1922 (cf. [39]). Loi de Darcy. Il s'agit d'une relation mathématique découverte en 1856 par l'ingénieur français Henri Darcy. Cette loi régit les écoulements d'un fluide à travers des matériaux poreux et elle permet de calculer la quantité d'eau s'écoulant dans ces milieux. Cette relation découle de la conservation du moment cinétique.

Compte tenu des dimensions importantes d'un aquifère par rapport à la taille caractéristique de la structure poreuse du sous-sol, nous considérons une description continue du milieu poreux. La vitesse effective v de l'écoulement est donc liée à la pression P par la loi de Darcy donnée ci-dessous

v = - κ K 0 µ (∇P + ρg ∇z), (1.2.1)
où z désigne la hauteur, ρ et µ sont respectivement la masse volumique et la viscosité du fluide, K 0 est la perméabilité du sol, κ est la conductivité relative et g est la constante d'accélération gravitationnelle.

Dans l'expression (1.2.1), le terme en ∇P représente l'énergie due à la pression du fluide sur les pores alors que le terme en ∇z correspond à l'énergie potentielle gravitationnelle. L'énergie cinétique est négligée en raison du mouvement souvent très lent des eaux souterraines.

La conductivité hydraulique caractérise la facilité avec laquelle un fluide peut s'écouler à travers un matériau. Cela dépend de nombreux facteurs (tels que la viscosité et la masse volumique du fluide mais aussi des caractériques du sous-sol). Elle reflète donc principalement la perméabilité du sol.

On introduit la charge hydraulique qui correspond à l'énergie mécanique disponible pour entraîner le débit d'un volume donné d'eau. Elle est définie comme suit :

H (t , x, z) = P (t , x, z) ρg + z. (1.2.2)
La loi de Darcy peut alors s'écrire v = -K ∇H -κ K 0 µ (ρρ 0 )g ∇z, (1.2.3) ρ 0 désignant la masse volumique de référence du fluide. Dans le cas où le fluide est supposé incompressible, l'équation (1.2.3) se simplifie en

v = -K ∇H , avec K = k r K o (k r = κ ρ 0 g µ ), (1.2.4)
où l'on a noté K la conductivité hydraulique.

Conservation de la masse. La conservation de la masse pendant le déplacement est donnée par l'équation suivante

∂(φρs) ∂t + div(ρv) = Q s (1.2.5)
où Q s désigne un terme source générique (correspondant à la production et/ou au réapprovisionnement). La fonction θ est la teneur volumétrique en humidité définie par θ = φs, où φ est la porosité du milieu et s la saturation. Si on suppose que l'air présent dans la zone non saturée a une mobilité infinie, la saturationt s et la fonction θ sont donc considérées comme des fonctions monotones dépendant de la pression comme nous le détaillerons plus loin.

Tenseur de perméabilité K 0 . La perméabilité K 0 (x, z) est un tenseur (3 × 3) symétrique défini positif décrivant la conductivité du sol saturé au point (x, z) ∈ Ω. On introduit de plus les notations K xx ∈ M 22 (R), K zz ∈ R * et K xz ∈ M 21 (R) de telle sorte que ρ = ρ 0 e α P (P -P 0 ) .

K 0 = K xx K xz K T xz K zz . ( 1
(1.2.7)

Le réel α P ≥ 0 est le coefficient de compressibilité du fluide et P 0 est la pression de référence. En supposant que α P = 0 nous récupérons le cas incompressible.

Équation de Richards

Enumérons maintenant les hypothèses sur les caractéristiques des fluides et des milieux, mais aussi sur celles des écoulements qui sont significatives dans le contexte de notre problème.

Hypothèses sur le fluide et sur le milieu.

-Compressibilité du sol. Nous négligeons dans le modèle les effets de la compressibilité de la roche, la porosité du milieu φ ne dépend pas des variations de pression et elle est donc supposée être constante. -Compressibilité du fluide. Nous supposons que le fluide (à savoir ici l'eau douce) est faiblement compressible. Cela signifie que α P 1.

(1.2.8)

Obtention de l'équation de Richards. Exploitons cette hypothèse. Dans des conditions naturelles et en particulier dans un aquifère, on observe une faible mobilité des fluides (définie par le rapport κ/µ). La première conséquence de la faible compressibilité du fluide combinée à la faible mobilité du fluide va se traduire dans l'équation du moment cinétique. Nous effectuons un développement de Taylor par rapport à P de la masse volumique ρ dans le terme de gravité de l'équation de Darcy. En négligeant les termes pondérés par α P κ/µ 1 dans (1.2.3), on obtient :

v = -K ∇H , K = κ(P ) ρ 0 g µ K 0 . (1.2.9)
La seconde conséquence est que ∇ρ•v 1 ce qui conduit à la simplification suivante dans l'équation de la conservation de la masse (1.2.5) :

ρ∂ t θ + θ ∂ t ρ + ρ div v = ρQ s .
Négliger ainsi la variation de densité dans la direction de l'écoulement est parfois considéré comme une hypothèse supplémentaire appelée hypothèse de Bear (cf [2]). En incluant (1.2.7), i.e. ∂ t ρ = ρα P ∂ t P dans la précédente équation, on obtient

ρ∂ t θ + ρ θα P ∂ t P + ρ div v = ρQ s .
Après simplification puisque ρ > 0, on obtient finalement

∂ t θ + θα P ∂ t P + div v = Q s . (1.2.10)
De manière équivalente, en utilisant la charge hydraulique (1.2.2) et la loi de Darcy (1.2.9), (1.2.10) devient

∂ t θ + S 0 ∂ t H -div(K ∇H ) = Q s où S 0 = ρ 0 g φ α P . (1.2.11)
On remarque que si le fluide est supposé incompressible, α P = 0, alors l'équation (1.2.10) est l'équation de Richards en formulation pression sous sa forme classique. Une définition adéquate de la teneur volumétrique en humidité θ ainsi que de celle de la mobilité κ constituent la clef du modèle.

Hypothèse de Richards. Le modèle de Richards repose en outre sur l'hypothèse que la pression atmosphérique dans le sous-sol est égale à la pression atmosphérique, ce n'est donc pas une inconnue du problème. On suppose que la saturation et la conductivité relative de sol sont des fonctions dépendant de la pression du fluide P , notées respectivement s = s(P ) et κ = κ(P ). Nous introduisons la pression de saturation P s qui est un nombre réel fixe. La partie complètement saturée du support correspond à la région {(x, z) ∈ Ω , P (•, x, z) > P s }, alors qu'il est partiellement saturé dans la frange capillaire {(x, z) ∈ Ω , P d < P (•, x, z) ≤ P s }.

La partie sèche est définie par l'ensemble {(x, z) ∈ Ω , P (•, x, z) ≤ P d }.

La teneur en humidité est telle que θ =    φ (zone saturée) si P (•, x, z) > P s , θ(P ) (avec 0 ≤ θ(P ) ≤ φ et θ (P ) > 0) si P d < P (•, x, z) ≤ P s , 0 (zone sèche) si P (•, x, z) ≤ P d .

(1.2.12)

La mobilité hydraulique relative associée est alors définie par κ(θ) =    1 (zone saturée) si P (•, x, z) > P s , κ(θ(P )) (avec 0 ≤ κ(θ(P )) ≤ 1 et (κoθ) (P ) > 0) si P d < P (•, x, z) ≤ P s , 0 (zone sèche) si P (•, x, z) ≤ P d .

(1.2.13) Il existe un grand choix de modèles disponibles pour s et κ. Les exemples les plus classiques sont les modèles de Van Genuchten [24], sans dépendance explicite de la pression de saturation mais avec des paramètres d'ajustement, et le modèle de Brooks et Corey [13] En particulier, cela signifie que la pression de l'eau est supérieure à la pression de saturation P s si et seulement si le sol est complètement saturé.

Les graphes des fonctions s et κ fournies par le modèle de Brooks-Corey (qu'on a utilisé pour les simulations numériques ) sont représentées à la Figure 1.1.

Géométrie de l'aquifère

On considère des aquifères de forme cylindrique correspondant au domaine Ω ⊂ R 3 . La projection de Ω sur le plan horizontal est un domaine ouvert Ω x ⊂ R 2 de frontière ∂Ω x . Les bases inférieures et supérieures de Ω sont définies respectivement par les graphes des fonctions h bot (x) et h soil (x)

h soil (x) > h bot (x) , ∀x ∈ Ω x . (1.2.16)
Le domaine est alors donné par :

Ω = (x, z) ∈ Ω x × R | z ∈ h bot (x), h soil (x) .
(1.2.17)

On décompose la frontière ∂Ω de Ω en trois parties (bas, haut et latéral) comme suit

∂Ω = Γ bot Γ soil Γ ver , Γ bot := (x, z) ∈ Ω | z = h bot (x) , Γ soil := (x, z) ∈ Ω | z = h soil (x) , Γ ver := (x, z) ∈ Ω | x ∈ ∂Ω x .
Dans le présent travail, nous dérivons une classe de modèles caractérisés par la position h d'une interface virtuelle dans le réservoir. Pour notre construction, cette fonction doit prendre ses valeurs dans l'intervalle semi-ouvert [h bot , h soil ). Commençons par présenter deux sous-régions auxiliaires de Ω dans lesquelles l'écoulement présentera un comportement très différent. La définition de ces sous-régions est basée sur celle de la fonction h = h(t , x) qui pourra être l'une des inconnues de notre modèle. La sous-région supérieure à h = h(t , x) est notée Ω + h (t ) et la sous-région inférieure est notée Ω - h (t ). Elles sont définies comme suit :

Ω + h (t ) := (x, z) ∈ Ω | z > h(x, t ) and Ω - h (t ) := (x, z) ∈ Ω | z < h(x, t ) , Γ h := (x, z) ∈ Ω | z = h(x, t ) .

Approximation de Dupuit dans la zone saturée

Hypothèse sur le fluide L'hypothèse suivante est introduite pour moyenner le problème 3d en un problèm 2d dans la partie saturée inférieure du domaine.

Approximation de Dupuit (approche hydrostatique)

L'hypothèse de Dupuit consiste à considérer la charge hydraulique constante dans la direction verticale (équipotentielles verticales). C'est légitime puisque l'on observe en réalité des déplacements quasi-horizontaux lorsque l'épaisseur de l'aquifère est faible par rapport à ses dimensions horizontales et lorsque l'écoulement est loin des puits et des sources.

Procédure de mise à l'échelle

S 0 ∂ t H + div v d z = h h bot Q s d z.
On note par B f = hh bot l'épaisseur de la zone saturée et par Q le terme source représentant l'alimentation en eau douce distribuée en surface de l'aquifère :

Q = 1 B f h h bot Q s d z.
En appliquant la règle de Leibnitz au premier terme du membre gauche de l'égalité, on obtient :

h h bot S 0 ∂ t H d z = S 0 ∂ ∂t h h bot H d z -S 0 H |z=h ∂ t h + S 0 H |z=h bot ∂ t h bot .
On introduit H la moyenne verticale de la charge hydraulique

H = 1 B f h h bot H d z. L'approximation de Dupuit entraine que H (x 1 , x 2 , z) H (x 1 , x 2 ), x = (x 1 , x 2 ) ∈ Ω, z ∈ (h bot , h). Nous avons alors h h bot S 0 ∂ t H d z = S 0 B f ∂ t H .

De la même façon, nous avons

h h bot div v d z = div x (B f ṽ ) + v |z=h -• ∇(z -h) -v |z=h + bot • ∇(z -h bot ), où ∇ = (∂ x 1 , ∂ x 2 ), v = (v x 1 , v x 2 ) et la vitesse moyenne de Darcy ṽ = 1 B f h h bot v d z est donnée par ṽ = - 1 B f h h bot K ∇ x H d z = - 1 B f h h bot K ∇ x H d z = -K ∇ x H , K = 1 B f h h bot K 0 ρ 0 g µ d z,
(rappelons que κ(P ) = 1 pour z ∈ (h bot , h)). Finalement, la loi de conservation de la masse moyennée de l'eau douce dans la zone saturée s'écrit donc -Expression du flux à travers l'interface : L'interface est définie par l'équation cartésienne

S 0 B f ∂ t H = div x (B f K ∇ x H ) + v |z=h + bot • ∇(z -h bot ) -v |z=h -• ∇(z -h) + B f Q. ( 1
F (x 1 , x 2 , z, t ) = 0 c'est à dire z -h(x 1 , x 2 , t ) = 0, le vecteur normal unitaire à l'interface n est donc colinéaire à ∇(z -h).
Si aucun transfert de masse ne se produit entre les deux zones, la composante normale de la vitesse effective est continue à l'interface z = h. La relation régissant la continuité de la composante normale de la vitesse s'écrit donc

v |z=h φ -q • n = 0,
où q désigne la vitesse de l'interface. On précise que l'interface est mobile donc au lieu d'avoir un flux de matière à travers l'interface, on considère que l'interface peut monter ou descendre. Elle satisfait

-∂ t h + q • ∇(z -h) = 0. Ainsi v |z=h + -v |z=h -• n = 0 ⇔ v |z=h + • ∇(z -h) = v |z=h -• ∇(z -h) = φ∂ t h
P (t , x, z) = ρ 0 g H (t , x) -z pour t ∈ [0, T [ , (x, z) ∈ Ω - h (t ). (1.2.21)
Par ailleurs, la pression est continue à travers l'interface Γ h , il s'en suit que

P (t , x, h -) = P (t , x, h + ) = P s ⇐⇒ H = P s ρ 0 g + h. ( 1
S 0 B f ∂ t h -div x (B f K ∇ x h) = B f Q -v |z=h + • ∇(z -h) in Ω x , (1.2.23) avec B f = (h -h bot ), K = 1 B f h h bot K 0 ρ 0 g µ d z et S 0 = ρ 0 g φ α P . (1.2.24)

Remarques

On note que dans le cas où le fluide sera supposé incompressible, le problème de Richards 3d dans la géométrie décrite dans la partie 1.2.3 prendra la forme suivante. Trouver la pression P et la vitesse v tels que :

             φ ∂s(P ) ∂t + div(v) = 0 dans ]0, T [×Ω v = -k r (P ) K 0 (∇P + e 3 ) dans ]0, T [×Ω α P + β v • n = F sur ]0, T [×Γ soil v • n = 0 sur ]0, T [×(Γ bot ∪ Γ ver ) (1.2.25)
où e 3 désigne le vecteur unitaire vertical orienté vers le haut. Cette situation sera considérée dans les Chapitres 2 et 3. Dans le Chapitre 4, cette compressibilité du fluide sera conservée. D'autre part, dans les Chapitres 2 et 3, on considèrera que la perméabilité K 0 (x, z) est le tenseur (3 × 3) défini par (1.2.6). En revanche, dans la section 1.5, le tenseur K 0 sera réduit à un scalaire pour ne pas alourdir les démonstrations.

UN MODÈLE COUPLÉ POUR DÉCRIRE L'ÉCOULEMENT DANS DES

AQUIFÈRES PEU PROFONDS

Dans cette section, nous allons présenter brièvement le résultat principal de cette thèse. Celui-ci consiste en l'obtention d'une classe de modèles qui décrivent l'écoulement de l'eau dans un aquifère très mince et très large. Chaque modèle de cette classe approche l'équation de Richards 3d introduite dans (1.2.1). L'objectif principal étant de donner un modèle proche de l'original tout en étant plus simple à traiter numériquement. Dans chaque modèle de cette classe, nous considérons un aquifère occupant un domaine géométrique peu profond par rapport à ses dimensions horizontales. Ce type d'hypothèses étant essentiellement assez peu restrictives puisque souvent vérifiées dans la nature.

Il s'avère que deux types d'écoulement dominant se superposent dans ce type d'aquifères peu profonds. Le premier est une composante rapide de l'écoulement et a lieu principalement dans une direction verticale. Le deuxième est une composante plus lente de l'écoulement et a lieu globalement en direction horizontale. Il s'avère aussi que le profil de pression correspondant à ce dernier type d'écoulement est proche du profil hydrostatique. En particulier, l'écoulement vertical apparaît comme étant instantané.

La classe de modèles que nous proposons est basée sur le couplage de ces deux types d'écoulement. Les modèles couplés obtenus conservent la masse et sont justifiés à la fin du Chapitre 2 à l'aide d'arguments d'analyse asymptotique.

Description de la classe de modèles couplés

On se place dans la suite de cette partie dans le cas d'un aquifère ayant la géométrie décrite dans la partie 1.2.3. On reprend également le cadre donné dans la partie 1.2.1.

Pour décrire le problème couplé, on commence par considérer deux sous-régions de l'aquifère Ω dans lesquelles l'écoulement présentera des comportements très différents. Ces sous-régions sont basées sur une fonction h = h(t , x) pouvant éventuellement être une inconnue du modèle. La partie au dessus du graphe de h est notée Ω + h (t ) et celle en dessous est notée Ω - h (t ). Elles sont données par

Ω + h (t ) := (x, z) ∈ Ω | z > h(x, t ) et Ω - h (t ) := (x, z) ∈ Ω | z < h(x, t ) , (1.3.1)
Γ h := (x, z) ∈ Ω | z = h(x, t ) . (1.3.2)
Comme nous le verrons dans la suite, il est crucial que la fonction h soit telle que

h bot (x) ≤ h(t , x) ≤ max min H (t , x), h max (x) , h bot (x) . (1.3.3)
La hauteur h max (x) est définè par

h max = h soil -δ, 0 < δ 1. (1.3.4)
En pratique, nous allons utiliser une caractérisation explicite de h en fonction de certaines autres inconnues du problème. Pour cela on introduit une fonction positive R (éventuellement dépendante de H ) ainsi que la fonction Q = Q(x, H ) suivante :

Q(x, H ) = max min H (t , x) -R, h max (x) , h bot (x) . (1.3.5)
Bien que de nombreux choix soient possibles pour R, nous allons nous concentrer dans la suite sur le seul cas d'une fonction constante R ≥ 0. On peut également remarquer que l'on a Q(x, H (t , x)) = h bot (x) pour tout (t , x) ∈]0, T [×Ω x si R est choisi suffisamment grand. On introduit la matrice M 0 de taille 3 × 3 suivante :

M 0 = S 0 0 0 0 , S 0 = K xx - 1 K zz K xz K zx .
Celle-ci jouera le rôle d'un tenseur de perméabilité effective. On introduit également le tenseur de conductivité moyen K ( H ) défini dans ]0,

T [×Ω x par K ( H )(t , x) = h soil (x) h bot (x) k r ρ g ( H (t , x) -z) M 0 (x, z) d z. (1.3.6)
Pour simplifier les notations, on utilisera la convention suivante

∇ x = (∂ x 1 , ∂ x 2 , 0) T pour la partie horizontale du gradient et div x (v) = ∇ x • v = ∂ x 1 v 1 +∂ x 2 v 2 pour la diver- gence horizontale de v ∈ R 3 .
La classe de modèles. Les modèles couplés consistent à caractériser l'écoulement au travers de sa pression P , sa vitesse v ainsi que des inconnues auxiliaires u, w, H et h telles que :

-La vitesse v est définie dans

Ω par          v = u + w pour t ∈]0, T [ , (x, z) ∈ Ω u = -k r (P ) ∂P ∂z + 1 K 0 e 3 pour t ∈]0, T [ , (x, z) ∈ Ω w = -k r H -z M 0 ∇ x H pour t ∈]0, T [ , (x, z) ∈ Ω (1.3.7) -Dans Ω + h (t ), l'équation de Richards verticale est satisfaite              φ ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 pour t ∈]0, T [ , (x, z) ∈ Ω + h (t ) α P + β u • e 3 = F pour (t , x) ∈]0, T [×Γ soil P t , x, h(t , x) = H (t , x) -h(t , x) pour (t , x) ∈]0, T [×Ω x P (0, x, z) = P init (x, z) pour (x, z) ∈ Ω + h (0) (1.3.8) -Dans Ω - h (t ) la pression P admet un profil hydrostatique P (t , x, z) = H (t , x) -z pour t ∈ [0, T [ , (x, z) ∈ Ω - h (t ) (1.3.9) -La charge hydraulique vérifie dans Ω x        div x K ( H ) ∇ x H = (u • e 3 ) Γ + h pour (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 pour (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) pour x ∈ Ω x (1.3.10) où (u • e 3 ) Γ + h
est la trace de u • e 3 sur Γ h depuis le dessus.

-Le niveau z = h en dessous duquel l'écoulement vertical est supposé instantanée est défini par

h(t , x) = Q(x, H (t , x)). (1.3.11)
En particulier, le problème couplé (1.3.7)-(1.3.11) dépendant de notre choix de fonction h(t , x) (c'est à dire de Q). Bien que tous les choix intermédiaires soient possibles, nous nous concentrerons dans la suite sur les cas extrémaux

Q(t , x) = h bot (x),
(1.3.12)

Q(t , x) = max min H (t , x), h max (x) , h bot (x) := h s (t , x), (1.3.13)
ainsi que sur le choix intermédiaire

Q(t , x) = max min H (t , x) -R, h max (x) , h bot (x) . (1.3.14)
La dépendance de la solution vis-a-vis de ces choix est discutée dans la partie numérique (Chapitre 3). La première remarque concernant les modèles précédents est qu'ils ont deux avantages par rapport au modèle de Richards 3d dont ils sont issus. Premièrement ils engendreront des problèmes numériques plus rapides à résoudre. En effet le problème initial 3d est remplacé par le couplage d'un problème 2d avec une multitude de problèmes 1d indépendants (et qui pourront être résolus en parallèle).

Deuxièmement, le problème couplé et le problème original de Richards 3d présentent les mêmes comportements dominants lorsque le ratio ε=profondeur/largeur de l'aquifère est petit. En effet, on prouve dans la partie 1.3.4 que ces deux problèmes admettent exactement les mêmes comportements asymptotiques lorsque ε → 0. De plus ces comportements effectifs sont identiques pour tous les choix d'échelles de temps considérés et indépendamment du choix de R ≤ 0 dans (1.3.5).

Commentaires sur le modèle couplé dans le cas (1.3.12)

Dans cette partie, on suppose que R est suffisamment grand (et/ou qu'il y ait suffisamment peu d'eau dans la nappe phréatique) pour que h = h bot . Il s'en suit donc que

Ω + h = Ω, Ω - h = et que Γ h = Γ bot (
         v = u + w pour t ∈]0, T [ , (x, z) ∈ Ω u = -k r (P ) ∂P ∂z + 1 K 0 e 3 pour t ∈]0, T [ , (x, z) ∈ Ω w = -k r H -z M 0 ∇ x H pour t ∈]0, T [ , (x, z) ∈ Ω (1.3.15)              φ ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 pour t ∈]0, T [ , (x, z) ∈ Ω α P + β u • e 3 = F pour (t , x, z) ∈]0, T [×Γ soil P = H -h bot pour (t , x, z) ∈]0, T [×Γ bot P (0, x, z) = P init (x, z) pour (x, z) ∈ Ω (1.3.16)      div x K ( H ) ∇ x H = (u • e 3 )| Γ bot pour (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 pour (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) pour x ∈ Ω x (1.3.
h(t , x) = Q(x, H (t , x)) = max min H (t , x) -R, h max (x) , h bot (x) . (1.3.19)
Comme précédemment, la vitesse de l'écoulement v résulte de la contribution de la composante rapide u et de la composante lente w. On note de plus que dans ce cas, l'ensemble Ω - h n'est plus vide en général. Ainsi, il est nécessaire de décrire l'écoulement dans les deux parties

Ω + h (t ) et Ω - h (t
P t , x, h(t , x)      = R si h bot (x) < h(t , x) < h max (x), ≥ R si h(t , x) = h max (x), ≤ R si h(t , x) = h bot (x).
(1. Composante lente de l'écoulement. Comme dans le cas précédent, nous introduisons la même pression auxiliaire P a (t , x, z) = H (t , x) -z. À présent on a P = P a dans tout Ω - h (t ). De plus comme précédemment, la continuité de P sur Γ h fait que l'on a asymptotiquement P = P a dans tout Ω lorsque l'aquifère est très fin et que la durée de l'expérience est suffisamment longue.

L'évolution de (P a , w) est caractérisée par l'évolution de H donné dans (1.3.10). Il s'agit du même type d'équation que dans le cas précédent. La difficulté ici étant que le second membre de cette équation fait intervenir le flux sur l'interface Γ h venant de la partie supérieure, interface qui dépend elle-même de la solution. En particulier une description précise de sa dépendance vis-à-vis de l'inconnue H sera nécessaire pour l'approximation numérique de ce problème (voir Chapitre 3).

Développements asymptotiques formels

L'objectif de cette section est d'obtenir et de comparer les problèmes effectifs issus du problème original de Richards 3d (1.2.25) et issus du problème couplé (1.3.7)- (1.3.11). Ces problèmes effectifs donnent une information précise sur les comportements dominants des solutions lorsque l'aquifère est très peu profond et très large et pour différentes échelles de temps.

On commence par introduire les problèmes effectifs issus du problème de Richards 3d donné dans (1.2.25). L'idée pour les obtenir est d'utiliser des arguments d'analyse asymptotique lorsque le ratio ε=profondeur/largeur de l'aquifère tend vers 0. On commence par introduire une version adimensionnée des équations.

On introduit les nombres positifs T , L x , L z qui représentent respectivement les temps, largeur et profondeur caractéristiques. Alors on introduit les variables sans dimension suivantes :

x = x L x , z = z L z , t = t T .
Puisque l'aquifère est supposé très large et peu profond, on va supposer que pour chaque réel 0 < ε 1, on a L z = 1 et L x = 1/ε de sorte que L z /L x = ε. Le problème de Richards 3d adimensionné est donné par

1 T φ ∂s(P ) ∂t -ε 2 div x k r (P ) K 0 ∇ x P - ∂ ∂z k r (P ) K 0 ∂P ∂z + 1 = 0 dans ]0, 1[×Ω.
(1.3.22) Il y a plusieurs choix convenables pour le temps caractéristique T . Le premier est T = 1 et va permettre de caractériser les phénomènes rapides. Le second est T = 1/ε 2 et va décrire les phénomènes lents. D'autres choix sont possibles mais conduisent à des modèles effectifs plus ou moins triviaux (noter que dans le chapitre 2 le cas T = 1/ε est aussi traité).

Pour γ ∈ {0, 2} on note P γ ε la solution du problème (1.3.22) pour T = ε -γ (on renvoie au chapitre 2 pour les conditions au bord rééchelonnées). On suppose que ces solutions vérifient les développements asymptotiques formels suivants :

P γ ε = P γ 0 + ε P γ 1 + ε 2 P γ 2 + . . . (1.3.23)
L'idée à ce stade est de caractériser formellement le terme dominant P γ 0 de ces développements. On les obtient comme solutions des problèmes effectifs suivants : Pour la composante rapide de l'écoulement

         φ ∂s(P 0 0 ) ∂t + ∂v 0 0 • e 3 ∂z = 0 dans ]0, 1[×Ω v 0 0 = -k r (P 0 0 ) ∂P 0 0 ∂z + 1 K 0 e 3 dans ]0, 1[×Ω (1.3.24)
et pour la composante lente de l'écoulement

     P 2 0 (t , x, z) = H 0 (t , x) -z dans ]0, 1[×Ω -div x K (H 0 ) ∇ x H 0 = - F 2 β - ∂ ∂t h soil h bot φ s(P 2 0 ) d z dans ]0,
                 φ s(P n ) -s(P n-1 ) ∆t + ∂u n 3 ∂z = 0 dans (x, z) ∈]h n , h soil [ u n 3 = -k r (P n ) K zz 1 ρ g ∂P n ∂z + 1 dans (x, z) ∈]h n , h soil [ α P n + β u n 3 = F n sur Γ soil P n x, h n (x) = H n (x) -h n (x) sur Ω x (1.4.1)            u n 3 h soil + 1 ∆t h soil (x) h bot (x) φs(P n )d z - h soil (x) h bot (x) φs(P n-1 )d z -div x K ( H n ) ∇ x H n = 0 dans Ω x K ( H n ) ∇ x H n • n = 0 sur ∂Ω x (1.4.2) P n (x, z) = H n (x) -z dans (x, z) ∈ Ω - h n (t ) (1.4.3) h n (x) = Q(x, H n (x)) dans Ω x . (1.4.4)
Le premier résultat concernant cette discrétisation implicite est qu'elle engendre un schéma qui conserve la masse. On renvoie au chapitre 3 pour une justification plus précise.

Une autre remarque concerne le choix lui-même d'un schéma implicite, qui de fait, demandera plus de travail numérique pour sa résolution. En fait c'est une grande difficulté de traiter l'équation (1.4.2) qui peut devenir mal posée lorsque ces trois premiers termes sont considérés constants par rapport à l'inconnue H n . Cette situation pouvant justement apparaître lorsque qu'un schéma explicite ou semi-implicite est considéré.

La première étape dans la résolution de problème (1.4.

                     φ ∆t s(P ) -s(P ) + ∂u 3 ∂z = 0 dans ]h, h soil (x)[ u 3 = -k r (P ) K zz ∂P ∂z + 1 dans ]h, h soil (x)[ α P + β u 3 = F pour z = h soil (x) P (t , x, z) = H (t , x) -z dans [h bot , h] h(t , x) = Q(x,
     P n (x, •), u n 3 (x), h n (x) solution de R x, H n (x), P n-1 (x, •) a.e. dans Ω x Θ x, H n (x), P n-1 (x, •) -div x K ( H n )∇ x H n = 0 dans Ω x K ( H ) ∇ x H n • n = 0 sur ∂Ω x ( 1
H n 0 = H n-1 -Hérédité : pour tout k ∈ N * , on pose (P n k , u n k 3 , h n k , H n k ) solution du problème linéaire            P n k (x, •), u n k 3 (x), h n k solution de R x, H n k-1 (x), P n-1 (x, •) a.e. dans Ω x Θ(x, H n k-1 , P n-1 ) + Λ(x, H n k-1 , P n-1 )( H n k -H n k-1 ) -div x K ( H n k-1 )∇ x H n k -div x K ( H n k-1 )( H n k -H n k-1 )∇ x H n k-1 = 0 sur Ω x K ( H n k-1 ) ∇ x H n k • n = 0 sur ∂Ω x (1.4.8) où K est la dérivée de K par rapport à P . Par construction, si cette suite (P n k , u n k 3 , h n k , H n k ) converge, sa limite sera précise- ment la solution de (1.4.7) (P n , u n , h n , H n ).
Le point clé de cette procédure réside dans le choix d'une fonction de stabilisation Λ = Λ(x, H , P ) convenable. On note que, par construction, si la régularité de Θ par rapport à H est suffisante, la procédure (1.4.8) est exactement la méthode de Newton en choisissant Λ(x, H , P ) = ∂Θ ∂ H . On propose dans le chapitre 3 deux expressions différentes pour Λ. Chacune de ces expression sera bien adaptée à chacune des deux situations suivante -Lorsque l'interface Γ h n'est pas contrainte par le niveau du sol ni par celui du fond de l'aquifère -La situation inverse où l'on a h n k-1 (x) ∈ {h bot (x), h soil (x)}. En effet la dépendance de Θ vis-a-vis de H diffère lorsque l'une où l'autre de ces situations à lieu.

Résultats numériques : comparaison entre les modèles

L'objectif de cette section est d'utiliser le schéma introduit dans la section précédente pour comparer numériquement le modèle original de Richards 3d (1.2.25) avec ceux de la classe de modèles (1.3.7)-(1.3.11) pour différents choix de valeurs R dans (1.3.5). Toutes ces comparaisons se baseront sur l'évolution de l'interface h sat (t , x) (définie ci-dessous) représentant le niveau de la nappe phréatique que l'on notera Ω - h sat (t ). Ces objets sont définis pour une pression P = P (t , x, z) donnée par :

h sat (t , x) := sup I t ,x , I t ,x := z ∈ [h bot (x), h max (x)] | P (t , x, z ) > P s , ∀z ∈ [h bot (x), z[ , ( 1
Ω - h sat (t ) := (x, z) ∈ Ω | z < h sat (t , x) . (1.4.10)
En particulier, on a par construction et si P est continue que

P t , x, h sat (t , x)      = 0 si h bot < h sat < h max ≥ 0 si h sat = h max ≤ 0 si h sat = h bot .
De plus P (t , x, z) ≥ 0 pour tout z ∈]h bot , h sat ]. On note également que le sol est complètement saturé dans Ω - h sat (t ) pour tout t ∈]0, T [. La fonction h sat (t , x) est l'isopression P (t , x, h sat (t , x)) = 0 si et seulement si l'aquifère ne déborde ni n'est vide à la position x ∈ Ω. Cette définition est telle que l'on a h(t , x) = h sat (t , x) lorsque R = 0 dans (1.3.14).

Dans chaque simulation nous allons nous restreindre à un aquifère 2d occupant un domaine Ω de type (1.2.17) avec Ω x =]0, L x [ pour L x > 0. On note encore h soil (x) et h bot (x) = h bot < 0 les fonctions caractérisant les frontières hautes et basses de l'aquifère. La valeur précise de L x sera modifiée afin de voir son influence 6sur l'écoulement.

Bien que de nombreux modèles existent pour décrire la saturation et la conductivité, nous nous restreignons dans cette thèse au modèle de Brooks et Corey [13] suivant. s(P ) = (P s /P ) λ , k r (P ) = (P s /P ) 2+3λ , (P s , λ) = (-1.5, 3),

ρ = 1, φ = 0.1, K 0 = 0.1 I 3 .
où I 3 est la matrice identité de taille 3. Nous présentons deux cas tests de recharge de la nappe phréatique. Le premier avec un réservoir d'eau déjà dans le milieu poreux mais au dessus de la position initiale de la nappe phréatique. Le second dans lequel l'apport en eau provient de la surface par l'intermédiaire d'une condition de Dirichlet au niveau du sol.

Dans chacun des cas, on considère une situation initiale à t = 0, pour laquelle la fonction h 0 sat dans (1.4.9) représente le niveau supérieur de la nappe phréatique. Ce niveau est de plus choisi constant h sat (0, x) = h 0 sat ∈]h bot , h soil [. La pression initiale correspondant est dans l'état stationnaire donné par P (0, x, z) = h sat (0, x) -z pour tout (x, z) ∈ Ω hormis près des réservoirs d'eau. Dans chacun des tests le niveau de la couche imperméable au fond de l'aquifère sera donné par h bot = -5. une rivière, un lac... On choisi un niveau du sol donné par

h soil (x) = - 3 2 e - L x 20 (x- L x 2 ) 2 .
(1.4.12)

Ce n'est pas l'objectif cette thèse de décrire précisément l'écoulement de l'eau de surface. On se concentre sur l'écoulement sous-terrain et on considérera un modèle de pression stationnaire pour l'eau de surface. On introduit pour cela la fonction h riv = h riv (x) qui caractérise le niveau supérieur de la rivière. L'eau de surface occupe donc la région

Ω riv := (x, z) ∈ Ω x × R | z ∈]h soil (x), h riv (x)[ . ( 1
Γ D soil = (x, z) ∈ Γ soil | h soil (x) ≤ h riv et Γ N soil = (x, z) ∈ Γ soil | h soil (x) > h riv .
La région sous la rivière est Γ D soil et on va y choisir une condition de Dirichlet pour la pression. En vu de l'équation (1.4.14) on choisit donc P (x, z) = h riv -h soil (x) sur Γ D soil . D'autre part, pour simplifier la modélisation on considère une condition de Neumann homogène sur Γ N soil . On représente la situation initiale dans le tracé de gauche de la Figure 1.3 Comme précédemment, l'échelle de gris dans ce tracé correspond à la saturation du sol, le plus sombre correspondant à s 1.

Dépendance par rapport à R. On compare ici les solutions obtenues à l'aide du modèle original de Richards 2d et celles obtenues avec le modèle couplé (1.3.7)- (1.3.11). En particulier on souhaite montrer l'influence du paramètre R sur l'écoulement.

On note h 2d sat l'iso-pression P = 0 obtenue du problème de Richards original. On note de plus pour tout R > 0, h R sat l'évolution de cette même iso-pression obtenue du problème couplé (

0 1 2 3 Case R = 0 Case R = 3 Case R = 7 Case R = ∞ FIGURE 1.5 -Erreur en espace et en temps 1 L x h 2d sat -h R sat L 1 ([0,T ]×Ω)
vs le ratio longueur/ profondeur de l'aquifère. En haut pour le cas test 1 et en bas pour le cas test 2.

Voici quelques commentaires sur les résultats obtenus (on renvoie au Chapitre 3 pour plus de détails) :

-Dans le cas R = 0, on a h = h 0 sat . L'écoulement vertical dans toute la zone saturée Ω - h sat = Ω - h est considérée comme étant instantanée. Lorsque l'eau provenant des réservoirs atteint la nappe phréatique, le flux (u • e 3 )| Γ h augmente rapidement. Ainsi la charge hydraulique H correspondante, obtenue de (1.3.10), évolue rapidement également. L'écoulement horizontal associé s'avère être plus rapide que celui de référence. Ceci continue durant toute l'expérience.

-Dans le cas R = +∞, on a h = h bot (voir (1.3.14)). L'écoulement vertical est décrit par l'équation de Richards 1d dans tout le domaine, même dans la partie saturée sous le niveau z = h ∞ sat . C'est la situation la plus éloignée de la situation précédente. Dans celle-ci, l'écoulement horizontal est plus lent (même trop lent) que celui de référence.

-Dans le cas R = 7, on a h bot ≤ h ≤ h 7 sat . Il s'agit d'un situation intermédiaire à celles ci-dessus. L'écoulement correspondant présente ainsi un comportement intermédiaire à ceux précédents. Le comportement de référence peut alors être mieux approché par ce type de modèle. La valeur optimale de R est cependant un problème encore ouvert. Pour conclure, nous somme ici dans des cas où le ratio ε =profondeur/largeur de l'aquifère vaut 5/40 = 1/8, ratio qui n'est donc pas très proche de 0. Ainsi, on note

ANALYSE MATHÉMATIQUE D'UN MODÈLE COUPLÉ DUPUIT-RICHARDS À INTERFACE LIBRE.35

que même pour ce type de ratio, le modèle couplé approche bien le problème de référence (cf Figure 1.4) pour un bon choix de R. On s'intéresse dans la partie suivante à l'évolution de l'erreur d'approximation lorsque ce ratio diminue.

Erreur totale en espace-temps v.s. ratio. On s'intéresse ici a quantifier l'erreur d'approximation du problème de Richards original par le problème couplé lorsque l'aquifère devient de plus en plus fin. D'après les résultats obtenus du développement asymptotique, on s'attend à ce que cette erreur diminue lorsque le ratio ε diminue.

Pour cela, on fixe la profondeur de l'aquifère à h bot = -5 et on fait varier L x dans [20,160]. On s'intéresse également à l'influence de R ∈ {0, 3, 7, ∞}. On mesure l'erreur entre les solutions par la valeur de la quantité 1 L x h 2d sat -h R sat L 1 ([0,T ]×Ω) . On trace cette erreur dans la Figure 1.5.

Comme attendu, l'erreur décroît lorsque le ratio largeur/profondeur augmente. Le choix optimal pour la fonction R étant dans ces situations obtenue pour R = 7. Il est à noté cependant que cette "optimalité" dépends du ratio dans le second cas, rendant plus difficile la détermination de ce choix a priori. Une autre remarque est que le choix R = 0 qui est classiquement utilisé lorsque l'hypothèse de Dupuit est faite, peut presque toujours être amélioré par l'utilisation d'une valeur R > 0.

ANALYSE MATHÉMATIQUE D'UN MODÈLE COUPLÉ DUPUIT-RICHARDS À INTERFACE LIBRE.

L'étude mathématique du modèle est particulièrement délicate déjà en raison de la présence de la frontière libre entre les deux sous-domaines. De plus, il existe une difficulté mathématique constante dans la structure du système des EDP modélisant la dynamique des eaux souterraines. En effet, dans le cas d'une nappe phréatique libre, nous devons faire face à la disparition progressive de l'eau dans la zone de désaturation et donc à la disparition d'une des principales inconnues du problème (même dans le modèle simplifié de Richards). Mais la difficulté principale réside certainement dans le couplage entre les deux zones qui s'exprime en général par des termes de flux à l'interface. La définition même de ces termes va nécessiter une régularité importante de la pression dans tout le domaine mais aussi de la fonction décrivant la profondeur de l'interface. Pour ces raisons, pour simplifier quelque peu l'analyse mathématique du problème, nous supposerons tout d'abord, que la contribution en eau de la zone de désaturation est prise en compte par un terme source variable (en temps et en espace) et non par le terme flux impliquant directement la pression définie dans la zone de désaturation. Quoiqu'il en soit, même avec cette approximation, nous aurons toujours besoin d'établir des résultats forts pour la régularité de la pression. Cette régularité découlera des hypothèses faites sur les paramètres caractérisant l'équation de Richards Il existe une littérature abondante concernant les équations classiques de Richards. Mentionnons les travaux incontournables de Alt et al ([4, 5]) ainsi que les articles de [14,25,[START_REF] Showalter | Partially saturated flow in a poroelastic medium[END_REF] dédiés à l'étude de l'équation "dégénérée" en temps ∂ t θ(p) -∆p = 0, où θ(p) désigne la teneur volumétrique en humidité. Citons également dans le cas unidimensionnel le travail de Yin ([51]) concernant l'existence d'une solution faible pour le problème totalement dégénéré

∂ t θ(p) -∂ x (κ(θ(p))∂ x p) = 0, où l'auteur suppose que θ et κ > 0.
Dans le contexte des problèmes de transport réactifs, impliquant l'équation Richards 3D couplée à une équation hyperbolique, on peut citer le papier de Choquet [16] où la saturation et la mobilité sont fortement couplées par la pression mais aussi les travaux antérieurs de Amirat et al [6] dans lequel le couplage est moindre, plus précisément où on a simplement que θ = θ(x) and κ = κ(x).

Le modèle (M ) que nous allons considérer ici est donc celui couplant l'équation de Richards 3D (pour la description de l'écoulement dans la frange capillaire) à l'équation obtenue après une moyennisation verticale dans la partie saturée de l'aquifère. Plus précisément, -Dans Ω + h (t ) on a l'équation de Richards 3D suivante

           ∂ t θ(P ) + θα P ∂ t P + ∇ • v = Q s dans ]0, T [×Ω + h (t ), v • n = 0 sur ]0, T [× Γ soil ∪ Γ ver , P t , x, h(t , x) = P s dans ]0, T [×Ω x , P (0, x, z) = P 0 (x, z) dans Ω + h (0). la vitesse effective v étant donnée par v = -K ∇( P ρ o g + z), K = κ(θ(P )) K 0 ρ 0 g µ .
-Dans Ω - h (t ), la pression P satisfait

P (t , x, z) = ρ 0 g P s ρ 0 g + h -z dans ]0, T [×Ω - h (t ).
-La profondeur de l'interface Γ h , h, vérifie dans

Ω x      S 0 B f ∂ t h -∇ • (B f K ∇ h) = B f Q -v |z=h + • ∇(z -h) dans ]0, T [×Ω x , K ∇ h • n = 0 sur (0, T ) × ∂Ω x , h(0, x) = h 0 (x) dans Ω x .

ANALYSE MATHÉMATIQUE D'UN MODÈLE COUPLÉ DUPUIT-RICHARDS À INTERFACE LIBRE.37

Nous rappelons que Ω est un ouvert borné de R 3 et que Ω x correspond à la projection de Ω sur le plan horizontal. On désigne par n la normale unitaire extérieure pointant à l'extérieur de Ω. La frontière de Ω, supposée C 1 , est notée par Γ et Γ = Γ soi l ∪ Γ bot ∪ Γ ver . L'intervalle de temps est (0, T ), T étant un réel positif et on posera Ω T = (0, T ) × Ω. Pour des raisons de simplicité, on supposera que la topographie est définie par une fonction constante telle que h bot ∈ R.

Nous allons à présent introduire quelques notations et rappeler des résultats généraux mathématiques utiles pour la suite de l'étude.

Résultats auxiliaires.

Soit Ω un ouvert borné de R 3 . Pour des raisons de brièveté, nous noterons

H 1 (Ω ) = W 1,2 (Ω ) et V = H 1 0 (Ω ), V = H -1 (Ω ), H = L 2 (Ω ).

Nous rappelons que les injections

V ⊂ H = H ⊂ V sont denses et compactes. Pour tout T > 0, soit W (0, T, Ω ) désignant l'espace W (0, T, Ω ) := ω ∈ L 2 (0, T ;V ), ∂ t ω ∈ L 2 (0, T ;V ) muni de la norme Hilbertiennee • W (0,T,Ω ) = • 2 L 2 (0,T ;V ) + ∂ t • 2 L 2 (0,T ;V )
W (0, T, Ω ) ⊂ C ([0, T ]; [V,V ] 1 2 ) = C ([0, T ]; H ) tandis que l'injection W (0, T, Ω ) ⊂ L 2 (0, T ; H ) (1.5.1)
est compacte (lemme d'Aubin, voir [START_REF] Simon | Compact sets in the space lp(o,t ; b)[END_REF]).

Il sera utile d'introduire l'espace

X (0, T, Ω ) = L ∞ (0, T ; H 1 (Ω )) ∩ H 1 (0, T ; L 2 (Ω )) muni de la norme u X (0,T,Ω ) = u L ∞ (0,T ;H 1 (Ω )) + u H 1 (0,T ;L 2 (Ω )) .
On utilisera le résultat suivant dû à F. Mignot (see [23]).

Lemme 1.5.1 Soit f : R → R une application croissante et continue telle que

lim sup |λ|→+∞ f (λ)/λ < +∞. Soit ω ∈ L 2 (0, T ; H ) tel que ∂ t ω ∈ L 2 (0, T ;V ) et f (ω) ∈ L 2 (0, T ;V ). Alors 〈∂ t ω, f (ω)〉 V ,V = d d t Ω ω(•,y) 0 f (r ) d r d y in D (0, T ). Ainsi, pour tout 0 ≤ t 1 < t 2 ≤ T t 2 t 1 〈∂ t ω, f (ω)〉 V ,V d t = Ω ω(t 2 ,y) ω(t 1 ,y) f (r ) d r d y.
Nous allons aussi rappeler un lemme donnant un résultat de régularité dans le cas parabolique. On définit X p = L p (0, T ;W 1,p 0 (Ω)), muni de la norme

( T 0 ||v(t )|| p W 1,p 0 (Ω) d t ) 1/p = ||∇v|| L p (Ω T ) n . On introduit Y p = L p (0, T ;W -1,p (Ω)). Nous soulignons que l'application v → d i v x v envoie (L p (Ω T )) n dans L p (0, T ;W -1,p (Ω)
). On peut alors énoncer le résultat suivant (cf. [34]) :

Lemme 1.5.2 Soit P = ∂ ∂t -∆, l'opérateur associé aux conditions aux limites de Dirichlet homogènes. Alors, étant donné F ∈ Y p , il existe une unique solution u ∈ X p telle que :

Pu = F i n Ω T , u(0) = u 0 . De plus u X p ≤ ĝ (p) F Y p + u 0 W 1,p 0 (Ω) , (1.5.2
)

où ĝ (p) = ||P -1 || L (Y p ;X p ) .
Nous rappelons que ĝ (2) = 1.

Résultats principaux

Nous souhaitons donner un résultat d'existence de solutions faibles, physiquement admissibles pour le modèle (M ) complété par les conditions initiales et les conditions aux limites.

Nous introduisons les fonctions x

+ := sup(0, x) et T l définies par T l (u) = (u -h bot ) ∀u ∈ (h bot , h soi l ).
La fonction T l est étendue continûment et de manière constante en dehors de l'intervalle (h bot , h soi l ). T l (h) représente l'épaisseur de la zone saturée d'eau douce dans le réservoir. Nous soulignons que la fonction T l agit aussi sur le terme source Q empêchant ainsi le pompage dans les régions où il n'y aurait pas d'eau.

Posant

Ω + h = Ω x × (h, h soi l ) et Ω - h = Ω x × (h bot , h
), on considère alors le système suivant :

∂ t θ(P ) + θ(P )α P ∂ t P + ∇ • v = Q s , v = -K (θ(P ))∇( P ρ o g + z), dans (0, T ) × Ω + h , (1.5.3) P (t , x, z) = ρ 0 g P s ρ 0 g + h -z dans (0, T ) × Ω - h , (1.5.4
)

φ∂ t h -∇ • T l (h) K ∇ h = T l (h) Q dans (0, T ) × Ω x . ( 1
T δ = T l + δ et l'équation régularisée suivante φ∂ t h δ -∇ • T δ (h δ ) K ∇ h δ = T l (h δ ) Q dans (0, T ) × Ω x . (1.5.6)
Le système (1.5.3)-(1.5.5) est complété par les conditions aux limites et les conditions initiales suivantes :

P Γ h = P s dans (0, T ), ∇P • n = 0 sur (0, T ) × (Γ soi l ∪ Γ ver ), P (0, x, z) = P 0 (x, z) dans Ω + h 0 . (1.5.7) ∇h • n = 0 sur (0, T ) × ∂Ω x , h(0, x) = h 0 (x) dans Ω T , (1.5.8)
La fonction P s est supposée constante par rapport au temps et à l'espace. La fonction P 0 ∈ H 2 (Ω) satisfait la condition de compatibilité

P 0 (x, h 0 ) = P s in Ω + h 0 .
On suppose aussi que h 0 ∈ L ∞ (Ω x ) est telle que h 0 ≥ h bot a.e. dans Ω x . Enfin, le terme source Q est une fonction donnée de l'espace L 2 (0, T ; H ).

Nous allons maintenant détailler les hypothèses mathématiques. Nous commençons par les caractéristiques de la structure poreuse. Nous limitons notre étude au cas isotropique ainsi le tenseur K 0 est supposé être un scalaire. Dans la partie saturée de l'aquifère, la conductivité hydraulique K est alors égale à la constante K 0 ρ 0 g µ . A partir de maintenant, nous noterons la masse volumique ρ 0 simplement ρ. Les fonctions θ and κ dépendent de la pression et nous supposons que

θ ∈ C 1 (R), 0 ≤ θ(x) ≤ θ + , θ (x) ≥ 0 ∀x ∈ R, (1.5.9) κ ∈ C (R), 0 ≤ κ(x) ≤ κ + ∀x ∈ R + . (1.5.10)
Pour le précédent système parabolique, nous énonçons et prouvons le résultat suivant Théorème 1.5.4 Supposons qu'il existent deux nombres réels θ -and κ -tels que

θ(x) ≥ θ -> 0 ∀x ∈ R, κ(x) ≥ κ -> 0 ∀x ∈ R + . ( 1
∈ L ∞ (0, T ; H 1 (Ω)) ∩ H 1 (0, T ; L 2 (Ω)) est solution de (1.5.3)-(1.5.7) et (1.5.4) ; (b) la fonction h ∈ L ∞ (0, T ; H 1 (Ω)) ∩ H 1 (0, T ; L 2 (Ω)) est solution de (1.5.6)-(1.5.8).
Remarque 1.5.5 Comme mentionné précédemment, l'équation (1.5.5) devient dégénérée lorsque l'interface touche le bas de l'aquifère. Pour la première étape de la preuve, il est possible de surmonter cette dégénérescence en introduisant une fonctionnelle "entropie". Mais pour établir plus de régularité pour h, nous devons alors supposer que l'épaisseur d'eau douce δ dans l'aquifère reste toujours strictement positive i.e. (hh bot ) ≥ δ > 0. Cela fournit une interprétation du coefficient de diffusion δ. Une autre interpretation de δ est de le voir comme l'épaisseur de l'interface entre la zone saturée et la zone insaturée ainsi que cela est fait dans [15], cela implique alors qu'il ne s'agit plus d'une interface nette.

Les fonctions θ et κ caractérisent le type mathématique du problème. Plus précisément, le problème (1.5.3)-(1.5.5) est de type parabolique si θ et κ sont des fonctions strictement positives et de type parabolique dégénéré si les fonctions θ and κ sont seulement supposées positives. Le second résultat de ce chapître est consacré au cadre parabolique du système. Dans ce cas, nous prouvons le résultat suivant Théorème 1.5.6 Supposons que les fonctions θ et κ sont positives et qu'elles satisfassent (1.5.9)-(1.5.10). Supposons de plus que il existe 0 > 0 tel que κ soit croissante sur (0, 0 ).

(1.5.12)

Alors le système (1.5.3)-(1.5.7), (1.5.4), (1.5.6)-(1.5.8) admet une solution faible (P, h) satisfaisant

1.5. ANALYSE MATHÉMATIQUE D'UN MODÈLE COUPLÉ DUPUIT-RICHARDS À INTERFACE LIBRE.41 (a) la fonction P ∈ W (0, T, Ω) ∩ L 2 (Ω T ) est solution de (1.5.3)-(1.5.7) dans l'espace L 2 (0, T ; H -1 (Ω)) et la vitesse de Darcy v ∈ L 2 (Ω T ) 3 ; (b) la fonction h ∈ L 2 (0, T ; (H 1 (Ω)) )∩L 2 (0, T ; H 1 (Ω)) est solution de (1.5.5)-(1.5.8).
La section suivante est dédiée à la preuve du Théorème 1.5.4.

Trame de la preuve du Théorème 1.5.4

Ainsi que nous l'avons mentionné dans l'introduction, le problème est caractérisé par la présence d'une interface libre entre les deux domaines, par les difficultés inhérentes aux équations de Richards et par le couplage entre les deux équations. En outre, nous devons faire face à la disparition progressive de l'eau dans la zone de désaturation et, donc à la disparition de l'une des principales inconnues du problème. Comme le système est fortement couplé, nous appliquons une approche par point fixe pour le résoudre. La clé est de résoudre d'abord l'équation en h. Cela est possible grâce à la continuité de la composante normale du flux de Darcy à travers l'interface qui permet d'exprimer ce flux en fonction de la dérivée temporelle de la profondeur de l'interface. Nous devons ensuite obtenir un résultat de régularité suffisant pour h et sa dérivée temporelle afin de "linéariser" (en quelque sorte) l'équation en pression. Cette régularité peut être obtenue grâce à une régularisation qui garantit une épaisseur d'eau douce toujours strictement positive dans l'aquifère δ > 0. Pour surmonter la difficulté liée aux fortes non-linéarités dans l'équation en pression, nous effectuons un changement de variable. Plus précisément, nous utilisons la transformation fondamentale de Kirchoff qui permet de linéariser le terme correspondant à la divergence de l'équation (1.5.3) sur un domaine variable dépendant de l'interface h calculée à l'étape précédente. Nous établissons enfin des estimations uniformes suffisantes pour la pression sur l'ensemble du domaine. En utilisant le théorème du point fixe de Schauder, nous prouvons le résultat d'existence pour le problème complet.

Nous énonçons ci-dessous le lemme correspondant au premier résultat d'existence pour h.

Lemme 1.5.7 Soit h

0 ∈ H 1 (Ω x ) et Q ∈ L 2 (0, T ; Ω x ), il existe une fonction h ∈ W (0, T, Ω x ) solution de (1.5.5)-(1.5.8) qui satisfait h L 2 (0,T ;H 1 (Ω x )) ≤ M et h L 2 (0,T ;(H 1 (Ω x )) ) ≤ M ,

où M and M dépendent seulement des données du problème. De plus, le principe du maximum suivant est vérifié

h bot ≤ h(t , x) p.p. x ∈ Ω x et pour tout t ∈ (0, T ).
Esquissons notre stratégie. La première étape consiste à utiliser un théorème du point fixe de Schauder pour prouver un résultat d'existence pour un problème régularisé auxiliaire. Plus précisément nous régularisons la fonctionT l avec le paramètre δ > 0. Une difficulté est que l'application utilisée pour l'approche par point fixe doit être continue dans L 2 (0, T ; H 1 (Ω x )). Nous montrons ensuite que la solution régularisée satisfait le principe maximum annoncé dans le lemme1.5.7. Nous établissons enfin une estimation uniforme suffisante (grâce à une fonctionnelle "entropie", cf. [3]) puis nous faisons tendre la régularisation δ vers zéro. Afin de traiter le couplage avec la zone non saturée, nous avons besoin de plus de régularité pour la profondeur de l'interface h plus spécialement pour la dérivée en temps de h. Nous ne pouvons pas obtenir cette régularité dans le cas où l'interface de saturation h touche le fond de l'aquifère, ce qui correspond au cas précédent dégénéré. Nous considérons donc le problème régularisé (1.5.6)-(1.5.8) au lieu du problème d'origine (1.5.5)-(1.5.8). Il est évident que l'existence de la solution pour ce problème régularisé résulte du lemme précédent. 2 et la constante M u dépend seulement des données du problème.

Lemme 1.5.8 Soit δ > 0 et h ∈ W (0, T, Ω x ) une solution de (1.5.6)-(1.5.8), elle satisfait u L ∞ (0,T ;H 1 (Ω)) + u H 1 (0,T ;L 2 (Ω)) ≤ M u , (1.5.13) où u = ((h -h bot )+δ)
Remarque 1.5.9 Si on prend en compte la contribution de l'eau provenant de la zone supérieure de l'aquifère, le terme source peut être exprimé par le flux v |z=h + .∇(z -h) où v |z=h + est le flux de Darcy dans la frange capillaire. Nous devons donc estimer la norme L 2 de ce flux, ce qui représente la principale difficulté de l'analyse mathématique. Cela implique en effet d' estimer la norme L 4 de v, ce qui pourrait être fait en appliquant le lemme 1.5.2 à une linéarisation de l'équation de Richards.

Les hypothèses (1.5.9)-(1.5.10) sont suffisantes pour définir la fonction primitive B telle que

B (P ) = P κ(θ(P )) (1.5.14)
L'application B est bijective par (1.5.11) et l'existence de p telle que p = B (P ) est équivalente à l'existence de P solution de (1.5.3). En appliquant la transformation de Kirchoff à l'équation (1.5.3), nous considérons maintenant "le problème transformé" dans la frange capillaire supérieure

τ(p)∂ t p -K ∆p -∇ • ρg K κ(θ(B -1 (p)) e 3 = 0 dans (0, T ) × Ω + h , (1.5.15) p Γ h = B (P s ) dans (0, T ), ∇p • n = 0 sur (0, T ) × (Γ soi l ∪ Γ ver ), p(0, x, z) = B (P 0 ) dans Ω + h 0 , (1.5.16) 1.5. ANALYSE MATHÉMATIQUE D'UN MODÈLE COUPLÉ DUPUIT-RICHARDS À INTERFACE LIBRE.43 où τ(p) = θ + α P θ B -1 (p) B -1 (p). Notez qu'il existe un réel non négatif τ -tel que 0 < τ -:= α P θ - κ + ≤ τ(p) ≤ τ + := α P θ + κ - .
(1.5.17)

Nous construisons à présent le cadre pour appliquer le Théroème du point fixe de Schauder (see [21,[START_REF] Zeidler | Nonlinear functional analysis and its applications, part 1[END_REF]). Pour la stratégie du point fixe, nous introduisons un sousensemble convexe K p de W (0, T, Ω) . Nous posons

K p = {v ∈ X (0, T, Ω); v X (0,T,Ω) ≤ M p }, la constante M p étant définie plus tard. Soit p ∈ K p et h ∈ W (0, T, Ω) une solution (1.5.6)-(1.5.8), nous résolvons le problème suivant p(t , x, z) = B ρ g P s ρg + h -z dans (0, T ) × Ω - h , (1.5.18) τ( p)∂ t p -K ∆p -ρg K κoθoB -1 ( p)∂ z p = 0 dans (0, T ) × Ω + h , (1.5.19
)

p Γ h = B (P s ) dans (0, T ), ∇p • n = 0 sur (0, T ) × (Γ soi l ∪ Γ ver ), p(0, x, z) = B (P 0 ) dans Ω + h(t =0) , (1.5.20)
On peut alors énoncer et établir le lemme.

Lemme 1.5.10 Soit δ > 0, h ∈ W (0, T, Ω x ) une solution de (1.5.6)-(1.5.8) et p ∈ K p . Il existe une fonction unique p ∈ W (0, T, Ω) solution de (1.5.19)-(1.5.20) satisfaisant p L ∞ (0,T ;H 1 (Ω)) + p H 1 (0,T ;L 2 (Ω)) ≤ M p , (1.5.21)
où la constante M p dépend seulement des données du problème.

Ce lemme permet de définir une application de K p dans lui-même, continue pour la norme L 2 (0, T ; H 1 (Ω)). Les estimations du Lemme 1.5.10 et le Théorème de Schauder nous permettent de conclure la démonstration du Théorème 1.5.4.

Trame de la preuve du Théorème 1.5.6

Nous visons maintenant à établir le résultat d'existence pour le problème dégénéré du Théorème 1.5.6. Dégénéré signifie que nous ne supposons plus l'existence d'une saturation résiduelle strictement positive dans la zone de désaturation. La saturation peut donc être nulle dans certaines zones de l'aquifère. Du Théorème 1.5.4, on peut affirmer qu'il existe une solution faible (p ε , h ε ) ∈ X (0, T, Ω) 2 du problème parabolique suivant, pour tout > 0

∂ t θ (P ) + θ (P )α P ∂ t P + ∇ • v = 0, v = -κ (θ (P ))∇( P ρg + z), dans (0, T ) × Ω + h , (1.5.22) P (t , x, z) = ρ 0 g P s ρ 0 g + h -z dans (0, T ) × Ω - h , (1.5.23) φ∂ t h -∇ • (T l (h ) + δ) K ∇ h = T l (h ) Q dans (0, T ) × Ω x . (1.5.24)
Le système (1.5.3)-(1.5.5) est complété par les conditions aux limites et les conditions initiales suivantes :

∇h • n = 0 sur (0, T ) × ∂Ω x , h (0, x) = h 0 (x) dans Ω T , (1.5.25
)

P Γ h = P s dans (0, T ), ∇P • n = 0 sur (0, T ) × (Γ soi l ∪ Γ ver ), P (0, x, z) = P 0 (x, z) dans Ω + h 0 , (1.5.26) où θ = θ + et κ = κ + , (1.5.27)
Les fonctions θ et κ satisfont (1.5.9) et (1.5.10).

En suivant l'idée développée dans [16], on prouve qu'il existe une sous-suite extraite de solutions du problème régularisé (1.5.22)-(1.5.26) qui converge faiblement vers une solution du problème original. Nous rappelons d'abord des estimations uniformes. En raison de (1.5.12), la conductivité κ est une fonction croissante sur (0, ) dès que < 0 . Les estimations établies pour p dans la section 1.5.2 écrites pour P deviennent

θ(P ) + θ (P ) + ε + κ θ (P ) + ∂ t P L 2 2(Ω T ) ≤ C (1.5.28) κ θ (P ) + ∇P (L ∞ (0,T ;L 2 2(Ω))) 3 ≤ C . ( 1
H (q) = q θ(s)(κoθ)(s) d s.
(1.5.30)

Au vu de (1.5.28)-(1.5.29), la fonction H (P ) est uniformément bornée dans H 1 (Ω T ). Nous définissons la fonction limite H telle que

H (P ) → H dans L 2 (Ω T ) et p.p. dans Ω T .
Nous définissons ensuite la fonction de troncature T P d par

T P d (x) = x si x ≥ P d , P d si x < P d .
(1.5.31)

ANALYSE MATHÉMATIQUE D'UN MODÈLE COUPLÉ DUPUIT-RICHARDS À INTERFACE LIBRE.45

Soit t = T P d (P ).

(1.5.32)

Par définition de θ, H et T P d , on remarque que

H (P ) = H (T P d (P )) = H (t ε )
et le résultat de convergence ci-dessus s'écrit

H (t ) → H dans L 2 (Ω T ) et p.p. dans Ω T .
La fonction H est bien sûr non bijective sur (-∞, P d ). On considère alors une extension continue et bijective H de

H |(P d ,∞) → R. Posant b = H -1 ( H ), (1.5.33) on a H (t ) = H (t ) → H = H (b) dans L 2 (Ω T ).
La fonction H étant continue et bijective, nous concluons que

t = T P d (P ε ) → b in L 2 (Ω T ) and a.e. in Ω T .
Les passages à la limite dans les équations (1.5.22)-(1.5.26) sont alors immédiats.

Bibliographie 1. INTRODUCTION

Contamination of soil and groundwater is a major concern that affects all populated areas. Many works are thus developed for studying the vulnerability of aquifers with regard to agricultural, industrial, or sewage pollutions. There is an abundant literature on each of the involved processes (geological, physical, chemical...), so that we can consider that the corresponding model is already available. Nevertheless there is a so wide variety of processes (chemical, hydrogeological, anthropic) acting in a so wide range of temporal and geometrical length scales that the assembly of the corresponding model bricks, if considered like toolboxes of a software, is, at best, computationally expensive.

In this multi-scale context, a particularly interesting issue is a proper and tractable model for the exchanges between the overland and the underground waters. Indeed, the challenge consists in capturing very different physical phenomena, the fast and essentially vertical leakage coming from the surface through an unsaturated soil and the slow and essentially horizontal displacement in the saturated part of the aquifer, that are classically modelled by mathematical systems with very different structures. The question is all the more important that an accurate study of the interaction between the water table and the overland water is essential for many concerns, concerns that disallow the use of classical time upscaling processes. It is in particular crucial for studying the transport of chemical components in the aquifer. Indeed, it turns out that many chemical reactions occur in the first meters of the subsoil, where oxygen is still very present. As a byproduct, the chemical species that reach the water table are not necessarily the same than those that have left the surface, and there is a large range of kinetics reaction times to handle with. There is actually no scale separation.

In the present paper, we focus on the hydrogeological question. We thus consider the displacement of a wetting phase (water) in the presence of a non-wetting fluid (air) in a porous medium. Assuming that the air present in the unsaturated zone has infinite mobility allows to use a model for immiscible fluid flow simplified by the Richards hypothesis. The saturation is thus considered as a monotone function depending of the pressure head and the so-called Richards model consists in a nonlinear three-dimensional equation of degenerate parabolic type. All the existing simplified models for the fluid displacement in aquifers are motivated by the characteristics of the flow in their saturated part. A form of stratification enables the definition of interfaces and the slowness of the natural dynamics ensures that these interfaces have a smooth and stable behaviour. Moreover the flows are essentially orthogonal to the walls (Dupuit's hypothesis). These points allow the vertical integration of the Richards equation in the saturated area and lead to the use of a family of 2D models developed since the 60's (see e.g. the works of Jacob Bear, [5,6]). A main weakness of the approach by vertical integration lies in its justification. It is only valuable for very precise length and time scales, the time scale in particular being completely different of the typical durations of chemical reactions (see once again [5] for empirical and qualitative arguments, see [13] for asymptotic computations). However, such 2D models are widely used, even out of their validity range and even if it turns out to be especially difficult to properly couple them with the flow in the unsaturated part of the underground. Only numerical attempts were done in this direction. We mention [11] where the integrated model is directly coupled with a surface model (see also the references therein). The unsaturated area of the aquifer is taken into account in [16] using a 1D-Richards equation coupled with a simplified model in the saturated part. However, the study is purely numerical and the model is not mathematically justified. In [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF], the latter kind of model is integrated into a computational code called "SHE" (for "European Hydrological System" and later became SHETRAN) in the case where the water table remains away from the ground level. See also [21], [14].

To the best of our knowledge, there is no mathematical justification for any "Dupuit-Richards" model specifying the hypotheses as well as the scales that allow its derivation from a more complete model (such as the 3D-Richards one).

Notice finally that the coupling of the surface and underground flows turns out to be more tractable when handling with a Richards equation (see e.g. [18] or [2] and [3] where the surface behavior is reduced to a Signorini boundary condition).

The goal of this work is to provide a simple model exploiting the low thickness of a confined or unconfined aquifer. In summary it consists in coupling purely vertical models (describing the flow at a small time scale) with a horizontal model (describing the flow at a long time scale). Clearly, given its construction, the model is simpler to manipulate numerically since the original 3D problem is replaced by the coupling of a 2D problem with several independent 1D-problems (which can be solved in parallel). Significant time savings are expected in the numerical processing.

This work could be viewed as another attempt using the numerically pragmatical methodology of [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF] and leading to a "Dupuit-Richards" model. Yet, our approach is quite different. First, we actually derive a class of models, each of them being characterised by the definition of some virtual interface which does not necessarily coincide with the water table (especially when trying to optimize the error). It follows that a model of this class does not necessarily contains a Dupuit component. The position of the virtual interface may even be an unknown of our model. Next, we aim at describing the flow in a large range of time scales, and, more precisely without any assumption of scale separation. The idea consists in always capturing both the fast and slow components of the flow given by Richards 3D equations, whatever the time scale. Their coupling is done through flux terms ensuring that the model is mass conservative (and thus avoiding the criticism done in [19]). Finally, the large validity range of the new class of models is justified by an asymptotic study. But, as already mentioned, no time scale separation is assumed in the present paper so that we adopt a new methodology for the asymptotic arguments. Let ε > 0 describe the ratio of the aquifer's deepness over its characteristic horizontal length. Assume that ε is small. The usual approach would consist in choosing a reference time for the study, introducing an asymptotic expansion of the solution of the 3D-Richards system and using the scale separation for identifying the equations governing the main order terms of this ansatz. This is the classical process for deriving an effective model. Here the asymptotic 2 analysis is not used for deriving an effective model for a given reference time. Rather, it is used for proving that each model of our new class and the 3D-Richards equation are associated with the same effective problem for any time scale. Basically:

1. At short times, the horizontal flow is very small and the vertical one satisfies a 1D-Richards problem. 2. At non-short times, the vertical flow appears instantaneous. The corresponding pressure profile satisfies the stationary 1D-Richards problem. Then the hydraulic head H does not depend on the vertical variable z. This corresponds to the so-called Dupuit hypothesis. 3. At large times, the horizontal flux is non-zero. It is ruled by a 2D-horizontal diffusion equation where the conductivity is the vertical average of the permeability tensor on the whole depth of the aquifer.

The paper is organised as follows: In Section 2, we describe the geometry of the problem, the physical parameters and unknowns. The classical 3D-Richards model is recalled. The main result and numerical simulations are given in Section 3. Namely, we present the systems coupling the vertical and the horizontal flows and we comment on the model. Finally, the formal asymptotic analysis of our models and of the 3D-Richards model are performed and compared in Section 4.

DESCRIPTION OF THE PROBLEM

This section is devoted to the description of the domain of study, of the physical parameters and of the unknowns which are chosen for characterising the flow through the Richards model.

Geometry

The aquifer corresponds to a cylindrical domain Ω ⊂ R 3 . For the sake of the simplicity, we assume vertical walls. The projection of Ω on any horizontal plane is an open domain Ω x ⊂ R 2 with boundary ∂Ω x . The lower and upper bases of Ω are respectively the graphs of real-valued functions h bot and h soil such that

h soil (x) > h bot (x) , ∀x ∈ Ω x . (2.1)
In summary the domain is given by:

Ω = (x, z) ∈ Ω x × R | z ∈ h bot (x), h soil (x) . (2.2)
We split the boundary ∂Ω of Ω in three parts (bottom, top and vertical)

∂Ω = Γ bot Γ soil Γ ver , Γ bot := (x, z) ∈ Ω | z = h bot (x) , Γ soil := (x, z) ∈ Ω | z = h soil (x) , Γ ver := (x, z) ∈ Ω | x ∈ ∂Ω x .
In the present paper, as already mentioned, we derive a class of models that are characterised by the position h of some virtual interface in the reservoir. For our construction, this function has to take its values in the semi-open interval [h bot , h soil ). For numerical implementation, an easy recipe consists in replacing the condition h < h soil by h ≤ h soilδ where δ is an arbitrary small positive real number. We thus introduce the auxiliary function h max defined by

h max = h soil -δ, 0 < δ 1. (2.3)

Three-dimensional Richards equation

We aim at deriving alternatives to the Richards equation. Let us briefly describe this classical model. In this paper we limit our study to a one-phase incompressible fluid which accordingly admits a constant density ρ ∈ R * + . First, in multiphase systems, observations have shown that an increase of the saturation of the non-wetting phase leads to an increase of the capillary pressure. The Richards model is moreover based on the assumption that the air pressure in the underground equals the atmospheric pressure, thus is not an unknown of the problem. One thus assumes that the saturation and the relative conductivity of the soil are given as functions of the fluid pressure P , 3 denoted respectively by s = s(P ) and k r = k r (P ). There is a large choice of available models for s and k r . The most classical examples for an air-water system are the van Genuchten model [20], with no-explicit dependance on the bubbling pressure but with fitting parameters, and the Brooks and Corey model [4], that we use in the simulations below:

+ x z z = h bot (x) z = h soil (x) z = h max (x)
s(P ) = (P s /P ) λ if P < P s 1 if P ≥ P s , k r (P ) = (P s /P ) γ if P < P s 1 if P ≥ P s , (2.4)
where λ > 0, γ = 2 + 3λ and P s < 0. Notice that our model would easily adapt to hysteretic soil properties ( [15], [17]). Since these methods, as of today, do not permit three-dimensional calculations, we guess that our 1D-2D models are even more interesting for their implementation than the 3D-Richards model. The important point is that these models are such that s(P ) = 1 ⇐⇒ P ≥ P s and k r (P ) = 1 ⇐⇒ P ≥ P s .

(2.5)

In particular, the water pressure is greater than the bubbling pressure P s if and only if the soil is completely saturated (P s being a fixed real number). The graphs of the functions s and k r given by the Brooks-Corey model used below for the numerical simulations are represented in Figure 2 (the parameters are given at the beginning of Subsection 3.4).

The soil transmission properties are characterised by the porosity function, φ = φ(x, z) ∈ (0, 1), and the permeability tensor, K 0 (x, z). The latter is a 3 × 3 symmetric positive definite tensor which describes the conductivity of the saturated soil at the position (x, z) ∈ Ω. We introduce

K xx ∈ M 22 (R), K zz ∈ R * and K xz ∈ M 21 (R) such that K 0 = K xx K xz K T xz K zz . (2.6) 4
The fluid is characterised by its pressure P and its velocity v solving the following Richards problem:

                 φ ∂s(P ) ∂t + div(v) = 0 in ]0, T [×Ω v = -k r (P ) K 0 1 ρg ∇P + e 3 in ]0, T [×Ω α P + β v • n = F on ]0, T [×Γ soil v • n = 0 on ]0, T [×(Γ bot ∪ Γ ver ) (2.7)
where g is the gravity constant and e 3 is the unitary vertical vector pointing up. The first equation describes the mass conservation of the constant density fluid in the case of an incompressible soil. The second equation is the Darcy's law associated with the nonlinear anisotropic conductivity k r (P ) K 0 . The boundary condition v • n = 0 on Γ bot corresponds to the impermeable layer at the bottom of the aquifer. The same is assumed on Γ ver to simplify the presentation. The condition at the soil level Γ soil is a Robin condition associated with given (α, β) ∈ (R + ) 2 \{0, 0} and F : Γ soil → R.

Remark 1 (Dominant behaviors in a shallow aquifer).

In Section 4 we investigate the behavior of the flow described by the 3D-Richards equations in the case of a thin aquifer and for various time scales. Let us summarise the conclusions of this asymptotic analysis. They might shed light on the comments about our models in the next section.

1. At any time scale, the dominant flow is the one in the vertical direction (see for example (4.9) in which the horizontal diffusion term appears multiplied by the small parameter ε). 2. In the short time scale (T ∼ 1), the horizontal flow is very small and the vertical one solves a classical 1D-Richards problem. 3. In non-short time scales (T ∼ ε -1 or T ∼ ε -2 ), the vertical flow appears as being instantaneous. The corresponding pressure profile satisfies a stationary 1D-Richards problem. Then the pressure is

P = ρ g (H -z)
where the hydraulic head H does not depend on the vertical variable z. The velocity is horizontal. This corresponds to the so-called Dupuit hypothesis. 4. In the long time scale (T ∼ ε -2 ), the horizontal flow is non-zero and it is ruled by a 2D-horizontal diffusion equation where the conductivity is the vertical average of the permeability tensor on the whole depth of the aquifer, from h bot to h soil .

MAIN RESULT AND NUMERICAL SIMULATIONS

Models coupling vertical 1d-Richards flow and Dupuit horizontal flow

Each of our models splits the description of the flow into two subregions of Ω (possibly time-dependent). These zones are defined by a function h = h(t , x) such that h bot ≤ h < h soil :

Ω - h (t ) := (x, z) ∈ Ω | z < h(x, t ) and Ω + h (t ) := (x, z) ∈ Ω | z > h(x, t ) , (3.1)
Γ h := (x, z) ∈ Ω | z = h(x, t ) . (3.2)
We emphasise that choosing the level h corresponds to the specification of one of the models of our class. The function h can even be an unknown of our problem, more precisely depending of an unknown of the problem (see condition (3.8) below).

On the other hand we introduce the following tensor M 0 which will act as an effective permeability tensor:

M 0 = S 0 0 0 0 , S 0 = K xx - 1 K zz K xz K zx . (3.3)
The 2 × 2 matrix S 0 is the Schur complement of the block K zz in the tensor K 0 . Since K 0 is a symetric positive definite matrix (see just before (2.6)), the same holds for S 0 . We then introduce the averaged conductivity tensor K defined in ]0, T [×Ω x for any function H = H (t , x) by

K ( H )(t , x) = h soil (x) h bot (x) k r ρ g ( H (t , x) -z) M 0 (x, z) d z. (3.4) 5
Finally, for the 2D part of the model, we introduce the notations

∇ x = (∂ x 1 , ∂ x 2 , 0) T for the horizontal gradient and div x (v) = ∇ x • v = ∂ x 1 v 1 + ∂ x 2 v 2 for the horizontal divergence of v ∈ R 3 .
The model. Our coupled model consists in finding the pressure P , the velocity v and the auxiliary unknowns u, w, H and h such that:

• In Ω + h (t ) the following 1D-Richards equation holds              φ ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) α P + β u • e 3 = F for (t , x) ∈]0, T [×Γ soil P t , x, h(t , x) = ρ g H (t , x) -h(t , x) for (t , x) ∈]0, T [×Ω x P (0, x, z) = P init (x, z) for (x, z) ∈ Ω + h (0) (3.5)
• In Ω - h (t ) the pressure P satisfies

P (t , x, z) = ρ g H (t , x) -z for t ∈ [0, T [ , (x, z) ∈ Ω - h (t ) (3.6) • The hydraulic head solves in Ω x        div x K ( H ) ∇ x H = (u • e 3 ) Γ + h for (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 for (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) for x ∈ Ω x (3.7)
where (u • e 3 ) Γ + h denotes the trace of u • e 3 on Γ h from above.

• The level z = h below which we consider the vertical flow to be instantaneous is set such that

h bot (x) ≤ h(t , x) ≤ max min H (t , x) - P s ρ g , h max (x) , h bot (x) , (t , x) ∈ [0, T [×Ω x . (3.8)
• The velocity v is defined in Ω by

         v = u + w for t ∈]0, T [ , (x, z) ∈ Ω u = -k r (P ) 1 ρ g ∂P ∂z + 1 K 0 e 3 for t ∈]0, T [ , (x, z) ∈ Ω w = -k r ρ g ( H -z) M 0 ∇ x H for t ∈]0, T [ , (x, z) ∈ Ω (3.9)
The coupled model (3.5)-(3.9) depends on the definition of the function h. Although all intermediate choices respecting (3.8) are allowed, we will focus in the next on the two extremal choices

h(t , x) = h bot (x), (3.10) h(t , x) = max min H (t , x) - P s ρ g , h max (x) , h bot (x) := h s (t , x), ( 3.11)
and on the intermediate one

h(t , x) = max min H (t , x) - P s + R ρ g , h max (x) , h bot (x) , (3.12)
where R is some positive function possibly depending on H .

The class of models (3.5)-(3.9) is an alternative to the 3D-Richards problem for describing the flow in a shallow aquifer in a large range of time scales. This model is designed to fulfill the two following properties:

• to be simpler to handle numerically than the 3D-Richards model • to behave like the 3D-Richards model for any time scale when the ratio ε of the deepness over the horizontal length of the aquifer is small1 . For example the behaviors presented in Remark 1 are respected.

The first property holds for (3.5)-(3.9) since the 3D original Richards problem is replaced by the coupling of a 2D-problem with a lot of independent 1D-problems which can be solved in parallel. Significant time savings are expected in the computations. The second property is justified in Section 4. The idea is to study the limit ε → 0 of the solution of the 3D-Richards equations and to derive formally the associated effective problem. The same asymptotic analysis is performed for the coupled models (3.5)-(3.9) and shows that the corresponding effective problems are exactly the same for every considered time scale and for every choice of h satisfying (3.8).

Remark 2. It is natural to think that it is possibly not so useful to couple two phenomena which does not hold at the same time scale, since by essence they can not interact with each other. But the notion of time scale is senseless for a fixed physical situation and we just employ this term to enlighten the interpretations. The notion of scale has a precise sense when a sequence of problems is considered, for example parametrised by a small parameter ε tending to zero with the reference time of study depending on ε. This is what we do in Section 4 where ε is the ratio deepness/length of the aquifer. This limit process shows that the two kinds of flow appear at different time scales and then do not interact with each other. Nevertheless, the coupled problem (3.5)-(3.9) is not an effective problem and holds without time scale separation assumption. The depth / width ratio of the aquifer is then a fixed positive number given by the geometry of the aquifer. In particular, "short" and "long" time scales flows can interact without either being negligible or instantaneous.

The remainder of this subsection is devoted to comments on the new models (3.5)-(3.9). Before splitting those comments according to the choice of the function h, we prove that the model is always mass conservative.

Mass conservation.

Let M tot (t ) the total mass of the water contained in domain Ω at time t . We denote by M + h (resp. M - h ) the mass of the water filling the domain Ω + h (resp. Ω - h ). We have

M + h (t ) = ρ Ω x h soil h(t ,x) φ s(P ) d z d x, M - h (t ) = ρ Ω x h(t ,x) h bot (x) φ d z d x, (3.13) M tot (t ) = M + h (t ) + M - h (t ). (3.14)
Proposition 3.1. The total mass satisfies for all t ∈ (0, T ):

∂ ∂t M tot = -ρ Ω x (u • e 3 )| Γ soil d x.
PROOF. By using relation (3.13) and (3.14) it comes

∂ ∂t M tot = ρ Ω x h(t ,x) h bot (x) φ ∂s(P ) ∂t d z d x + ρ Ω x h soil (x) h(t ,x) φ ∂s(P ) ∂t d z d x = ρ Ω x h soil (x) h(t ,x) φ ∂s(P ) ∂t d z d x, ( 3.15)
where the first equality is due to s(P ) = 1 in ]h bot (x), h(t , x)] (indeed P ≥ P s by (3.6) and (3.8)). Thanks to the first equation in (3.5) we deduce

Ω x h soil (x) h(t ,x) φ ∂s(P ) ∂t d z d x = Ω x (u • e 3 )| Γ + h d x - Ω x (u • e 3 )| Γ soil d x. (3.16)
Finally by (3.7) and after an integration by parts

Ω x (u • e 3 )| Γ + h d x = ∂Ω x K ( H )∇ x H • n = 0. (3.17)
The result is obtained by plugging (3.16) and (3.17) in (3.15).

Comments on the model in the case (3.10)

In this case, we have h = h bot , then

Ω + h = Ω, Ω - h = and Γ h = Γ bot (see (3.1)
). The coupled model (3.5)-(3.9) reduces in: finding the pressure P , the velocity v and the auxiliary unknowns u, w and H such that: Fast component of the flow: globally vertical. The unknown u represents the velocity associated with the pressure P by the one dimensional Darcy's law given in the second equation of (3.18). This one is deduced from the 3D law (see the second equation of (2.7)) by neglecting the horizontal components of the gradient of the pressure P . By construction the field u is vertical if the conductivity tensor K 0 introduced in (2.6) is such that K xz = 0 but it may admit a non-zero horizontal component in the anisotropic case. Furthermore the mass conservation equation (3.19) holds. The pressure P then satisfies the following vertical Richards equation where the horizontal variable x ∈ Ω x appears only as a parameter:

         v = u + w for t ∈]0, T [ , (x, z) ∈ Ω u = -k r (P ) 1 ρ g ∂P ∂z + 1 K 0 e 3 for t ∈]0, T [ , (x, z) ∈ Ω w = -k r ρ g ( H -z) M 0 ∇ x H for t ∈]0, T [ , (x, z) ∈ Ω (3.18)              φ ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω α P + β u • e 3 = F for (t , x, z) ∈]0, T [×Γ soil P = ρ g H -h bot for (t , x, z) ∈]0, T [×Γ bot P (0, x, z) = P init (x, z) for (x, z) ∈ Ω (3.19)      -div x K ( H ) ∇ x H = -(u • e 3 )| Γ bot for (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 for (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) for x ∈ Ω x (3.
φ ∂s(P ) ∂t - ∂ ∂z k r (P ) K zz 1 ρ g ∂P ∂z + 1 = 0 in ]0, T [×Ω. (3.21)
The original 3D-Richards problem reduces to the latter equation when the horizontal diffusion terms are neglected. In the short-time scale indeed, those turn to be non-dominant in shallow aquifers as announced in Remark 1 and shown in Section 4.

The boundary condition on Γ soil remains the same than in the 3D-Richards problem. But on the bottom Γ bot , the structure of the boundary condition changes and becomes of Dirichlet type, namely P t , x, h bot (t , x) = ρ g H (t , x) -h bot (t , x) . In fact, even if this Dirichlet condition holds, we do not allow the water flowing out the aquifer through the bottom boundary. Indeed the possibly non-zero flux (u • e 3 )| Γ bot appears as a source term in the first equation of (3.20), so that, as proved in Proposition 3.14, the coupled model is globally mass-conservative. The particular value P = ρ g ( H -h bot ) for the bottom Dirichlet condition, has been chosen so that the fast and slow flows are correctly coupled. This point is further explained in the next paragraph.

8

Slow component of the flow: globally horizontal. On the one hand, introduce the auxiliary pressure Q,

Q := ρ g ( H -z),
for which H plays the role of the hydraulic head. Since H does not depend on z, we have (ρ g ) -1 ∂ z Q + 1 = 0. The first consequence is that the unknown w satisfies (see (3.18))

w = -k r (Q) M 0 ∇ x H .
We recover here the velocity associated to Q by the classical Darcy's law for the conductivity k r (Q) M 0 . The second consequence is that Q is ruled by

∂ ∂z k r (Q) K zz 1 ρ g ∂Q ∂z + 1 = 0 in ]0, T [×Ω,
that is the stationary version of equation (3.21).

On the other hand, we expect P to solve the same stationary problem when the duration of the experiment and when the boundary conditions allow the 1D-Richards problem (3.21) to reach its stationary state. Notice that such a vertical affine profile is also expected in the 3D-Richards model in any non-short time scale (see Remark 1 and Section 4). When this situation occurs, the hydraulic head H := P /ρg + z is constant with respect to z. The Dirichlet boundary condition on h bot in (3.19) then implies that

H (t , x, z) = H t , x, h bot (x) = P t , x, h bot (x) ρg + h bot (x) = H (t , x).
Accordingly, in any non-short time scale, we get H H and then P Q in Ω. This is the reason of the particular choice P = ρ g ( H -h bot ) for the Dirichlet boundary condition on h bot in (3.19). Roughly speaking, the couple (Q, w) characterizes the flow in a long-time experiment in which the vertical flow seems instantaneous with respect to the horizontal one. Unlike the velocity u, the field w is horizontal both in the isotropic and anisotropic cases due to the definition of the tensor M 0 . The computations leading to the definition of M 0 are done in Section 4. Let us give here some qualitative arguments. For large times, w is the main order term of the flow which turns out to be horizontal. The velocity w is also related to some hydraulic head, say L, by the classical Darcy's law w = -k r K 0 ∇L (as in the Richards equation (2.7); see (4.45)). But since w is horizontal we have

0 = w • e 3 = -k r K 0 ∇L • e 3 = -k r K zx ∇ x L -k r K zz ∂L ∂z and then ∂L ∂z = -k r K zx K zz ∇ x L
if K zz = 0 as assumed in this paper, otherwise the question is trivial. Accordingly, in the expression of w = -k r K 0 ∇L, only the term ∇ x L appears and it follows w = -k r M 0 ∇ x L. Notice that the tensor M 0 reduces to K xx in the isotropic case K xz = K zx = 0. Moreover w depends on z only through the term k r (ρ g ( Hz))M 0 which decreases to 0 when z increases above H -P s /ρg . This decrease is fast in general depending on the soil characteristic function k r . Then, roughly speaking, the horizontal component of the flow is maximum in the saturated part and almost vanishing in the unsaturated one far from the capillary fringe.

The evolution of the "stationary pressure" Q is ruled by the first equation of (3.20). This is an horizontal massconservation equation associated with the average velocity w :

= -K ( H )∇ x H = h soil h bot w d z.
The right-hand side is the source term computed from the 1D-Richards problem and which transfers the mass from the vertical description to the horizontal one.

Notice that in this model (3.18)-(3.20), the Dupuit hypothesis is not considered. We precise this point in the next Subsection.

h(t , x) = max min H (t , x) - P s + R ρ g , h max (x) , h bot (x) , (3.22)
for a given non-negative function R possibly depending on H . In the numerical simulations at the end of this section, we consider the constant cases R = 0, corresponding to (3.11), and R = 3. Choosing (3.11) could be guessed as the most intuitive choice since it means in general splitting the domain along the water table, thus separating the flows in the saturated and in the unsaturated areas. But simulations show that it is not necessary the optimal choice for the quality of the 3D-Richards approximation.

Velocity of the flow. As previously, the velocity v of the flow results from the contribution of a fast component u and of a slow one w. The set Ω - h is no more empty in general and an additional brick is introduced in the model for describing the flow in this area. We start by giving some properties of the interface Γ h .

Interface discriminating the flow behaviors. As seen in (3.1), the sets Ω - h (t ) and Ω + h (t ) are characterised by h. In view of the constraint (3.8), the condition

h bot (x) ≤ h(t , x) ≤ h max (x) (3.23)
holds for all (t , x) ∈]0, T [×Ω x . Due to (3.6) and (3.8) the pressure at the level z = h(t , x) satisfies for all (t , x) ∈ ]0, T [×Ω x :

P t , x, h(t , x)      = P s + R if h bot (x) < h(t , x) < h max (x), ≥ P s + R if h(t , x) = h max (x), ≤ P s + R if h(t , x) = h bot (x). (3.24)
In particular, thanks to (2.5) and since R ≥ 0 we get

s P (t , x, z) = 1 if h bot (x) < z ≤ h(t , x), ( 3.25)
which means that the set Ω - h (t ) contains a saturated part of the aquifer for any choice of R ≥ 0. More precisely, the soil is fully saturated in Ω - h (t ) for every t ∈]0, T [ if R > 0, and if R = 0, that is for (3.11), Ω - h can be interpreted as the water table (see Remark 3 below for precisions).

By construction h(t , x) ≤ h max so that the interval ]h(t , x), h soil (x)[ remains non-empty for all (t , x) ∈]0, T [×Ω x . Then we do not have to explicit a direct coupling of the flow in Ω - h with the one in the overland. The coupling between Ω - h and Ω + h is sufficient.

Fast component of the flow: globally vertical, a part being instantaneous. We start by remarking that, as in the previous case, the velocity u is related to P by the vertical Darcy's law (3.9). Moreover the same 1D-Richards equation (3.5) holds, but now, only in the upper part of the aquifer. In particular, in the short-time scale, the dominant vertical flow in Ω + h (t ) remains well described. The main difference between cases h = h bot and h = h bot is related to the vertical flow in the saturated area Ω - h (t ). Indeed, the pressure profile (3.6) now holds in Ω - h and in particular u is zero in Ω - h . As said before, this affine profile is expected in the non-short time scale when the vertical flow appears instantaneous. Hence, the model (3.5)-(3.9) describes precisely the vertical flow in Ω + h and assumes that this flow is instantaneous in Ω - h . Such an assumption is classical in models of saturated shallow aquifers and is known as the Dupuit hypothesis. Then, the model (3.5)-(3.9) in the cases (3.12) can be seen as the coupling of a Dupuit horizontal flow in a saturated part at the bottom of the aquifer with many vertical 1D-Richards flows for a precise description of the leaking fluxes from the overland to the water table.

Notice that, even if h = h bot , the model (3.5)-(3.9) does approximate the 3D-Richards problem at every time scale when the ration ε =deepness / horizontal length tends to zero. Indeed, Proposition 4.1 below holds for any choice of function h such that (3.8) is satisfied. This is explained by the following points in short times: 10 • From the 3D-Richards problem, we expect a vertical description given by the 1D-Richards in the whole Ω, with a vanishing flux at the bottom of the domain (see (4.21)).

• From our model, we get 1D-Richards only in Ω + h with a zero flux in Ω - h (see proof of the short-time scale near equation (4.57)) and the continuity of the pressure.

In fact, these problems are exactly the same.

The field u is non-singular thanks to the continuity condition satisfied by P on Γ h (see (3.5) and (3.6)). As for h = h bot , the particular value of the Dirichlet condition on Γ h has been chosen for a proper coupling of the fast and slow components of the flow. This is further developped in the next paragraph. However if u •e 3 has a trace on the boundary Γ h of Ω + h , this one is non-zero in general whereas u • e 3 = 0 in Ω - h . This is a notable difference with the case h = h bot .

Slow component of the flow. Again, we introduce the auxiliary pressure Q = ρ g ( Hz) and we remark that now P = Q in Ω - h (t ) (even for short times). The fact that P Q in the whole Ω for any non-short times comes, as in the case h = h bot , from the Dirichlet condition P = ρ g ( Hz) which holds on Γ h .

The evolution of (Q, w) is characterized by the evolution of H given in (3.7). In this case where Ω - h (t ) in nonempty in general, we can explicit a little more the dynamic of H . This is detailed in the next paragraph.

Evolution of the hydraulic head. Rewrite the problem (3.7) using the first equation of (3.5) averaged on [h, h soil ]:

-div x K ( H ) ∇ x H = -u Γ soil • e 3 - h soil (x) h(t ,x) φ ∂s(P ) ∂t d z in ]0, T [×Ω x . (3.26)
Since

s(P ) = 1 for z ∈ [h bot , h], we get -div x K ( H ) ∇ x H = -u Γ soil • e 3 - ∂ ∂t h soil (x) h bot (x) φ s(P ) d z in ]0, T [×Ω x , (3.27)
or equivalently by using the Leibniz rule in (3.26) and s(P ) |z=h = 1:

φ| Γ h ∂h ∂t -div x K ( H ) ∇ x H = -u Γ soil • e 3 - ∂ ∂t h soil (x) h(t ,x) φ s(P ) d z in ]0, T [×Ω x . (3.28)
The hydraulic head H is characterized by the latter equation completed by the limit conditions in (3.8). This problem is a non-linear degenerate diffusion equation. Indeed, the diffusion tensor K ( H ) vanishes when H tends to -∞. If moreover (3.11) holds, in view of (3.8), the time derivative can be expressed as

∂h ∂t = C ( H ) ∂ H ∂t with C ( H ) = 1 if H -P s /ρ g ∈]h bot , h max [ 0 if not.
The right-hand side of the first equation in (3.7) plays the role of a source term and represents for each x ∈ Ω the evolution of the amount of water which flows in or out the column ]h(t , x), h soil (x)[ through its lower boundary h(t , x). As we have shown in Proposition 3.1 above, this source term ensures the mass conservation in the coupled model (3.5)-(3.9). Of course this term also depends (non linearly) on the solution H . However this dependence is more easy to handle than the one given in the first equation of (3.7). In particular, the expression (3.28) is well adapted to the numerical implementation of the coupled problem (3.5)- (3.9).

Notice that the level z = h s , defined in (3.11), represents the interface between the saturated and unsaturated part of the aquifer according to the auxiliary pressure Q := ρ g ( H (t , x) -z). In particular Q(t , x, h s (t , x)) = P s if h s (t , x) ∈ (h bot (x), h soil (x)) (regardless of the choice of R ≥ 0 in (3.22)). The conductivity tensor K ( H ) defined in (3.4) can be then decomposed into two parts:

K ( H )(t , x) = C0 + h soil (x) h s (t ,x) k r (Q) M 0 (x, z) d z (3.29)
11
where C0 is the averaged conductivity of the saturated soil, i.e.

C0 = h s (t ,x) h bot (x) M 0 (x, z) d z.
In classical models for the saturated part of an aquifer obtained by vertical integration under the Dupuit's assumption, the definition of the effective conductivity (see for example [6]) reduces to C0 instead of K ( H ), the latter being a little greater. The quantity C0 takes into account the horizontal flow in the saturated part but it ignores the (little) one in the unsaturated part, in particular close to the interface z = h s where the capillary effects lead to a nonnegligible saturation. In practice, the smaller h s , the more significant is the difference K ( H ) -C0 . In particular, if a part of the bottom of the aquifer is not saturated, that is h s = h bot , considering only the vanishing conductivity C0 whereas K ( H ) remains positive is physically incorrect.

Numerical simulations

In this section we compare numerically the original 3D-Richards model (2.7) and the coupled model (3.5)-(3.9) for several choices of h satisfying (3.8).

Physical parameters and geometry. All the simulations are done with the following set of data. Denoting I 3 the 3 × 3 identity matrix we set:

s(P ) = (P s /P ) λ , k r (P ) = (P s /P ) 2+3λ , (P s , λ) = (-1.5, 3), ρ = 1, φ = 0.1, K 0 = 0.1 I 3 .
To lighten the numerical results, we consider the simplified 2D aquifer Ω =] -5, 0[×Ω x , Ω x =]0, L x [. In the experiments illustrated in Figures 3 and4, the horizontal length is L x = 28. In those of Figure 5, L x ∈ [21,393]. The parameter δ in (2.3) is chosen as small as possible, that is equal to the size of one vertical mesh. We assume an impermeable layer at the bottom and the top of the aquifer.

Visualisation. For the visualization of the results, we introduce a function h sat representing in a lot of cases the top level of the saturated region at the bottom of the aquifer (i.e. the water table). Let h sat = h sat (t , x) and the set Ω - h sat (t ) be defined for a given pressure P = P (t , x, z) by

h sat (t , x) := sup I t ,x , I t ,x := z ∈ [h bot (x), h max (x)] | P (t , x, z ) > P s , ∀z ∈ [h bot (x), z[ , (3.30
) does not coincide with the saturated region of the soil at the bottom of the aquifer. Indeed a saturated region just over z = h sat is possible for example if P ≥ P s also in Ω \ Ω - h sat . The interface z = h sat then describes • either the interface between the saturated part at the bottom of the aquifer and the unsaturated part above in the simplest setting,

Ω - h sat (t ) := (x, z) ∈ Ω | z < h sat (t , x) . ( 3
• or a level between two saturated part when for example a saturated front flow down and reach Ω -

h sat ,
• or the bottom of the aquifer when h sat = h bot , that is when there is no saturated part at the bottom,

• or the maximum allowed height h sat = h max when, roughly speaking, the water table overflows.

Of course here, since h sat (t , x) ≤ h soilδ by (2.3), the set Ω - • In the case (3.12), the nonlinear coupling between the 1D-vertical Richards equations and the 1D-horizontal diffusion equation is performed by using a Picard's fixed-point method at each time step. This one alternatively solves (3.5) (for an explicit H and h) and (3.28) (for an explicit right-hand side).

In any case all the 1D-Richards equations remain independent at the discrete level and can be solved in parallel.

Reference flowing experiment. At time t = 0, we consider a setting where the function h sat introduced in (3.30) corresponds to the height of the water table. To show the influence of the deepness of the saturated area, we choose a function h sat (0, •) which goes smoothly from -4.5 on the left part of the aquifer to -2.5 on the right one:

h sat (0, x) = -4.5 + 2 e -( 15
Lx ) 2 (x-0.55 L x ) 2 in [0, 0.55L x ], -2.5 in ]0.55L x , L x ].
The initial pressure P is defined by P (0, x, z) = ρ g (h sat (0, x) -z) + P s for all (x, z) except near two rectangular regions above z = h sat where the pressure goes smoothly to the saturation value P s , corresponding to an infiltration process. These rectangles are

R 1 =]L x /10, 3L x /10[× ]-3.5, -1.7[ and R 2 =]7L x /10, 9L x /10[× ]-2, -0.2[ . (3.32)
This initial situation is drawn in the first picture of Figure 3. In every picture the gray scale corresponds to the saturation value, the maximal darkness corresponding to s 1. The total time of the experiment is 4 days. The solution of the classical Richards problem at time 0, 10, 20 and 96 hours respectively, is drawn in Figure 3. The graph of the visualization function h sat defined in (3.30) is also plotted. Its evolution will be used for comparing the original Richards model with the coupled model (3.5)-(3.9).

At time t = 10 the water initially in rectangles R 1 and R 2 started to flow down. In the right part, some water coming from R 2 have reached the saturated water table inducing an increase of its level. In the mean time, we see in the middle of the domain Ω x that the water moves to the left and that the function h sat is smoother than the initial one.

At time t = 20 the water initially in rectangle R 1 has continued to flow down and is about to reach the water table. It is important to notice that this flow was essentially along the vertical direction. In particular the water front which is very close to h sat is approximately horizontal as in the initial situation.

After some time almost all the water initially located in the rectangle supplies have reached the water table. Then the interface h sat becomes flat and is associated with a pressure admitting the stationary profile P (t , x, z) = P s + ρ g h sat (t , x) -z .

Comparison of the models. In this part we compare the solution of the classical Richards model with the one obtained by using the coupled model (3.5)-(3.9). We test three particular choices for the function h satisfying (3.8): the minimal one (3.10), the maximal one (3.11) and an intermediate one given by (3.12) for R = 3. All data remain the same as in the previous paragraph. In this paper, we focus on the evolution of the functions h sat defined by (3.30). As indicated in Remark 3, this function roughly represents the upper level of the water table. In the following we denote by h 2d sat the level coming from the reference 2d-Richards model and we denote by The functions h 2d sat , h a sat , h b sat and h c sat are plotted in Figure 4 at time t ∈ {10, 24, 48, 96} (in hours). We of course do not plot the initial situation which is the same for each model and is the one of the reference test case described in the previous paragraph. The curve h 2d sat is the reference one and is plotted with a black solid line in Figure 4.

13

Impermeable rock

h sat h soil , h bot t = 0 0 x 28 0 z -5 R 1 R 2 t =

15

Bear in mind that the function h characterizes the level below which the vertical flow is assumed to be instantaneous (instead of being described by the 1D-Richards equation). In every case, the horizontal flow is ruled by equation (3.7).

• In the case (3.10), h = h bot . The vertical flow is described by the 1D-Richards model in the whole domain, even in the saturated part below the level z = h a sat . The horizontal flow in this case seems to be slower than the one given by the Richards model (compare the gray dot-dashed line with the black solid one in Figure 4). Roughly the idea is that in this case the water have to travel along the whole vertical direction before reaching the level z = h = h bot . Then the flux (u • e 3 )| Γ bot at the bottom of the aquifer takes a lot of time to increase when the water coming from rectangles R 1 and R 2 reaches the water table. This flux being the source term in equation (3.7), the function H increases with some delay and the corresponding horizontal flow is slower.

• In the case (3.11), h = h b

sat . This case is opposite of the previous one in the sense that the vertical flow in the whole saturated zone Ω - h sat is considered to be instantaneous. Then, when the water coming from rectangles R 1 and R 2 reaches the water table, the flux (u • e 3 )| Γ h increases very quickly. So does the corresponding hydraulic head H and the horizontal flow is very and even too fast (see the black dotted line compared to the black solid line in Figure 4).

• In the case (3.12) for R = 3, h bot ≤ h ≤ h c sat . The corresponding flow should exhibit an intermediate behavior between the too previous ones. Here, the value R = 3 was chosen so that h c sat is very close to the reference one h 2d sat (see the gray dashed line).

Notice that in every situation, the error between h 2d sat and h κ sat , κ ∈ {a, b, c}, is smaller in the left part of the domain than in the right one. This is due to the fact that the saturated zone is thiner in this region. For a very thin saturated region, considering an instantaneous vertical flow or the one given by the vertical 1D-Richards problem gives similar results. Conversely, the thicker the saturated water table is, the more the results issued from the two extremal situations (3.10) and (3.11) differ from the reference one. Basically, h b sat is expected to move too fast while h a sat moves too slowly. In this kind of deep situation and if the ratio between the deepness and the length of the aquifer is not so small, one of the intermediate choices (3.12) is obviously more appropriate.

Error made by the coupled model versus the ratio deepness/largeness. In the previous simulations, where Ω = ]0, 28[×] -5, 0[, the ratio ε=deepness/length of the aquifer is such that 1/ε = 5.6. It is important to notice that even in this case of large ratio ε the error between the original Richards model and the coupled model (3.5)-(3.9) in the case (3.12) is particularly small (see the dashed plot in Figure 4). This supports the fact that the coupled model may be considered for approaching the Richards model also in an aquifer which is not so shallow. This guess is confirmed by the results plotted in Figure 5. The evolution of the error h 2d sat -h κ sat L 1 (]0,T [×Ω x ) for κ ∈ {a, b, c} is drawn in terms of the ratio 1/ε.

As expected all the errors decrease with ε. Moreover, the intermediate case (3.12) is always the best, mainly in the case of a "large" value of ε. After comes the maximal choice. The worst choice is the maximal one (3.10) but with an error which decreases a lot with ε.

Remark 4.

The accuracy of the model depends on the choice of R in (3.22), e.g. for minimizing the error h 2d sat -

h κ sat L p (]0,T [×Ω x )
. This optimization process is postponed to a forthcoming work.

FORMAL ASYMPTOTIC EXPANSION

In this section, the 3D-Richards problem (2.7) and the coupled model (3.5)-(3.9) are compared using asymptotic analysis arguments. We prove that these models behave the same, whatever the time scale, when the ratio between the characteristic deepness and the length of the shallow aquifer tends to zero.

c}). Function h 2d

sat is the iso-pressure P = P s in the original 2d-Richards problem and h κ sat is the one associated with the coupled problems for three different choices of h satisfying (3.8). The test case is the one of figure 3.

Dimensionless form of the 3D-Richards and coupled problems

Introduce a fixed dimensionless reference domain Ω of type (2.2) and a dimensionless real number T > 0. Fix Ω x , h soil and h bot such that

Ω = (x, z) ∈ Ω x × R | z ∈ h bot (x), h soil (x) .
To obtain a rescaled version of equations (2.7) and (3.5)-(3.9) in the domain ]0, T [×Ω, we introduce positive reference numbers L x , L z , T . Then, keeping the same notations as in Section 3, we have:

• The physical variables are given by

x = L x x, z = L z z, t = T T t .
• The corresponding physical domain Ω is given as in (2.2) with

Ω x = L x Ω x , h soil (x) = L z h soil (x), h bot (x) = L z h bot (x).
• The unknowns are such that

P (t , x, z) = P (t , x, z), v(t , x, z) = v(t , x, z), u(t , x, z) = u(t , x, z), w(t , x, z) = w(t , x, z), L z H (t , x) = H (t , x), L z h(t , x) = h(t , x).
• The reference subdomains are

Ω - h (t ) = (x, z) ∈ Ω x × R | z ∈ h bot (x), h(t , x) , Ω + h (t ) = (x, z) ∈ Ω x × R | z ∈ h(t , x), h soil (x) .
• The reference boundaries are

Γ bot := {(x, z) ∈ Ω | z = h bot (x)}, Γ soil := {(x, z) ∈ Ω | z = h soil (x)} and Γ ver := {(x, z) ∈ Ω | x ∈ ∂Ω x }.
• The reference exterior normals are

n(x, z) =                e 3 - L z L x ∇ x h soil (x) L 2 z L 2 x |∇ x h soil (x)| 2 + 1 -1/2 on Γ soil L z L x ∇ x h bot (x) -e 3 L 2 z L 2 x |∇ x h bot (x)| 2 + 1 -1/2 on Γ bot n(x, z) on Γ ver
where the vector n is horizontal and does not change during the rescaling. 17 It means that the reference saturation and relative permeability are of order one. Indeed P and P take the same values, independently of the scale change.

• For the conductivities, we set

K 0 (x, z) = K 0 (x, z), M 0 (x, z) = M 0 (x, z), (4.2) K (H )(t , x) = L z h soil (x) h bot (x) k r ρ g (H (t , x) -z M 0 d z. (4.3)
We choose (4.2) for the sake of simplicity in the presentation. Indeed, we could also introduce K and M such that K K 0 (x, z) = K 0 (x, z) and M M 0 (x, z) = M 0 (x, z) and then perform the same study assuming that

K /L x = O (ε), M /L x = O (ε) and K /L z = O (1).
• The source term is

F (t , x) = F (t , x)
Dimensionless Richards problem. Introducing the latter quantities in (2.7), we get the following set of rescaled equations:

T T φ ∂s(P ) ∂t + 1 L x div x (v) + 1 L z ∂v ∂z = 0 in ]0, T [×Ω, (4.4
) v = -k r (P ) K 0 1 L x 1 ρg ∇ x P + 1 L z 1 ρg ∂P ∂z + 1 e 3 in ]0, T [×Ω, (4.5)
v • L z L x ∇ x h bot -e 3 = 0 on ]0, T [×Γ bot , (4.6
)

α P L 2 z L 2 x ∇ x h soil 2 + 1 1/2 + β v • e 3 - L z L x ∇ x h soil = F L 2 z L 2 x ∇ x h soil 2 + 1 1/2 on ]0, T [×Γ soil , (4.7) v • n = 0 on ]0, T [×Γ ver . (4.8)
Since the aquifer is assumed to be very thin with respect to its horizontal width, the quantity L z /L x is very small. We choose to consider an aquifer with a fixed height of order L z = 1 and a large horizontal dimension L x = 1/ε for ε 1. We get • the mass conservation equation which depends on the time scaling choice T :

T T φ ∂s(P ) ∂t + ε div x (v) + ∂v • e 3 ∂z = 0 in ]0, T [×Ω (4.9)
• associated with the following Darcy's law and boundary conditions: • The velocity problem:

                   v = -k r (P ) K 0 ε ρg ∇ x P + 1 ρg ∂P ∂z + 1 e 3 in ]0, T [×Ω α P ε 2 ∇ x h soil 2 + 1 1/2 + β v • e 3 -ε ∇ x h soil = ε 2 ∇ x h soil 2 + 1 1/2 F on ]0, T [×Γ soil v • n = 0 on ]0, T [×Γ ver v • ε∇ x h bot -e 3 = 0 on ]0, T [×Γ bot
         v = u + w for t ∈]0, T [ , (x, z) ∈ Ω u = -k r (P ) 1 ρ g ∂P ∂z + 1 K 0 e 3 for t ∈]0, T [ , (x, z) ∈ Ω w = -ε k r ρ g (H -z) M 0 ∇ x H for t ∈]0, T [ , (x, z) ∈ Ω (4.11)
• The 1D-Richards equation in the transition zone:

               φ T T ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) α P + β u • e 3 = F for (t , x) ∈]0, T [×Γ soil P t , x, h(t , x) = ρ g H (t , x) -h(t , x) for (t , x) ∈]0, T [×Ω x P (0, x, z) = P init (x, z) for (x, z) ∈ Ω + h (0) (4.12)
• The pressure problem in the water table:

P (t , x, z) = ρ g H (t , x) -z for t ∈ [0, T [ , (x, z) ∈ Ω - h (t ) (4.13)
• The hydraulic head problem:

       ε 2 div x K (H ) ∇ x H = u Γ + h • e 3 for (t , x) ∈]0, T [×Ω x K (H )∇ x H • n = 0 for (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) for x ∈ Ω x (4.14)
Equivalently, by using (3.27), the first equation of (4.14) admits the formulation: for (t ,

x) ∈]0, T [×Ω x ε 2 div x K (H ) ∇ x H = u Γ soil • e 3 + T T ∂ ∂t h soil (x) h bot (x)
φ s(P ) d z (4.15)

• The definition of the interface separating the two different kind of flows:

h bot (x) ≤ h(t , x) ≤ max min H (t , x) - P s ρ g , h max (x) , h bot (x) for (t , x) ∈ [0, T [×Ω x (4.16)

Effective problems

We are interested in the asymptotic behavior of the flow, thus of the models, for both short, intermediate and large times. For the asymptotic analysis, the question is related to the behavior of the dimensionless models above. More precisely, we want to describe the effective flow obtained for the short time T = T , the intermediate time T = ε -1 T and the long time scales T = ε -2 T . Asymptotic expansion. We introduce the following formal asymptotics for the pressure and the velocity:

P γ ε = P γ 0 + ε P γ 1 + ε 2 P γ 2 + . . . v γ ε = v γ 0 + ε v γ 1 + ε 2 v γ 2 + . . . ( 4

.17)

We emphasize that no arbitrary scaling is imposed, in particular we do not suppose as in [10] that the vertical velocity is much smaller than the horizontal one when the ratio is very small. We assume also the existence of formal asymptotics for the auxiliary variables appearing in (3.5)-(3.9) Effective problems at the main order. Let us introduce the following effective problems:

u γ ε = u γ 0 + ε u γ 1 + ε 2 u γ 2 + . . . w γ ε = w γ 0 + ε w γ 1 + ε 2 w γ 2 + . . . H γ ε = H 0 + ε H γ 1 + ε 2 H γ 2 + . . . h γ ε = h γ 0 + ε h γ 1 + ε 2 h γ 2 + . . . , (4.18
F ε = F 0 + ε F 1 + ε 2 F 2 + . . . . ( 4
• related to the short time scale (T = T ),

                   φ ∂s(P 0 ) ∂t + ∂v 0 • e 3 ∂z = 0 in ]0, T [×Ω v 0 = -k r (P 0 ) 1 ρg ∂P 0 ∂z + 1 K 0 e 3 in ]0, T [×Ω α P 0 + β v 0 • e 3 = F 0 on ]0, T [×Γ soil v 0 • e 3 = 0 on ]0, T [×Γ bot (4.21)
• related to the non-short cases (T = ε -1 T or T = ε -2 T ),

P 0 (t , x, z) = ρ g H 0 (t , x) -z in ]0, T [×Ω v 0 = 0 in ]0, T [×Ω (4.22)
• related to the non-short cases

(T = ε -1 T or T = ε -2 T ) if α = 0 H 0 (t , x) = F 0 (t , x) α ρ g + h soil (t , x) in ]0, T [×Ω x (4.23)
• related to the intermediate time scale (T = ε -1 T ) if α = 0 (and then β = 0)

ρ g h soil h bot φ s (P 0 ) d z ∂H 0 ∂t = - F 1 β in ]0, T [×Ω x (4.24)
• related to the long time scale (4.25) and concerning the first order of the velocity

(T = ε -2 T ) if α = 0      -div x K (H 0 ) ∇ x H 0 = - F 2 β - ∂ ∂t h soil h bot φ s(P 0 ) d z in ]0, T [×Ω x K (H 0 ) ∇ x H 0 • n = 0 on ]0, T [×Γ ver
v 1 = -k r (P 0 ) M 0 ∇ x H 0 in ]0, T [×Ω (4.
= F 1 = 0 if α = 0. Moreover v 2 1 satisfies (4.26) if α = 0.
We emphasize that the intermediate variable h which characterizes the coupled model (4.11)-(4.15) does not appear in any of the main order effective problems (4.21)-(4.25). This agrees with the fact that the whole class of models given by (3.5)-(3.9) for any h satisfying (3.5) can approximate the reference Richards model. 20

Proof of Proposition 4.1 for the Richards model

The proof of Proposition 4.1 consists in substituting the formal asymptotic expansion (4.17)-(4.20) in the rescaled 3D-Richards problem (4.9)-(4.10). A cascade of equations follows by identifying the powers of ε. Then we characterize the main order terms in the expansion (4.17). In order to reduce ratings in this section, we do not write the exponent γ on the variables name.

General relations.

Let us start by obtaining the first relations holding in every time scale (i.e. for all γ ∈ {0, 1, 2}). By plugging the asymptotic expansion (4.17) in the first equation of (4.10) we get the following relations holding in ]0,

T [×Ω          v 0 = -k r (P 0 ) 1 ρ g ∂P 0 ∂z + 1 K 0 e 3 , v 1 = - k r (P 0 ) ρ g K 0 ∇ x P 0 + ∂P 1 ∂z e 3 -k r (P 0 ) P 1 1 ρ g ∂P 0 ∂z + 1 K 0 e 3 .
(4.27)

The same process in the three last equations of (4.10) yields the following relations in ]0, T [:

• on Γ soil    α P 0 + β v 0 • e 3 = F 0 , α P 1 + β v 1 • e 3 -v 0 • ∇ x h soil = F 1 , α P 2 + 1 2 ∇ x h soil 2 P 0 + β v 2 • e 3 -v 1 • ∇ x h soil = 1 2 ∇ x h soil 2 F 0 + F 2 ; (4.28) • on Γ bot , for all k ∈ N * v 0 • e 3 = 0, v k-1 • ∇ x h bot = v k • e 3 ; (4.29)
• on Γ ver , for all k ∈ N v k • n = 0. The latter equation completed with the first equations of (4.27), (4.28) and (4.29) gives exactly the system (4.21). The first claim of Proposition 4.1 is proven.

Intermediate time case.

In this part, we prove the second claim of Proposition 4.1 which is associated with the intermediate time scale T = ε -γ T for γ = 1. Equation of (4.9) is now

εφ ∂s(P ) ∂t + ε div x (v) + ∂v • e 3 ∂z = 0. (4.33)
We introduce the asymptotic expansion (4.17) in the previous equation and we identify the main order terms. We obtain ∂v 0 • e 3 ∂z = 0 on ]0, T [×Ω. (4.34) This constant vertical velocity is actually zero due to (4.29). Moreover, with the first equation of (4.27) and since k r and (K 0 ) 33 are non-vanishing (K 0 is positive definite), we get in ]0, T [×Ω

∂P 0 ∂z + ρ g = 0 and v 0 = 0. (4.35) 21
The existence of H 0 = H 0 (t , x) such that

P 0 (t , x, z) = ρ g H 0 (t , x) -z in ]0, T [×Ω (4.36)
follows. Next, since v 0 = 0, the first equation of (4.28) is

α P 0 = F 0 on Γ soil . (4.37)
We now have to differentiate the computations depending on whether α = 0 or not. If α = 0, then for all (t , x) ∈]0, T [×Ω x we have P 0 (t , x, h soil (t , x)) = F 0 (t , x)/α. Accordingly, thanks to (4.36), it holds

H 0 (t , x) = F 0 (t , x) α ρ g + h soil (t , x).
This ends the proof of the second claim of Proposition 4.1 in the case α = 0. If α = 0 (then β = 0), equation (4.37) only implies that F 0 = 0. In particular, H 0 remains as a degree of freedom and we have to exploit the next order terms in the asymptotic expansion for the closure of the effective problem. Identifying the coefficients associated with ε 1 in equation (4.33) we have φ ∂s(P 0 )

∂t + ∂v 1 • e 3 ∂z = 0 in ]0, T [×Ω. (4.38)
To eliminate v 1 , we integrate vertically on ]h bot , h soil [ the equation above. After using the fact that ∂ t (s(P 0 )) = ρ g s (P 0 ) ∂ t H 0 (consequence of (4.36)) we have

ρ g h soil h bot φ s (P 0 ) d z ∂H 0 ∂t + (v 1 | h soil -v 1 | h bot ) • e 3 = 0. (4.39)
Thanks to the second equations of (4.28) and (4.29) in the case where α = 0 and v 0 = 0, it follows:

v 1 • e 3 = F 1 /
     P 0 (t , x, z) = ρ g H 0 (t , x) -z in ]0, T [×Ω ρ g h soil h bot φ s (P 0 ) d z ∂H 0 ∂t = - F 1 β in ]0, T [×Ω x (4.41)
which correspond to the second claim of Proposition 4.1 in the case α = 0.

Long time case.

In this part, we prove the third claim of Proposition 4.1 which is associated with the intermediate time scale T = ε -γ T for γ = 2. Equation (4.9) takes the form

ε 2 φ ∂s(P ) ∂t + ε div x (v) + ∂v • e 3 ∂z = 0. (4.42)
We substitute the asymptotic expansion(4.17) in the previous equation. The main order part of the equation is ∂ z (v 0 • e 3 ) = 0 which leads, as before, to (4.22) for some function H 0 which does not depends on z. The same relation (4.37) holds and the characterization of H 0 depends on the values of α. As before, if α = 0 we have (4.23).

22

It remains to deal with the case α = 0 and to exhibit the equations of system (4.25). In this case, the compatibility condition F 0 = 0 holds as before because of (4.37). The characterization of H 0 needs to go at the next order in the asymptotic expansion. In equation (4.42) we get

0 = div x (v 0 ) + ∂v 1 • e 3 ∂z = ∂v 1 • e 3 ∂z (4.43)
where the second equality is due to v 0 = 0. Moreover, the second equations of (4.28) and (4.29) for k = 1 lead to (since α = 0)

β v 1 • e 3 = F 1 on Γ soil and v 1 • e 3 = 0 on Γ bot . (4.44)
Then, the vertical component of the velocity (which is constant by (4.43)) v 1 • e 3 is zero. Moreover the second compatibility condition F 1 = 0 holds true thanks to (4.44). Using the second equation of (4.27) and bearing in mind that (ρ g ) -1 ∂ z P 0 + 1 = 0, we obtain

v 1 = - k r (P 0 ) ρ g K 0 ∇ x P 0 + ∂P 1 ∂z e 3 . ( 4

.45)

Since v 1 • e 3 = 0, using the same notation for K 0 than in (2.6), we compute ∂ z P 1 by

∂P 1 ∂z = - 1
K zz K 0 ∇ x P 0 • e 3 .
Finally, substitution in the equation above with the relation

P 0 = ρ g (H 0 -z) give v 1 = -k r (P 0 ) M 0 ∇ x H 0 with M 0 = I 2 -K xz K zz 0 0 K 0 = S 0 0 0 0 (4.46)
where I 2 is the 2d identity matrix and S 0 = K xx -K -1 zz K xz K zx . On the other hand, the equation (4.30) for k = 1 leads to v 1 • n = 0 on Γ ver . Since k r (P 0 ) does not vanish, we obtain the last equation of (4.25). After identifying the coefficients associated with ε 2 in equation ( 4

v 2 • e 3 -v 1 • ∇ x h soil = F 2 /β, v 2 • e 3 -v 1 • ∇ x h bot = 0 on Γ bot . ( 4
h bot div x v 1 d z + v 1 | h soil • ∇ x h soil + F 2 β -v 1 | h bot • ∇ x h bot = 0.
We use the Leibniz rule in the second integral and we get

∂ ∂t h soil h bot φ s(P 0 ) d z + div x h soil h bot v 1 d z = - F 2 β . ( 4

.49)

Using the averaged conductivity K defined in (4.3), we get, with the first equation of (4.46),

h soil h bot v 1 d z = - h soil h bot k r (P 0 ) M 0 ∇ x H 0 = -K (H 0 ) ∇ x H 0 .
The above equation associated with equation (4.49) is exactly the system (4.25). This ends the proof of the last claim of Proposition (4.1). 23

Proof of Proposition 4.1 for the coupled models

The strategy of the proof is exactly the same than in the previous subsection.

General relations. Let γ ∈ {0, 1, 2}. Using the expansion (4.17)-(4.20), we identify powers of ε in all the equations in (4.11)-(4.16) that does not depend on the time scale T . We obtain from the second equation of (4.11)

         u 0 = -k r (P 0 ) 1 ρ g ∂P 0 ∂z + 1 K 0 e 3 in ]0, T [×Ω, u 1 = - k r (P 0 ) ρ g ∂P 1 ∂z K 0 e 3 -k r (P 0 ) P 1 1 ρ g ∂P 0 ∂z + 1 K 0 e 3 in ]0, T [×Ω, (4.50)
from the third equation of (4.11)

w 0 = 0, w 1 = -k r ρ g (H 0 -z) M 0 ∇ x H 0 in ]0, T [×Ω, (4.51)
from the first equation of (4.11)

     v 0 = u 0 + w 0 = u 0 = -k r (P 0 ) 1 ρ g ∂P 0 ∂z + 1 K 0 e 3 in ]0, T [×Ω, v 1 = u 1 + w 1 in ]0, T [×Ω. (4.52)
It follows from (4.13) that for t ∈]0, T [ and (x, z) ∈ Ω - h 0 (t ) For the boundary conditions, we infer from the second and third equations of (4.12) and from the second equation of (4.14) that, for all k ∈ N,

P 0 (t , x, z) = ρ g H 0 (t , x) -z , P k (t , x, z) = ρ g H k (t , x) ∀k > 0. ( 4
     α P k + β u k • e 3 = F k on ]0, T [×Γ soil , P 0 t , x, h 0 (t , x) = ρ g H 0 (t , x) -h 0 (t , x) for t ∈]0, T [, x ∈ Γ h (t ), K (H 0 )∇ x H 0 • n = 0 on ]0, T [×Γ ver .
u 1 = 0 in Ω - h 0 (t ). ( 4

.56)

Short time case. In this part, T = T , that is γ = 0. The first equations of (4.12) and (4.14) become

       φ ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ), ε 2 div x K (H ) ∇H = (u 0 • e 3 )| Γ h for (t , x) ∈]0, T [×Ω x .
(4.57)

We identify the main order terms appearing when the asymptotics (4.17 From (4.52) and (4.53) we also compute u 0 = 0 in Ω - h 0 (t ). In addition, from (4.54) we get s(P 0 ) = 1 in Ω - h (t ) so that (P 0 , u 0 ) satisfies (4.58) also in Ω - h (t ). The continuity of u 0 • e 3 being ensured by (4.59), (P 0 , u 0 ) satisfies (4.58) in the whole Ω. By using (4.52) and (4.55) we obtain the system (4.21) and then the first claim of Proposition 4.1 holds once again.

Intermediate time case.

In this part, T = ε -1 T , γ = 1. The first equation of (4.12) and the equation (4.15) become

           φε ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) -ε 2 div x K (H ) ∇ x H = -(u • e 3 )| Γ soil -ε ∂ ∂t h soil (x) h bot (t ,x) φ s(P ) d z for (t , x) ∈]0, T [×Ω x (4.60)
The corresponding main order relations are It follows that the constant vertical component of the velocity u 0 • e 3 equals zero in Ω + h (t ). We deduce from the first equation of (4.50) that the pressure P 0 is affine with respect to the z variable with the slope -ρ g in Ω + h (t ). Accordingly, thanks to the first equation of (4.53) and the continuity condition given in (4.55), the first equation of (4.22) holds. Using relation (4.52) we obtain the second equation of (4.22). Next, thanks to u 0 = 0 and to the first equation of (4.55) for k = 0, we get α P 0 = F 0 .

u 0 • e 3 = 0 on ]0,
If α = 0 then for all (t , x) ∈]0, T [×Ω x we have P 0 (t , x, h soil (t , x)) = F 0 (t , x)/α. Accordingly, thanks to the first equation of (4.22), we have

H 0 (t , x) = F 0 (t , x) α ρ g + h soil (t , x).
The second claim of Proposition 4.1 in the case α = 0 is proved. If α = 0, the compatibility condition F 0 = 0 is imposed. After identifying the coefficients associated with ε 1 in the second equation of (4.60) we have

0 = -(u 1 • e 3 )| Γ soil - ∂ ∂t h soil (x) h bot (x)
φ s(P 0 ) d z and, with the first equation of (4.22),

ρ g h soil (x) h bot (x) φ s (P 0 ) d z ∂H 0 ∂t = -(u 1 • e 3 )| Γ soil .
The first equation of (4.55) for k = 1 implies, since α = 0, that (u 1 •e 3 )| Γ soil = F 1 /β. This ends the proof of the second claim of Proposition 4.1 in the case α = 0.

Long time case.

In this part, T = ε -γ T , γ = 2. The first equation of (4.12) and equation (4.15) are now It remains to deal with the case α = 0. In this case we first remark that the compatibility condition F 0 = 0 holds (see (4.55) for k = 0). Furthermore, since P 0 = ρ g (H 0 + z) we get from (4.56) and (4.64) that u 1 = 0 in ]0, T [×Ω. Thus, using (4.51) and (4.52) we get v 1 = w 1 = -k r (P 0 )M 0 ∇ x H 0 . Moreover the first equation of (4.55) for k = 1 gives F 1 = 0 (since α = 0). It remains to get the first relation of system (4.25). By plugging asymptotics (4.17)-(4.20) in the second equation of (4.63) and by identifying the coefficients associated with ε 2 we get ((4.50)) The first purpose of this work is to describe a numerical conservative scheme to approximate the non-linear coupled model. The scheme is fully implicit in time. It is based on a reformulation of the problem in which the non-linear coupling is reduced in a Dirichlet to Neumann operator. The resulting equation is roughly speaking a 2d mass conservation equation (non-linear and degenerative) associated to a non-linear Darcy's law. The numerical resolution of this problem is done at each time step thanks to a quasi-Newton method. Two different formulations of the Newton direction are given to deal with general flow situations leading to possibly overflowing or empty aquifers.

         φ ε 2 ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) -ε 2 div x K (H ) ∇ x H = -(u • e 3 )| Γ soil -ε 2 ∂ ∂t h soil (x) h bot (x) φ s(P ) d z for (t , x) ∈]0, T [×Ω x
-div x K (H 0 ) ∇ x H 0 = -(u 2 cd ot e 3 )| Γ soil - ∂ ∂t h soil (x) h bot (x) φ s(P 0 ) d z for (t , x) ∈]0, T [×Ω x . ( 4
The second goal of this paper is to use numerical simulations to compare the coupled model with the original Richards one. We know by [2] that those models behaves the same when the aquifer admits a small ratio ε=deepness/longitudinal length. The purpose here is to quantify the difference beetween them for non-vanishing ε. We show in particular that in most situations, a not so shallow geometry is needed to have a good approximation.

GEOMETRY AND GOVERNING MODEL

Geometry. The aquifer consists in thin and large cylindrical domain Ω ⊂ R 3 in the vertical direction e 3 . The projection of the Ω on an horizontal plane is Ω x ⊂ R 2 and the boundary of Ω is ∂Ω. The bottom of the domain and the upper level of the domain are the graphs of functions h bot = h bot (x) and h soil = h soil (x) respectively. These functions are definied from Ω x to R and satisfy

h soil (x) > h bot (x) ∀x ∈ Ω x . (1.1)
The domain is given by:

Ω = (x, z) ∈ Ω x × R | z ∈]h bot (x), h soil (x)[ . (1.2)
The boundary ∂Ω is splitted in three zones as follows

∂Ω = Γ bot Γ soil Γ ver with Γ bot := (x, z) ∈ Ω | z = h bot (x) , Γ soil := (x, z) ∈ Ω | z = h soil (x) , Γ ver := (x, z) ∈ Ω | z ∈ ∂Ω x . ( 1

.3) NUMERICAL ASPECTS OF A CLASS OF MODELS COUPLING SLOW AND FAST FLOW IN SHALLOW AQUIFERS

2 We introduce for given 0 < δ 1

+ x z z = h bot (x) z = h soil (x) z = h max (x)
h max (x) := h soil (x) -δ. (1.4)
We represent the typycal geometry in Figure 1 in the case of 2d domain where Ω x is an interval.

Unknowns and constitutive parameters. Each model of the class proposed in [2] characterizes the flow in a shallow aquifer. The main unknowns are the water pressure denoted by P = P (t , x, z) and the corresponding velocity v = v(t , x, z). The models will also determine the auxiliary unknowns H , h, u 3 and w corresponding respectively to

• H : the hydraulic head of the water table

• h: the interface between the two different descriptions of the flow • u: the fast component of the velocity • w: the slow component of the velocity

We consider the saturation and the relative permeability of the soil characterized as non decreasing functions of the pressure P . They are denoted respectively by s = s(P ) and k r = k r (P ). Moreover we assume that they are such that s(P ) = 1 ⇐⇒ P ≥ P s and k r (P ) = 1 ⇐⇒ P ≥ P s (1.5) where P s ∈ R is given. For example those conditions are satisfied if the classical Books and Correy model is considered to describe the soil:

s(P ) =      P s P λ if P < P s 1 if P ≥ P s , k r (P ) =    P s P γ if P < P s 1 if P ≥ P s (1.6)
where λ > 0, γ = 2 + 3λ and P s < 0. In every numerical simulations of section 3 this model will be consider.

The soil porosity is denoted by φ = φ(x, z) ∈]0, 1[. The permeability of the saturated soil is a 3 × 3 symmetric positive definite tensor K 0 = K 0 (x, z). We introduce

K xx ∈ M 22 (R), K zz ∈ R * and K xz ∈ M 21 (R) such that K 0 = K xx K xz K T xz K zz . (1.7)
On the other hand we introduce the following tensor M 0 which will be used to characterize an effective permeability tensor in the horizontal direction:

M 0 = S 0 0 0 0 , S 0 = K xx - 1 K zz K xz K zx . (1.8)
The 2 × 2 matrix S 0 is the Schur complement of the block K zz in the tensor K 0 . Since K 0 is a symmetric positive definite matrix, the same holds for S 0 . We then introduce the averaged permeability tensor K defined in ]0, T [×Ω x for any function H = H (t , x) by

K (H )(t , x) = h soil (x) h bot (x) k r ρ g (H (t , x) -z) M 0 (x, z) d z.
(1.9)

We assume that there is only one incompressible fluid which flow into the aquifer. Its corresponding fluid density is then a constant parameter ρ > 0. We denote by g the gravitational constant.

Boundary conditions. The aquifer lies on an impermeable bed rock at the level z = h bot . We then consider an homogeneous Neumann boundary condition on Γ bot . On the vertical walls Γ ver we consider also an homogeneous Neumann condition to simplify the presentation. Conversely, at the soil level we consider a general Robin condition of parameter (α, β) ∈ R 2 \{0} associated with the source term F = F (t , x). In particular Neumann and Dirichlet conditions can be considered.

Two different kind of flows. We now introduce two auxiliary subregions of Ω in which the flow will present very different behavior. These subregions are based on a function h = h(t , x) which could be one of the unknowns of our model. The subregion above h = h(t , x) is denoted by Ω + h (t ) and the one below is denoted by Ω - h (t ).They are given as follows

Ω + h (t ) := (x, z) ∈ Ω | z > h(x, t ) and Ω - h (t ) := (x, z) ∈ Ω | z < h(x, t ) , (1.10)
Γ h := (x, z) ∈ Ω | z = h(x, t ) . (1.11)
The only condition that this level h have to satisfy is the following

h bot (x) ≤ h(t , x) ≤ max min H (t , x) - P s ρ g , h max (x) , h bot (x) .
(1.12)

As we will see in the model, this condition imposes that the soil is saturated in the lower domain Ω - h (t ) for all t . In practice we will use an explicit characterization of h with respect to the other unknowns. For that we consider a non-negative constant (or function) R and the function

Q = Q(x, H ) Q(x, H ) = max min H - P s + R ρ g , h max (x) , h bot (x) . (1.13)
Although a lot of choices are possible for R, we will focus in the next only on the case of a constant R ≥ 0. We can also remark that it holds Q(x, H (t , x)) = h bot (x) for every (t , x) ∈ ]0, T [×Ω x when R is chosen large enough.

Class of models coupling fast and slow component of the flow.

Each model of the class roughly speaking couples an only vertical flow described by 1d-Richards equations in the upper part Ω + h with a 2D horizontal flow in Ω - h . In this latter part, the vertical flow is assumed to be instantaneous following the classical Dupuit hypothesis.

The coupled model consists in finding pressure head P and velocity v and the auxiliary variables u, w, H and h such that:

• The velocity v in the whole Ω take the form

         v = u + w for t ∈]0, T [ , (x, z) ∈ Ω u = -k r (P ) 1 ρ g ∂P ∂z + 1 K 0 e 3 for t ∈]0, T [ , (x, z) ∈ Ω w = -k r ρ g ( H -z) M 0 ∇ x H for t ∈]0, T [ , (x, z) ∈ Ω (1.14) • In the Ω + h (t ) the following vertical 1d-Richards equation holds              φ ∂s(P ) ∂t + ∂ ∂z u • e 3 = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) α P + β u • e 3 = F for (t , x) ∈]0, T [×Γ soil P = ρ g ( H -h) for t ∈]0, T [ , (x, z) ∈ Γ h (t ) P (0, x, z) = P init (x, z) for (x, z) ∈ Ω + h (0) (1.15)
• In the water table the pressure P satisfies

P (t , x, z) = ρ g H (t , x) -z for t ∈ [0, T [ , (x, z) ∈ Ω - h (t ) (1.16) • The hydraulic head satisfies in Ω x        div x K ( H ) ∇ x H = (u • e 3 )| Γ + h for (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 for (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) for x ∈ Ω x (1.17)
where (u • e 3 )| Γ + h denotes the trace of u • e 3 on Γ h from above. • The interface verify h(t , x) = Q(x, H (t , x)).

(1.18)

The first remark about the model above is that the soil is saturated in Ω - h (t ). More precisely the solution P of problem (1.14)-(1.18) satisfies for every t ∈]0, T [

P ≥ P s in Ω - h (t ). ( 1
P (t , x, z) = ρ g ( H (t , x) -z) ≥ ρ g ( H (t , x) -h(t , x)) ≥ P s .
The second remark is that we have the following result concerning the water flux at the interface Γ h .

(u • e 3 )| Γ + h = (u • e 3 )| h soil + ∂ ∂t h soil (x) h bot (x)
φs(P )d z.

(1.20)

Proof. By integrating first equation of (1.15)

(u • e 3 )| Γ + h = (u • e 3 )| h soil + h soil (x) h bot (x) φ ∂s(P ) ∂t d z = (u • e 3 )| h soil + ∂ ∂t h soil (x) h bot (x)
φs(P )d z.

Thanks to the proposition above, we rewrite (1.17) to get

         (u • e 3 )| h soil + ∂ ∂t h soil (x) h bot (x) φs(P )d z -div x K ( H ) ∇ x H = 0 for (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 for (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) for x ∈ Ω x (1.22)
This formulation implies directly the following result and will also be well adapted to obtain a conservative discrete scheme.

Proposition 1.2. Problem (1.14)-(1.18) is mass-conservative in the sense that:

∂ ∂t Ω φs(P )d x d z = - Ω x (u • e 3 )| h soil d x. (1.23)
In view of reduce ratings we will consider the following problem in which we search only the pressure and the vertical velocity (P, u 3 ) and the auxilliary unknowns H and h satisfying

• In the Ω + h (t ) the following vertical 1d-Richards equation holds                      φ ∂s(P ) ∂t + ∂u 3 ∂z = 0 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) u 3 = -k r (P ) K zz 1 ρ g ∂P ∂z + 1 for t ∈]0, T [ , (x, z) ∈ Ω + h (t ) α P + β u 3 = F for (t , x, z) ∈]0, T [×Γ soil P = ρ g ( H -h) for t ∈]0, T [ , (x, z) ∈ Γ h (t ) P (0, x, z) = P init (x, z) for (x, z) ∈ Ω + h (0) (1.24)
• In the water table the pressure P satisfies

P (t , x, z) = ρ g H (t , x) -z for t ∈ [0, T [ , (x, z) ∈ Ω - h (t ) (1.25)
• The hydraulic head satisfies in

Ω x          u 3 | h soil + ∂ ∂t h soil (x) h bot (x) φs(P )d z -div x K ( H ) ∇ x H = 0 for (t , x) ∈]0, T [×Ω x K ( H ) ∇ x H • n = 0 for (t , x) ∈]0, T [×∂Ω x H (0, x) = H init (x) for x ∈ Ω x (1.26) • The interface verify h(t , x) = Q(x, H (t , x)). (1.27)
It is indeed clear that the solution of (1.14)-(1.18) satisfies (1.24)-(1.27). In the other hand if we have the solution (P, u 3 ) of (1.24)-(1.27) it is possible to get the one of (1.14)-(1.18) by defining v, u and w as solution of (1.14). In this case it holds u • e 3 = u 3 in Ω + h .

NUMERICAL SCHEME

FULLY IMPLICIT SCHEME

In this section, we present a possible conservative numerical scheme for the approximation of problem (1.24)-(1.27) with function Q given in (1.13). The non-trivial task consists in the linearization of the non linear coupling in time. Therefore, we focus in the next only on the time scheme and we keep general the space one.

Time discretization. We begin by performing the time discretization of system (1.24)-(1.27).

We introduce for a fixed M ∈ N * the discrete time t n = n ∆ t with n = {0, ...., M } and ∆ t = T /M . The discrete unknowns at time t n for fixed n ∈ {0, ...., M } are

P n (x, z) P (t n , x, z), H n (x) H (t n , x), h n (x) h(t n , x) and u n (x, z) u(t n , x, z).
By using a backward Euler method to approximate the time derivatives in (1.24)-(1.27), we get the following discrete problem holding for a given n ∈ {1, ...., M } and for known function

P n-1 : finding (P n , u n 3 , H n , h n ) such that                  φ s(P n ) -s(P n-1 ) ∆t + ∂u n 3 ∂z = 0 in (x, z) ∈]h n , h soil [ u n 3 = -k r (P n ) K zz 1 ρ g ∂P n ∂z + 1 in (x, z) ∈]h n , h soil [ α P n + β u n 3 = F n on Γ soil P n x, h n (x) = ρ g H n (x) -h n (x) on Ω x (2.1)            u n 3 | h soil + 1 ∆t h soil (x) h bot (x) φs(P n )d z - h soil (x) h bot (x)
φs(P n-1 )d z

-div x K ( H n ) ∇ x H n = 0 in Ω x K ( H n ) ∇ x H n • n = 0 on ∂Ω x (2.
2)

P n (x, z) = ρ g H n (x) -z in (x, z) ∈ Ω - h n (t ) (2.3) h n (x) = Q(x, H n (x)) in Ω x . (2.4)
The first result about this implicit scheme is following mass-conservation property.

Proposition 2.1. Let P n-1 a given function and

(P n , u n 3 , H n , h n ) the solution of problem (2.1)- (2.4). Then it holds in ]0, T [×Ω x : Ω φ s(P n ) -s(P n-1 ) d x d z = -∆ t Ω x u n 3 | h soil d x. (2.5)
Proof. The result is obtained by integrating the first equation of (2.2) over Ω x and by using an integration by part formula and the boundary condition in (2.2). Notice that the choice of a fully implicit scheme in the discretization of problem (1.24)-(1.27) is not only done to obtain a conservative scheme. In fact, a real difficulty is to deal with equation (2.2) which can become an ill-posed if the first three terms are treated as given source term constant with respect to H n . We precise this point in Remark 2.3 bellow.

Reformulation with a Dirichlet to Neumann operator.

To treat the non-linear coupling of systems (2.1)-(2.4) we start by giving a more convenient formulation. Then we introduce the following function Θ(x, H , P ) defined for x ∈ Ω x , H ∈ R and P :]h bot (x),

h soil (x)[ → R by Θ(x, H , P ) = u 3 | h soil + 1 ∆t h soil (x) h bot (x)
φs(P )d z -

h soil (x) h bot (x) φs(P )d z , (2.6)
where (P, u 3 ) is the solution of the following problem R(x, H , P ) :

                     φ ∆t s(P ) -s(P ) + ∂u 3 ∂z = 0 in ]h, h soil (x)[ u 3 = -k r (P ) K zz 1 ρ g ∂P ∂z + 1 in ]h, h soil (x)[ α P + β u 3 = F for z = h soil (x) P = ρ g (H -z) in [h bot , h] h = Q(x, H ). ( 2
, h n , H n ) of (2.1)-(2.4) satisfies equivalently      P n (x, •), u n 3 (x), h n (x) solution of R x, H n (x), P n-1 (x, •) a.e. in Ω x Θ x, H n (x), P n-1 (x, •) -div x K ( H n )∇ x H n = 0 in Ω x K ( H ) ∇ x H n • n = 0 on ∂Ω x (2.8)
In the above formulation, the second equation can be seen as a mass conservation equation associated with the velocity

w n := -K ( H n )∇ x
H n and with the evolution of the volume given by Θ. More precisely, the quantity 1 ∆t h soil (x)

h bot (x) φs(P n )d z - h soil (x)
h bot (x) φs(P n-1 )d z represents the evolution of the total volume of water contained in the column ]h bot (x), h soil (x)[ between the two time steps t n-1 and t n . This evolution is non vanishing when some water flows out, or flows in, the coloumn through the soil level h soil (x) or in the horizontal direction. Due to the presence of the term u n | h soil (x) in the the definition of Θ we have that the value Θ(x, h n , H n , P n-1 ) characterizes the water which flows out (or flows in) from the coloumn in the horizontal direction only.

Remark 2.2. Let (P n , u n , h n , H n ) solution of (2.1)-(2.4) and x ∈ Ω x . By subsituting equation (2.1) in the definition of Θ(x, H n (x), P n-1 (x, •)), we can easily notice that it holds:

Θ(x, H n (x), P n-1 (x, •)) = u n 3 | h n + h n (x) h n-1 (x)
φ(1 -s(P n-1 ))d z.

(2.9)

In particular in the case where h

n (x) = h n-1 (x) =: h(x) we have Θ(x, H n (x), P n-1 (x, •)) = u n 3 | h n . In this case Θ(x, H n (x), P n-1 (x, •))

is the Dirichlet to Neumann operator which associate the flux u n

3 | h at z = h(x) obtained from the Richards equation with the Dirichlet condition at z = h(x) given by H n . This append for exemple when h n (x) = h n-1 (x) = h bot (x) (if R is large enough and/or if there is not a lot of water at the bottom of the aquifer) and when h n (x) = h n-1 (x) = h max (x) (which can occur if the aquifer overflows and R small enough).

Another remark about formula (2.9) is that it is the discrete version of equality (1.1) holding for the continuous problem. Notice in particular that appears in the discrete formula the correcting term h n (x) h n-1 (x) φ(1 -s(P n-1 ))d z which is necessary to have a conservative scheme.

Fixed point strategy. At this stage, the discrete problem consists in finding (P n , u n , h n , H n ) solution of the nonlinear problem (2.8). We use for that a fixed-point method to linearize the problem. It consists in building the following sequence (P n k , u n k 3 , h n k , H n k ) which, when converges, tends to (P n , u n , h n , H n ). For fixed n ∈ N * and known pressure P n-1 we define

• Initialization: H n 0 = H n-1 • Heredity: for all k ∈ N * , we set (P n k , u n k 3 , h n k , H n k ) as solution of the linear problem            P n k (x, •), u n k 3 (x), h n k solution of R x, H n k-1 (x), P n-1 (x, •) a.e. in Ω x Θ(x, H n k-1 , P n-1 ) + Λ(x, H n k-1 , P n-1 )( H n k -H n k-1 ) -div x K ( H n k-1 )∇ x H n k -div x K ( H n k-1 )( H n k -H n k-1 )∇ x H n k-1 = 0 on Ω x K ( H n k-1 ) ∇ x H n k • n = 0 on ∂Ω x (2.
10) where K is the derivative of K with respect to P .

To use this procedure we need to chose a convenient stabilization function Λ = Λ(x, H , P ). Notice that, by construction, if Θ is enough regular with respect to H we can make the procedure (2.10) to be exactly the Newton method. In this case we only have to set Λ(x, H , P ) := ∂Θ ∂H (x, H , P ). Of course there is no a priori reason for Θ to be enough regular and a good choice for Λ(x, H , P ) will be crucial to make the procedure (2.10) convergent. This is the object of the next subsection.

Remark 2.3. The simplest idea to deal with the procedure (2.10) is to make the choice Λ = 0. It turns out that it is not possible since problem (2.10) become ill-posed in this situation (many possible solutions due to the Neumann boundary condition). In the other hand, if the boundary condition in (2.10) is changed to a Dirichlet one, the problem become well posed and the sequence (H n k ) k can be considered. Nevertheless this sequence is non-converging in practice and cannot be used to get H n solution of (2.8).

CHARACTERIZATION OF FUNCTION Λ

In this subsection we investigate the dependence of the function Θ with respect to H . The goal being to propose a function Λ making the procedure (2.10) to converge.

In the next of this subsection, we consider a function ,H ). We introduce also (P, u) the solution of problem R(x, H , P ). We start by giving the following new formulation of Θ(x, H , P ) given in (2.6): we split the first integral into two parts and by using the fact that s(P ) = 1 in ]h bot (x), h[ (see (1.19)) we get

P : [h bot (x), h soil (x)] → R, x ∈ Ω x , H ∈ R and h = Q(x
Θ(x, H , P ) = u 3 | h soil + h h bot (x) φ ∆t d z + h soil (x) h φ ∆t s(P )d z - h soil (x) h bot (x) φ ∆t s(P )d z.
(2.11)

In particular we remark that the definition of Θ is obtained as soon as we know (P, u) in the upper part ]h, h soil (x)[. We then consider (P, u) as solution of the following vertical Richards equation (discretized in time).

C (x, H , P )                  φ ∆t s(P ) -s(P ) - ∂u 3 ∂z = 0 in ]h, h soil (x)[ u 3 = -k r (P ) K zz 1 ρ g ∂P ∂z + 1 in ]h, h soil (x)[ P = ρ g (H -h) for z = h α P + β u 3 = F for z = h soil (x)
(2.12)

At this stage we see from (2.11) that the dependence of Θ with respect to H involves the one of the solution (P, u 3 ) of problem C (x, H , P ) above. Moreover, we see in (2.12) two different kind of dependence with respect to H (and h := Q(x, H )):

• one holding through the value of the boundary condition on z = h • one holding through the position of this boundary itself These two cases are treated separately in the two next paragraphs.

Stabilization term in the constrained situation. In this part, we assume that (x, H , P ) is such that

H < h bot (x) + P s + R ρ g or H > h max (x) + P s + R ρ g . (2.13)
Accordingly, thanks to (1.13), we have h := Q(x, H ) ∈ {h soil (x), h max (x)}. In particular for any ε > 0 small enough we have Q(x, H ± ε) is constant with respect to ε. We then denote in the next of this part h cons = Q(x, H ) to insist on the fact that Q(x, H ) is constant with respect to small variations of H . In this situation, the problem (2.12) becomes C cons (x, H , P ) :

                 φ ∆t s(P ) -s(P ) - ∂u 3 ∂z = 0 in ]h cons , h soil (x)[ u 3 = -k r (P ) K zz 1 ρ g ∂P ∂z + 1 in ]h cons , h soil (x)[ P = ρ g H -h cons for z = h cons α P + β u 3 = F for z = h soil (x) (2.14)
and Θ satisfies

Θ(x, H , P ) = (u 3 )| h soil + h cons h bot (x) φ ∆t d z + h soil (x) h cons φ ∆t s(P )d z - h soil (x) h bot (x) φ ∆t s(P )d z.
It is important to notice that in this constrained situation, the dependence of the solution (P, u 3 ) of problem (2.

                 φ ∆t s (P )P - ∂u 3 ∂z = 0 in ]h cons , h soil (x)[ u 3 = -K zz k r (P )P 1 ρ g ∂P ∂z + 1 + k r (P ) 1 ρ g ∂P ∂z in ]h cons , h soil (x)[ α P + β u 3 = 0 on Γ soil P = ρ g on Γ bot (2.16)
In the constrained situation where (2.13) holds an thanks to (2.15) and (2.16), we will define Λ = Λ(x, H , P ) by

Λ(x, H , P ) = u 3 | h bot + 1 ∆t h soil (x) h bot (x)
φs (P )P d z (2.17)

where (P , u ) solution of R (x, H , P ) and P solution of R(x, H , P ).

Stabilization term in the free situation. In this part, we consider x ∈ Ω x , a function P :

[h bot (x), h soil (x)] → R and H ∈ R such that h bot (x) < H - P s + R ρ g < h soil (x). (2.18)
Accordingly, thanks to (1.13), we have h

:= Q(x, H ) ∈]h soil (x), h max (x)[ and h = H - P s + R ρ g . ( 2.19)
In this situation, the problem (2.12) becomes

                 φ ∆t s(P ) -s(P ) - ∂ ∂z (u 3 ) = 0 in ]h, h soil (x)[ u = -k r (P ) 1 ρ g ∂P ∂z + 1 K 0 e 3 in ]h, h soil (x)[ P = P s + R for z = h α P + β u 3 = F for z = h soil (x) (2.20)
Conversely to the previous case, the dependence of (P, u) solution of problem below holds only through the position of the interface z = h. The Dirichlet condition at z = h being now the constant P s + R. This dependence is more difficult to handled. The strategy here is to use a quasi-Newton procedure in which Λ only describe a part of the dependence of Θ with respect to H (Θ defined in (2.6)). We introduce

Θ 1 (x, H , P ) := h h bot (x) φ ∆t d z (2.21) Θ 2 (x, H , P ) := h soil (x) h φ ∆t s(P )d z - h soil (x) h bot (x) φ ∆t s(P )d z + u 3 | h soil It comes Θ(x, H , P ) = Θ 1 (x, H , P ) + Θ 2 (
x, H , P ) and we choose Λ in (2.10) defined formally by

Λ(x, H , P ) = ∂Θ 1 ∂H (x, H , P ).
Finally, Λ in this situation will be given by

Λ(x, H , P ) = φ| h ∆t . ( 2.22)
With this choice, we remark that Λ(x, H , P ) is in fact constant with respect to P and is explicitly obtained from H .

Definition and comments on the stabilization term. In view of (2.13), (2.17), (2.18) and (2.22) we have the following characterization of Λ.

Definition 2.4. Let x ∈ Ω x , H ∈ R, P a function of z ∈]h bot (x), h soil (x)[ and h := Q(x, H ).
We denote (P, u 3 ) the solution of problem (2.12) and (P , u 3 ) the solution of problem (2.16). The function Λ in (2.10) is defined by

Λ(x, H , P ) :=        u 3 | h soil (x) + 1 ∆t h soil (x) h φs (P )P d z if h ∈ h bot (x), h soil (x) φ| h ∆ t if h ∈ h bot (x), h soil (x) (2.23)
The first remark about the definition of Λ above concerns the free situation h ∈ h bot , h max . As seen in the previous paragraph, the proposed value is obtained by taking into account only the dependence of the component Θ 1 with respect to H (Θ 1 defined in (2.21)). Accordingly, the dependence of Θ 2 with respect to H is neglected.

This type of choice is classical and the corresponding procedure (2.10) is a quasi-Newton one. Generally, the main interest is to simplify the characterization of the stabilization term Λ at the price of a possibly reduction of the order of convergence. The lost of convergence is in general minimized when the neglected part (here Θ 2 ) has a non-dominant dependence with respect to H .

It is nevertheless important to notice that in our case, the dependence of Θ 2 with respect to H has no a priori reason to be non-dominant. In fact, in several cases it turns out that the function Θ 2 almost compensates the dependence of Θ 1 with respect to H . In this case, the function Θ is almost constant with respect to H and the procedure (2.10) does not converge.

ALGORITHMS Algorithms for the 1d vertical problems.

In this part we construct the functions which computes the solutions of problems (2.12)-(2.16) and the values of Θ and Λ given in (2.6) and (2.23) respectively. We then fix x ∈ Ω x , a function P :]h bot (x), h soil (x)[ → R and a real number H .

We give in Algorithm 1 the function RICHARD1D. It computes the solution of problem R given in (2.7) by using a classical Newton method. The tolerance is ε and the maximum number of iterations for the Newton method is N max . This Algorithm is the same for any values of h independently of the case h ∈]h bot (x), h max (x)[ or h ∈ {h bot (x), h max (x)}. Notice that in practice the solution of (2.24) is done in a discrete way. The associated vertical discretization has in general a fixed space step (except eventually for the lowest). Accordingly the number of space steps depends on the level h and then decrease when h increase. The extremal situation being when h reach the maximum value h = h max and in which the number of space step is the smallest. In particular this number can be equal to one if the parameter δ in (1.4) coincide with the highest space step. If a low order approximation method is used for the

Algorithm 1 function RICHARDS1D(x, H , P ) h ← Q(x, H ), P ← P | (h,h soil (x)) δ ← ε + 1, n ← 0 while (δ > ε) & (n < N max ) do P o ← P , n ← n + 1 (P, u 3 ) unique solution of the linear problem                  φ s(P o ) -s(P ) + φ s (P o )(P -P o ) -∆t ∂u 3 ∂z = 0 in ]h, h soil (x)[ u 3 K zz = -k r (P o ) 1 ρ g ∂P ∂z + 1 -k r (P o ) 1 ρ g ∂P o ∂z + 1 (P -P o ) in ]h, h soil (x)[ P = ρ g (H -h) for z = h α P + β u 3 = F for z = h soil (x)
(2.24)

δ ← P -P o ]h,h soil (x)[ end while if (δ > ε) then Er r ← 1 end if P | (h,h soil (x)) ← P , P | (h bot (x),h) ← ρ g (H -z), u 3 | (h,h soil (x)) ← u 3 , u 3 | (h bot (x)
,h) ← 0 return (P, u 3 , Er r ) end function vertical approximation (for example P 1 finite element method), the solution of (2.24) can be explicit. In the other hand it is classical that the Newton method is not necessarily convergent for solving the highly degenerate Richards problem (even with high order method). In particular the returned variable Er r indicate if yes or no the method has converged. The strategy in the case of a non-converging sequence consists in dividing the time step by 2. This point is precised in Algorithm 6.

The function THETA is described in Algorithm 2. Notice that in practice we chose the inputs as being (x, P, u 3 , P ) instead of (x, H , P ) as it is the case in (2.6). Indeed, the dependence of Θ with respect to H and P holds through the solution (P, u 3 ) of (2.7). It turns out that this solution (P, u 3 ) will be not exclusively used to compute the value Θ and, since (P, u 3 ) is computationally expensive to obtain, it is crucial to solve (2.7) only when it is required. Algorithm 2 function THETA(x, P , u 3 , P )

Θ ← u 3 | h soil + 1 ∆t h soil (x) h bot (x) φs(P )d z - h soil (x) h bot (x) φs(P )d z ,

return Θ end function

We give in Algorithm 3 the function DERIV_RICHARDS which computes the solution (P , u 3 ) of problem (2.16). In particular those solutions are only defined on ]h cons , h soil (x)[. The procedure is more simple than the one in Algorithm 1 (does not need Newton fixed-point iterations) since the problem is linerar. The inputs are the horizontal position x ∈ Ω x , the constrained level h cons and the pressure P obtained as solution of problem (2.7) by function RICHARDS1D. Of course in practice, as for (2.24), the solution of (2.25) is search in a discrete Algorithm 3 function DERIV_RICHARDS(x, h cons , P ) (P , u 3 ) unique solution of the linear problem

                 φ s (P ) P -∆t ∂u 3 ∂z = 0 in ]h cons , h soil (x)[ u 3 = -k r (P ) K zz ρ g ∂P ∂z -k r (P ) K zz 1 ρ g ∂P ∂z + 1 P in ]h cons , h soil (x)[ P = ρ g for z = h cons α P + β u 3 = F for z = h soil (x) (2.25)
return (P , u 3 ) end function way. The resolution of this problem is then remplaced by the solution of a finite linear system. It is important to notice that some matrices (e.g. associated to s (P ) and k r (P )) which characterize this linear system are in fact already computed during the solving of (2.24). The natural strategy to save time in the numerical resolution is then to make the function RICHARDS1D to return also those matrices.

We introduce in Algorithm 4 the function LAMBDA which returns the value Λ(x, H , P ) following the definition (2.23). Notice that in practice we chose the inputs as being (x, P, h) instead of (x, H , P ) proposed in (2.23). Indeed, the dependence of Λ with respect to H and P holds through the solution P of (2.7). Since this solution P is computationally expensive to obtain, it is crucial to solve (2.7) only when it is required.

Algorithm 4 function LAMBDA(x, P , h) if h bot (x) < h < h soil (x) then Λ ← - φ| Γ h ∆ t else P , u 3 ← DERIV_RICHARDS x, h, P Λ ← u 3 | h soil (x) + 1 ∆t h soil (x) h φs (P )P d z end if return Λ end function
Algorithm for the 2d horizontal problem. In this part, we precise the numerical strategy to obtain the solution (P n , u n , H n ) of problem (2.8) which reads as

     P n (x, •), u n (x), h n (x) solution of R x, H n (x), P n-1 (x, •) a.e. in Ω x Θ x, H n (x), P n-1 (x, •) -div x K ( H n )∇ x H n = 0 in Ω x M 0 ∇ x H n • n = 0 on ∂Ω x Algorithm 5 function HORIZONTAL_PROBLEM(H ol d , P ol d ) H ← H ol d , δ ← ε + 1, k ← 0 I ← Ω x N Er r ← 0 while (δ > ε) & (k < N max ) do k ← k + 1 H p ← H for all x ∈ I do
This loop can be done parallelly

P (x), u 3 (x), Er r ← RICHARDS1D x, H p (x), P ol d (x) N Er r ← N Er r + Er r Θ(x) ← THETA x, P (x), u 3 (x), P ol d (x) Λ(x) ← LAMBDA x, P (x),Q(x, H p (x)) end for H solution of the linear problem      Θ + Λ(H -H p ) -div x K (H p )∇ x H -div x K (H p )(H -H p )∇ x H p = 0 on Ω x M 0 ∇ x H • n = 0 on ∂Ω x δ ← H -H p Ω x I ← x ∈ Ω x , |H (x) -H p (x)| > ε 3 end while if (δ > ε) or (N Er r > 0) then Er r ← 1 else
Er r ← 0 end if return (H , P, Er r ) end function and where P n-1 (x, •) is a given function of ]h bot (x), h soil (x)[ for every x ∈ Ω x . The the corresponding procedure is given Algorithm 5 whose the inputs are functions H ol d : Ω x → R and

P ol d : Ω → R.
It is important to notice the presence of the evolutive set I is this function. In practice this set collects the horizontal abscissa for which it is needed to update the value of Θ and Λ. In facts in many cases, the convergence will be obtained quickly in a part of the domain Ω x and would need more iteration to be reached in another part. In particular, by continuity, if for any x ∈ Ω it holds H (x) ∼ H p then we expect that Θ(x, H (x), P ) ∼ Θ(x, H p (x), P ) and the same for Λ. Is this case we choose to not compute the new value Θ(x, H (x), P ) and to keep Θ(x, H p (x), P ).

Another strategy to save computation time is to take advantage of the fact that every computation in the for loop over Ω x do not depend each other. A good idea, of simple implementation, is to make those computations in parallel.

In the other hand, as it is the case in Algorithm 1, the tracking of possibly non converging sequence in this fixed point procedure is crucial. Indeed, we a priori do not have any result proving theoretically this kind of convergence. The procedure then return the variable Er r which will be equal to 0 if the sequence effectively converges and 1 in the converse case (if the fixed point and/or if any 1d-Richards fail). As previously, if Er r = 1 at the end of this function, the strategy is again based on the reducing of the time step (see Algorithm 6 below).

Algorithm 6 function EVOLUTIVETIME(H 0 , P 0 ) for all n ∈ {1, . . . , M } do (H , P ) ← (H n-1 , P n-1 ) ∆ e t = ∆ t count ← 0 n t ← 1 while count < n t do count ← count + 1 (H , P, Er r ) ← HORIZONTAL_PROBLEM(H , P, ∆ e t ) if Er r = 1 then n t ← 2 n t ∆ e t = ∆ t /n t count ← 0 (H , P ) ← (H n-1 , P n-1 )
end if end while (H n , P n ) ← (H , P ) end for return (H , P ) end function Algorithm for the time approximation. We present in Algorithm 6 the procedure which construct iteratively the discrete solution (P n , H n ) of problem (2.1)-(2.4). The strategy is obious and consist in looping over the time n and calling at each step the function HORI-ZONTAL_PROBLEM which computes the corresponding solution. The only problem is when an non converging sequence is obtained during the solving of any 1d-Richards problem or during the fixed point procedure in the HORIZONTAL_PROBLEM. As said before a classical strategy to overcome this kind of problem is to reduce the time step. It is what do the while loop. In particular, the time step could be divided by two more than one time if needed.

Notice that this strategy of evolutive time step is not optimal. Indeed in the presented Algorithms, if only one 1DRICHARDS is non converging, the time step will be reduce to increase the chance of convergence but not only for this problematic 1DRICHARDS but for all of them. In particular the fail of one 1DRICHARDS leads to the computation (in the best case) 2 times Ω x one. It is then quite obvious that a better idea is the use this evolutive time step strategy also at the level of each 1DRICHARDS procedure.

NUMERICAL RESULTS

In this section we compare quantitatively the flow characterized by models of the class (1.14)-(1.18) with a reference flow obtained from a more gerneral and precise model. The later one is the following 3d-Richards problem. Its unknowns are the pressure P and the velocity vector field v. By considering the same geometry and physical description of the soil as in the section 1, the problem is to find (P, v) solution of

                 φ ∂s(P ) ∂t + div(v) = 0 in ]0, T [×Ω v = -k r (P ) K 0 1 ρg ∇P + e 3 in ]0, T [×Ω α P + β v • n = F on ]0, T [×Γ soil v • n = 0 on ]0, T [×(Γ bot ∪ Γ ver ) (3.1)
Notice that in fact the class of model (1.14)-(1.18) has been obtained form this general 3done. The idea being to characterize the dominant flow arising at different time scale when the aquifer has a large horizontal dimension with respect to its deepness. This is the purpose of [2]. Accordingly, since the coupled model is an approximation of the 3d-Richards one, it is crucial to quantify the difference between them. It is the purpose of this section.

In this article we want to compare the evolutions of the water table obtained respectively through problem (3.1) and (1.14)-(1.18) for several choice of Q. What we understand by water table is the saturated part of the soil that relies on the impermeable layer at the bottom of the aquifer. To describe it we introduce the following function h sat and the set Ω - h sat (t ) defined for a given pressure P = P (t , x, z).

h sat (t , x) := sup I t ,x , I t ,x := z ∈ [h bot (x), h max (x)] | P (t , x, z ) > P s , ∀z ∈ [h bot (x), z[ , ( 3.2)
Ω - h sat (t ) := (x, z) ∈ Ω | z < h sat (t , x) . (3.3)
By construction and if P is continuous we have

P t , x, h sat (t , x)      = P s if h bot < h sat < h max ≥ P s if h sat = h max ≤ P s if h sat = h bot and P (t , x, z) ≥ P s for all z ∈]h bot , h sat ].
In particular the soil is fully saturated in Ω - h sat (t ) for every t ∈]0, T [. In other words, h sat is the iso-pressure {P = P s } as soon as the aquifer is not overflowing (in this case h sat = h max ) and non empty (in this case h sat = h soil ) at the position x ∈ Ω.

Notice also that in the model (1.14)-(1.18) the function h represent by construction the iso-pressure {P = P s + R} as soon as the aquifer is not overflowing (in this case h sat = h max ) and non empty (in this case h sat = h soil ) at the position x ∈ Ω. In particular it holds h sat = h if and only if R = 0 in (1.13).

To simplify the numerical computations and the visualizations of the results, all the experiments of this paper will be done in a two-dimensional setting (Ω x being an interval as in Figure (1)).

We will consider two situations describing the water recharge of the aquifer. The difference between them will be the origin of recharging water. In the first experiment we will consider a water reservoir above the water table. In the second one the water will come from the overland (from a river, lake ...) through the boundary Γ soil and characterized by the Dirichlet condition. We decompose the study of those different cases in the two next subsections.

Physical parameters.

In each experiment, the soil will be characterized by the following set of data. In particular the Brooks and Correy model (1.6) is used (see [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF]). Denoting I 3 the 3×3 identity matrix we set:

s(P ) = (P s /P ) λ , k r (P ) = (P s /P ) 2+3λ , (P s , λ) = (-1.5, 3), ρ = 1, φ = 0.1, K 0 = 0.1 I 3 .
Moreover, the parameter δ in (1.4) is chosen as small as possible, that is equal to the size of one vertical mesh.

Comparaison between models.

We are going to focus on the two following parameters which will influance the flow and then the precision of (1.14)-(1.18) with respect to (3.1). Those parameters are

• The geometrical ratio horizontal length/deepness of the aquifer. In particular, as explian in [2], we expect that problems (1.14)-(1.18) approximate well (3.1) when this ratio is small. The remaining question being how small is needed to have a good approximation. • The parameter R which characterize the model in the class (1.14)-(1.18). Therefore, for each experiment, we will consider a large number of situations by varying these parameters.

FIRST EXPERIMENT: UNCONFINED AQUIFER

Geometry, initial and boundary conditions. In this first experiment, the geometry of the aquifer is a simple rectangle Ω :=]0, L x [×]h bot , 0[ where L x > 0 is the horizontal length and where the top and bottom boundaries are characterized by constant functions h soil (x) = 0 and h bot (x) = h bot < 0. The precise values for L x and h bot will be changed to observe their influences on the flow (described by (1.14)-(1.18) or (3.1)), and more precisely, to quantify the precision of model (1.14)- (1.18).

At time t = 0, we consider a setting where the function h 0 sat introduced in (3.2) corresponds to the top level of the water table. We choose a constant function h sat (0, x) = h 0 sat ∈]h bot , h soil [. The corresponding initial pressure P (0, •, •) is choosen in the stationnary state P (0, x, z) = ρ g (h sat (0, x) -z) + P s for all (x, z) except near a rectangular region above z = h sat where the pressure goes smoothly to the saturation value P s . The saturated part in this rectangle plays the role of a reservoir which will infiltrates downward due to the gravity. The rectangle is

ω =]L x /10, 3L x /10[× ]-3.5, -1.7[ . (3.4)
This initial situation is drawn in the first picture of Figure 2. In every picture the gray scale corresponds to the saturation value, the maximal darkness corresponding to s 1. This one is obtained in the particular case:

h 0 sat = -3.6, L x = 40, and
h bot = -5.
In this experiment we are only interested on the infiltration coming from the supply ω. For this reason we consider homogeneous Neumann condition at the soil level (α = F = 0 in (1.15)). At time t = 1.425 the water initially in supply ω has started to flow down due to the gravity. We notice that the water flows gobally in the vertical direction. In particular the initially horizontal front is still present. This front in about to reach the water table.

At time t = 3.25 the infiltrating water has touched the water table. The top level of the latter, that is the level z = h sat , has increased. From this moment, we observe a more significative horizontal flow since the level z = h sat increase also for x which are not directly below the supply ω.

After some time almost all the water initially located in the rectangle supply have reached the water table. Then the interface h sat becomes flater and is associated with a pressure admitting the stationary profile P (t , x, z) = P s + ρ g h sat (t , x) -z . Error v.s. ratio: in the unsaturated part. We start by comparing the flow in the unsaturated part of the aquifer obtained through the 2d-richard problem or through the coupled model. In particular, we want to see its dependence with respect to the ratio L x |h sat | . For this purpose, we are interested on the the time it takes for the infiltrating water to reach the water table. This quantity will gives an precision about the vertical flow hypothesis which is done in Ω + h in the coupled problem (1.14)- (1.18). Notice that for this test, it is not usefull to precise the value R in (1.13) since this one do not influance the flow in Ω + h . We introduce t r the time for wich the water initially in the water supply ω is about to reach the water table. It is characterized in this situation by the time for which the level of the saturated part h sat start to increase. We introduce then t r 2d and t r c respectiveley for this reached time obtained by using the 2d-Richards problem or the coupled problem.

We represent in the left drawing of Figure 3 the evolution of relative difference |t r 2d -t r c |/t r 2d in terms of the ratio L x |h sat | . The horizontal length of aquifer L x is chosen such that L x ∈ 20, 30, . . . , 100, 120, 140, 160 , so that the ratio L x |h bot | vary in [4,32]. We see that this difference is quite small even for the ratio of 4 which is far from being a small one. Moreover, as expected, this difference globally decrese when the ratio increase.

To precise this result we represent the saturation profile which is obtained exactly at this time t r c . It is done for L x = 40 and t r c = 1.425 days in the drawing in the right of Figure 3. This profile is to compare with the corresponding one in Figure 2 physical parameters. A we can see, those two profile are very close each other and present for exemple the same horizontal front. Moreover we see that those front are nearly as close for the saturated water table (in agreement with the small error |t r 2d -t r c |/t r 2d ). In conclusion, it seems that the horizontal flow in this situation is neclectable even for not so shallow aquifer. The hypothesis of only vertical flow in Ω + h seems to be well adapted in this situation of a free infiltration case. To complete the investigation of this case, we compare the evolution of the water table after that the infiltrating water flows inside it. The influance of parameter R is also crucial to took into account.

Dependence with respect to R. In this part we compare the solution of the classical Richards model with the one obtained by using the coupled model (1.14)- (1.18). We test four particular

L x h 2d sat -h R sat L 1 ([0,T ]×Ω)
versus the ratio length/deepness of the aquifer (R ∈ {0, 4, 8, ∞}). Function h 2d sat is the iso-pressure P = P s in the original 2d-Richards problem and h κ sat is the one associated with the coupled problems for different choices of R characterizing h by (1.13). The test case is the one of Figure 2 with h 0 sat = -2.2.

h soil is given by

h soil (x) = - 3 2 e -Lx 20 (x-Lx 2 ) 2 .
This hollow in the soil level may corresponds to a riverbed or a basin/lake filled on water. We represent the geometry in pictures of Figure 6. As in the previous case, the values of L x and h bot will be changed to observe their influences on the flow.

Overland water. We consider the hole in the soil level to be filled by water, e.g. due to the presence a river or a lake. We do not describe precisely the overland flow in this paper to focus on the underground one. We then consider the simplest model for it. Let us consider a function h riv = h riv (x) which characterizes the top level of the river. The overland water occupy the region

Ω riv := (x, z) ∈ Ω x × R | z ∈]h soil (x), h riv (x)[ . (3.5)
To simplify the description we assume that the surface of the overland water is horizontal. Accordingly h riv is constant with respect to x ∈ Ω x . In the other hand, the water is assumed to be in the hydrostatic state. The corresponding water pressure satisfies for all (x, z) ∈ Ω riv

P (x, z) = ρ g (h riv -z) (3.6)
Boundary conditions. To take into account the overland water, we split the boundary Γ soil into two the parts Γ D soil and Γ N soil with

Γ D soil = (x, z) ∈ Γ soil | h soil (x) ≤ h riv and Γ N soil = (x, z) ∈ Γ soil | h soil (x) > h riv .
The region below the river is then Γ D soil and we will consider a Dirichlet condition over it. According to the equation (3.6) we choose the value to be given by P (x, z) = ρ g (h riv -h soil (x)) on Γ D soil . In the other hand, and to simplify the modeling, we consider an homogeneous Neumann condition on Γ N soil . Initial conditions. At time t = 0, we consider a setting where the function h 0 sat introduced in (3.2) corresponds to the top level of the water table. We choose a constant function h sat (0, x) = h 0 sat ∈]h bot , h soil [. The corresponding initial pressure P (0, •, •) is chosen in the stationary state P (0, x, z) = ρ g (h sat (0, x)-z)+P s for all (x, z) except near the river Γ D soil where the pressure goes smoothly to the value of the Dirichlet condition ρ g (h riv -h soil (x)). Notice that the boundary condition at the soil level where there is no overland water is, as before, an homogeneous Neumann condition.

Impermeable rock

This initial situation is drawn in the first picture of Figure 6. As before, in every picture the gray scale corresponds to the saturation value, the maximal darkness corresponding to s 1. This one is obtained in the particular case: At time t = 1.215 the water of the river Ω riv has started to infiltrates and to flow down due to the gravity. In contrary of the previous case, the satureted infiltration zone is a little larger than the river. This indicates that, even if the vertical flow is dominant, the horizontal one is not as small as in the previous case. The hypothesis of vertical flow in Ω + h made by the model (1.14)- (1.18), may be less valid in this case. At this time, the inflitration front in about to reach the water table.

At time t = 2 the infiltrating water has reached the water table. The top level of the latter, that is the level z = h sat , has increased a lot. This increase is faster than is the first experiment since the whole region below the river is already saturated before reaching the water table.

From this moment, we observe a more significant horizontal flow since the level z = h sat increase also for x which are not directly below the river. Error v.s. ratio. As in the first experiment, we start by comparing the flow in the unsaturated part of the aquifer obtained through the 2d-Richards problem or through the coupled model.

We denote again t r the time for which the water infiltrating from the river reach the water table. As before it is characterized here by the instant when the level of the saturated part h sat start to increase. We introduce then t r 2d and t r c respectiveley for this reach time obtained by using the 2d-Richards problem or the coupled problem.

We represent in the top left drawing of Figure 7 . Evolution of the iso-pressure P = P s obtained from the classical Richards equation (h 2d sat ) and from the coupled model for three choices of h given by R = 0, R = 3 and h = h bot (h κ sat for κ ∈ {a, b, c} respectively). The test case is the one of Figure 6.

In contrary of the first experiment, there is an important difference for small ratios. This can be explained by the presence of the horizontal flow in the general case (see second drawing of Figure 6).

To precise this result we represent the saturation profile which is obtained exactly at this time t r c . It is done for L x = 40 and t r c = 1.215 days in the drawing in the top right of Figure

It is also important to notice that in this case, the total mass of water which enters the aquifer depends also on the choice R. This is due to the Dirichlet condition which allows different fluxes at the soil level. In particular it is not only the profile of the water table but also the mass transfer from the overland to the aquifer which depends on the choice of R. Of course, this dependence is less significant when the ration L x |h bot | increases as we will see in the next part.

Error v.s. ratio: whole space-time. We represent in Figure 9 the error 1

L x h 2d sat -h R sat L 1 ([0,T ]×Ω)
for R ∈ {0, 4, 7, ∞} and for L x ∈ [20,160]. The initial water table is given by h 0 sat = -3.6. Again, this error globally decrese when the ration L x /|h bot | increase and good choice of the parameter R provide a reduction of the error. Nevertheless, the situation is not as simple as in the first situation. Ideed a significant dependence of the optimal value R apears with respect to the ratio. In this case, the value R = 7 seems to be a good choice. This is the same value as the one obtained in the first experiment. To see if this optimal value is or not dependent on the situation, more experiment will be necessary. The final model consists of a coupled system of parabolic-type pde's that can be degenerated according to the degeneration of the moisture content. We prove a result of existence of weak solutions in both cases non-degenerate and degenerate.

INTRODUCTION

More and more often populated areas are affected by contamination of soil and groundwater. Many modeling approaches are developed for studying the vulnerability of aquifers with regard to agricultural pollutions, with a particular focus on the nitrates input. There is an abundant literature on each of the involved processes (mathematical, physical, ...), so that we can consider that a "realistic" modeling already exists. Nevertheless there is a wide variety of processes (chemical, hydrogeological, anthropic, ...) acting in a wide range of temporal and geometrical length scales. It follows that the assembly of the corresponding model bricks, if considered like toolboxes of a software, is, at best, computationally expensive.

It should be noted that the main concern for the derivation of the hydrogeological model is related to its good description of the flow between the ground level (the level of the anthropic processes) and the water table. This will be crucial when studying the transport of chemical components in the aquifer. Indeed, it turns out that many chemical reactions are expected in the first meters of the subsoil, where oxygen is still very present. In particular, chemical species that reach the water table are not necessarily the same as those that have left the surface. This yields different speeds of the reactive kinetics. As a result, for an efficient mathematical modeling, the time upscaling process in this zone must keep track of all the time scales.

In this paper, we focus on the hydrogeological question. Aquifers are often characterized by a form of stratification of flows which enables the definition of interfaces, the slowness of the natural dynamics which ensures that these interfaces have a smooth and stable behavior and the fact that the flows are essentially orthogonal to the wall (Dupuit's hypothesis). These points allow the vertical integration of the Richards equation in the saturated zone. In this spirit, a lot of 2D models are developed and used since the 60's (see e.g. the works of Jacob Bear, [11,12]). For more historical notes on the origin of groundwater modeling, we refer the interested readers to [17,18,21,28]. A main weakness of the approach by vertical integration lies in its justification. It is only valuable for very precise length and time scales, the time scale in particular being completely different of the typical durations of chemical reactions. However, such 2D models are widely used, even if it is particularly difficult to couple them correctly to the flow in the unsaturated part of the subsoil. Several numerical studies have been conducted in this direction. Let us mention the work of [24] where the integrated model is directly coupled with GLOBAL EXISTENCE RESULT FOR A MODEL OF SHALLOW AQUIFER IN INTERACTION WITH OVERLAND WATER. 2 a surface model. In [8], [9] and [19], the coupling of the surface and underground flows is done with a Richards equation associated with a Signorini boundary condition (for the surface behavior). A class of models is proposed in [13] which consists in coupling purely vertical models (for describing the flow at a small time scale) with an horizontal model (describing the flow at a long time scale). They admit the same behavior than the 3D-Richards model for any time scale when the aquifer present a small deepness compared to its large horizontal dimensions. They describe the essentially horizontal flow of a water table and the essentially vertical water supply flux from the surface through the unsaturated part between the groundwater and the ground level; In [32] we can find a presentation of a rather similar model coupling 1D-Richards equation with a simplified model in the saturated part. Finally, in [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF], this kind of model is integrated into a computational code called "SHE" (for "European Hydrological System" and later became SHETRAN) in the case where the water table remains away from ground level.

In this paper, we present a model belonging to the class of "Dupuit-Richards" models. Indeed the first part of the model (corresponding to the capillary fringe) consists in the 3D-Richards equation, while in the second part (corresponding to the saturated zone of the aquifer), we make a vertical average of the mass conservation law. We impose the transmission properties for the pressure and the normal fluxes at the saturation interface. This model differs from that one described in [13], already because we consider the complete equations of Richards 3D in the unsaturated part. But the main difference lies in the coupling between the two areas. In [13], the coupling is done through flux terms ensuring that model is mass conservative whereas in our case, flux terms result from the vertical average of mass conservation. The transmission property makes it possible to express them according to the speed of the depth of the interface.

The mathematical study of the model is particularly delicate because of the presence of the free boundary. Moreover, there is a constant mathematical difficulty in the structure of the system of PDEs modeling the dynamics of underground water. Indeed, in the case of a unconfined water table, we have to deal with the gradual disappearance of water in the desaturation zone and thus the disappearance of a main unknown of the problem (even in the simplified model of Richards). But the main difficulty is certainly the coupling between the two zones that express themselves by flux terms at the interface. First of all, to simplify somewhat the mathematical analysis of the problem, we will assume that the water contribution of the desaturation zone is taken into account by a variable ( in time and space) source term. But, in any case, this analysis requires to have a very regular pressure. This regularity will be deduced from assumptions on the parameters characterizing the Richards equation.

There is a huge literature regarding the classical Richards equations. Let us mention the unavoidable works of Alt et al ([4, 5]) and the papers [14,22,34] devoted to the study of the degenerate in time equation ∂ t θ(p) -∆p = 0, where θ(p) denotes the moisture content. We quote also in the one-dimensional case the work of Yin ([38]) concerning the existence of weak solution for the fully degenerate problem

∂ t θ(p) -∂ x (κ(θ(p))∂ x p) = 0,
when just assuming that θ , κ > 0. In the context of pressure-driven transport problems, that involves 3D-Richards equation coupled with a hyperbolic equation, we can mention the paper of Choquet [15] where the saturation and the mobility are strongly coupled through the pressure and that one of Amirat et al [6] with a lower coupling, that is θ = θ(x) and κ = κ(x).

PERMEABILITY TENSOR K 0

The non-linear hydraulic conductivity K is given by K = κ(P ) ρ 0 g µ K 0 . The soil transmission properties are characterized by the porosity function φ and the permeability tensor K 0 (x, z). The matrix K 0 is a 3 × 3 symmetric positive definite tensor which describes the conductivity of the saturated soil at the position (x, z) ∈ Ω. We introduce

K xx ∈ M 22 (R), K zz ∈ R * and K xz ∈ M 21 (R) such that K 0 = K xx K xz K T xz K zz .
(2.5) {K0} {K0}

HYPOTHESIS

Let us now list the assumptions on the fluid and medium characteristics but also on the flow which are meaningful in the context of our problem.

HYPOTHESIS ON THE FLUID AND ON THE MEDIUM Soil Compressibility

We neglect in the model the effects of the rock compressibility, the porosity of the medium φ do not depend on the pressure variations and it is thus assumed to be a constant.

Compressibility of the fluid First, we assume that the fluid (namely here fresh water) is weakly compressible. It means that

α P 1. (2.6) {AssComp} {AssComp}
Let us exploit this assumption. In natural conditions and especially in an aquifer, one observes small fluid mobility (defined by the ratio κ/µ). First consequence of the low compressibility of the fluid combined with the low mobility of fluid appears in the momentum equation. We perform a Taylor expansion with regard to P of the density ρ in the gravity term of the Darcy equation. Neglecting the terms weighted by α P κ/µ 1 in (2.2), we get:

q = -K ∇H , K = κ(P ) ρ 0 g µ K 0 . ( 2

.7) {newDarcy} {newDarcy}

Second consequence is ∇ρ • q 1 which leads to the following simplification in the mass conservation equation (2.3):

ρ∂ t θ + θ ∂ t ρ + ρ∇ • q = ρQ.
Neglecting in this way the variation of density in the direction of flow is sometimes considered as an extra assumption called Bear's hypothesis (cf [2]). Including (2.4), that is ∂ t ρ = ρα P ∂ t P in the latter equation, we get ρ∂ t θ + ρ θα P ∂ t P + ρ∇ • q = ρQ. After simplification by ρ > 0, we finally obtain

∂ t θ + θα P ∂ t P + ∇ • q = Q.
(2.8) {new_eq_P} {new_eq_P}

Equivalently, using the hydraulic head (2.1) and the Darcy law (2.7), (2.8) can be written

∂ t θ + S 0 ∂ t H -∇ • (K ∇H ) = Q where S 0 = ρ 0 g φ α P . ( 2

.9) {new_eq_H} {new_eq_H}

We notice that if the fluid is assumed incompressible, α P = 0, then Eq. (2.8) is the classical Richards equation in pressure formulation. An adequate definition of the volumetric moisture content θ and of the mobility function κ is the key of the model. Richards hypothesis. The Richards model is moreover based on the assumption that the air pressure in the underground equals the atmospheric pressure, thus is not an unknown of the problem. One thus assumes that the saturation and the relative conductivity of the soil are given as functions of the fluid pressure P , denoted respectively by s = s(P ) and κ = κ(P ). We introduce the saturation pressure P s which is a fixed real number. The fully-saturated part of the medium corresponds to the region {x; P (., x) > P s }, while it is partially-saturated in the There is a large choice of available models for s and κ. The most classical examples for an airwater system are the van Genuchten model [37] with no-explicit dependance on the bubbling pressure but with fitting parameters, and the Brooks and Corey model [10] where λ > 0, γ = 2 + 3λ and P s < 0. Notice that our model would easily adapt to hysteresis soil properties ( [31], [35]). The important point is that these models are such that s(P ) = 1 ⇐⇒ P ≥ P s and κ(P ) = 1 ⇐⇒ P ≥ P s .

(2.13) {sat} {sat}

In particular, the water pressure is greater than the bubbling pressure P s if and only if the soil is completely saturated . The graphs of the functions s and κ given by the Brooks-Corey model are represented in Figure 1.

HYPOTHESIS ON THE FLOW

The following assumption is introduced for upscaling the 3D problem to a 2D model in the saturated part of the domain.

Dupuit approximation (hydrostatic approach)

Dupuit assumption consists in considering that the hydraulic head is constant along each vertical direction (vertical equipotentials). It is legitimate since one actually observes quasi-horizontal displacements when the thickness of the aquifer is small compared to its width and its length and when the flow is far from sinks and wells.

GEOMETRY

The aquifer is represented by a three-dimensional domain Ω := Ω x ×(h bot , h soi l ), Ω x ⊂ R 2 , function h bot (respect. h soi l ) describing its lower (respect. upper) topography. The upper and lower surfaces are thus defined by the graph of the functions h bot = h bot (x) and h soi l = h soi l (x), x ∈ Ω x . We assume that h soil (x) > h bot (x), ∀x ∈ Ω x .

(2.14) {soilbot} {soilbot} More precisely the domain is given by:

Ω = (x, z) ∈ Ω x × R | z ∈ h bot (x), h soil (x) . ( 2

.15) {omega} {omega}

We always denote by ν the outward unit normal and e 3 is the unitary vertical vector pointing up. We decompose the boundary ∂Ω of Ω in three zones (bottom, top and vertical)

∂Ω = Γ bot Γ soil Γ ver , with
Γ bot := (x, z) ∈ Ω | z = h bot (x) , Γ soil := (x, z) ∈ Ω | z = h soil (x) , Γ ver := (x, z) ∈ Ω | x ∈ ∂Ω x .
Our model split the description of the flow in two subregions of Ω (possibly time-dependent) in each of which the flow present different behavior. We denote by h the depth of the free interface separating the freshwater layer and the unsaturated part of the aquifer. The definition of these zones is thus based on the function h = h(t , x) which is an unknown of our problem. We then introduce, for a given function h = h(t , x) such that h bot ≤ h ≤ h soil :

Ω - h (t ) := (x, z) ∈ Ω | z < h(x, t ) and Ω + h (t ) := (x, z) ∈ Ω | z > h(x, t ) , ( 2
          ∂ t θ + θα P ∂ t P + ∇ • q = Q for (t , x) ∈ (0, T ) × Ω x , q • ν = 0 for (t , x) ∈ (0, T ) × Γ soil ∪ Γ ver , P t , x, h(t , x) = P s for (t , x) ∈ (0, T ) × Ω x , P (0, x, z) = P init (x, z) for (x, z) ∈ Ω + h (0). ( 2

.18) {richards} {richards}

The effective velocity q is given by q = -K ∇( P ρ o g + z), K = κ(θ(P )) K 0 ρ 0 g µ .

• Dupuit horizontal flow in the saturated zone

Upscaling procedure

We now use the approximations introduced in 2.4 to vertically integrate equation (2.9), thus reducing the 3D problem to a 2D problem. We perform the vertical integration between depths h bot and h. Since θ(P ) = φ in the saturated zone, the vertical average (2.9) leads to

h h bot S 0 ∂ t H + ∇ • q d z = h h bot Q d z.
We denote by B f = hh bot the thickness of the saturated zone and by Q the source term representing distributed surface supply of fresh water into the free aquifer:

Q = 1 B f h h bot Q d z.
Applying Leibnitz rule to the first term in the left-hand side yields:

h h bot S 0 ∂ t H d z = S 0 ∂ ∂t h h bot H d z -S 0 H |z=h ∂ t h + S 0 H |z=h bot ∂ t h bot .
We denote by H the vertically averaged hydraulic head

H = 1 B f h h bot H d z. Because of Dupuit approximation, H (x 1 , x 2 , z) H (x 1 , x 2 ), x = (x 1 , x 2 ) ∈ Ω, z ∈ (h bot , h), we have h h bot S 0 ∂ t H d z = S 0 B f ∂ t H .
We also have

h h bot ∇ • q d z = ∇ • (B f q ) + q |z=h -• ∇(z -h) -q |z=h + bot • ∇(z -h bot ),
where ∇ = (∂ x 1 , ∂ x 2 ), q = (q x 1 , q x 2 ) and the averaged Darcy velocity q

= 1 B f h h bot q d z is given by q = - 1 B f h h bot K ∇ H d z = - 1 B f h h bot K ∇ H d z = -K ∇ H , K = 1 B f h h bot K 0 ρ 0 g µ d z,
(we remind that κ(P ) = 1 for z ∈ (h bot , h)). The averaged mass conservation law for the freshwater in the saturated zone thus finally reads

S f B f ∂ t H = ∇ • (B f K ∇ H ) + q |z=h + bot • ∇(z -h bot ) -q |z=h -• ∇(z -h) + B f Q. ( 2

.19) {dupuitH} {dupuitH}

In this equation, term B f K may be viewed as the dynamic transmissivity of freshwater layer. At this point, we have obtained an undetermined system of two pdes ((2.18)-(2.19)) with three unknowns P , H and h.

Fluxes and continuity equations across the interface

Our aim is now to include in the model the continuity and transfert properties across interface. As a consequence, we express the two flux terms appearing in (2. 19) and we reduce the number of unknowns.

• Flux across the saturation interface:

The saturation interface is characterized by F (x 1 , x 2 , z, t ) = 0 ⇔ zh(x 1 , x 2 , t ) = 0, the unit normal vector ν to the interface is thus co-linear to ∇(z -h).

If no mass transfert occurs between the two areas, the normal component of the effective velocity is continue at the interface z = h. The relation ruling continuity of the normal component of the velocity thus reads

q |z=h φ -v • ν = 0,
where v denotes the interface's velocity. It satisfies

-∂ t h + v • ∇(z -h) = 0.
Then

q |z=h + -q |z=h -• ν = 0 ⇔ q |z=h + • ∇(z -h) = q |z=h -• ∇(z -h) = φ∂ t h. ( 2

.20) {interface} {interface}

Since we assume that gravity and capillary pressure effects are neglected at the interface, the water contribution of the desaturation zone will be taken into account via the variable source term Q. This is of course an important simplification of the physical situation.

• Impermeable layer at z = h soi l Since the lower layer is impermeable, there is no flux across the boundary z = h bot :

q(h bot ) • ∇(z -h bot ) = 0. ( 2
S 0 B f ∂ t h -∇ • (B f K ∇ h) = B f Q -q |z=h + • ∇(z -h) i n Ω x , K ∇ h • ν = 0 on (0, T ) × ∂Ω x , (2.24) {dupuith} {dupuith} with B f = (h -h bot ), K = 1 B f h h bot K 0 ρ 0 g µ d z and S 0 = ρ 0 g φ α P . (2.25) {K_S0} {K_S0}
The homogeneous Neumann condition on ∂Ω x is assumed to simplify the presentation. The final model (M ) coupling 3D-Richards flow and Dupuit horizontal flow consists in system (2.18), (2.22) and (2.24), namely we have

• In Ω + h (t ) the following 3d-Richards equation holds            ∂ t θ(P ) + θα P ∂ t P + ∇ • q = Q in (0, T ) × Ω + h (t ), q • ν = 0 on (0, T ) × Γ soil ∪ Γ ver , P t , x, h(t , x) = P s in (0, T ) × Ω x , P (0, x, z) = P 0 (x, z)
in Ω + h (0). The effective velocity q is given by q = -K ∇( P

ρ o g + z), K = κ(θ(P )) K 0 ρ 0 g µ .
• In Ω - h (t ) the pressure P satisfies

P (t , x, z) = ρ 0 g P s ρ 0 g + h -z in (0, T ) × Ω - h (t ). • The depth of Γ h , h, satisfies in Ω x      S 0 B f ∂ t h -∇ • (B f K ∇ h) = B f Q -q |z=h + • ∇(z -h) in (0, T ) × Ω x , K ∇ h • ν = 0 on (0, T ) × ∂Ω x , h(0, x) = h 0 (x) in Ω x .

MASS CONSERVATION

We start by proving that the model is always mass conservative. Let M tot = M tot (t ) the total mass of the water contained in domain Ω. We denote by M + h = M + h (t ) the mass of the water filling the domain Ω + h , where the 3d-Richards problem stands and by M - h = M - h (t ) the mass of the water contained in domain Ω - h , where the Dupuit's approximation is considered. We have Thanks to the first equation of (2.18) we deduce

M + h (t ) = Ω x h soil h(t ,x) φ ρ s(P ) d z d x and M - h (t ) = Ω x h(t ,x) h bot (x) φ ρ d z d x, ( 2
Ω x h soi l (x) h(t ,x) ρ φ ∂s(P ) ∂t + φ s(P )α P ∂P ∂t d z d x = Γ + h ρ q |z=h + • ν d σ + Ω x h soi l (x) h(t ,x) ρQ d z d x, (here ν = ∇(z -h)/|∇(z -h)|). Besides, integrating (2.24) on Ω x leads to Ω x ρ S 0 B f ∂ t h d x = Ω x ρB f Q d x -ρ Ω x q |z=h -• ∇(z -h) d x, ⇔ Ω x h(t ,x) h bot (x) ρ S 0 ∂ t h d z d x = Ω x h(t ,x) h bot (x) ρQ d z d x -Γ + h ρ q |z=h + • ν d σ.
We conclude by summing up the two above equations.

MATHEMATICAL SETTING AND MAIN RESULTS

We remind that Ω is an open bounded domain of R 3 and Ω x corresponds to the projection of Ω on the horizontal plane. We denote by ν the outward unit normal pointing outward Ω. The boundary of Ω, assumed C 1 , is denoted by Γ and Γ = Γ soi l ∪ Γ bot ∪ Γ ver . The time interval of interest is (0, T ), T being any nonnegative real number, and we set Ω T = (0, T ) × Ω. For the sake of simplicity, we assume that the bottom surface of the aquifer is at constant depth, h bot ∈ R.

SOME AUXILIARY RESULTS

Let Ω an open bounded domain of R 3 . For the sake of brevity we shall write H 1 (Ω ) = W 1,2 (Ω ) and

V = H 1 0 (Ω ), V = H -1 (Ω ), H = L 2 (Ω ). The embeddings V ⊂ H = H ⊂ V are dense and compact. For any T > 0, let W (0, T, Ω ) denote the space W (0, T, Ω ) := ω ∈ L 2 (0, T ;V ), ∂ t ω ∈ L 2 (0, T ;V ) endowed with the Hilbertian norm • W (0,T,Ω ) = • 2 L 2 (0,T ;V ) + ∂ t • 2 L 2 (0,T ;V )
W (0, T, Ω ) ⊂ C ([0, T ]; [V,V ] 1 2 ) = C ([0, T ]; H ) while the embedding W (0, T, Ω ) ⊂ L 2 (0, T ; H ) (3.1)
is compact (Aubin's Lemma, see [36]).

It will also be useful to introduce the space

X (0, T, Ω ) = L ∞ (0, T ; H 1 (Ω )) ∩ H 1 (0, T ; L 2 (Ω ))
embedding with the norm u X (0,T,Ω ) = u L ∞ (0,T ;H 1 (Ω )) + u H 1 (0,T ;L 2 (Ω )) .

The following result by F. Mignot (see [23]) is used in the sequel.

Lemma 3.1. Let f : R → R be a continuous and nondecreasing function such that lim sup |λ|→+∞ f (λ)/λ < +∞. Let ω ∈ L 2 (0, T ; H ) be such that ∂ t ω ∈ L 2 (0, T ;V ) and f (ω) ∈ L 2 (0, T ;V ). Then

〈∂ t ω, f (ω)〉 V ,V = d d t Ω ω(•,y) 0 f (r ) d r d y in D (0, T ).
Hence for all 0

≤ t 1 < t 2 ≤ T t 2 t 1 〈∂ t ω, f (ω)〉 V ,V d t = Ω ω(t 2 ,y) ω(t 1 ,y) f (r ) d r d y.
The second auxiliary lemma is a parabolic extension of the Meyers regularity theorem [29]. The aim is to obtain a precise estimate of a solution of a parabolic system in X p = L p (0, T ;W 1,p 0 (Ω)), p ≥ 2, endowed with the norm

( T 0 ||v(t )|| p W 1,p 0 (Ω) d t ) 1/p =: ||∇v|| L p (Ω T ) N .
The proof may be recovered in Chapter 2, Theorem 2.2 and Chapter 1, Section 4 of [7]. This result is essentially founded on the fact that the application

v → divv from (L p (Ω T )) N into Y p = L p (0, T ;W -1,p (Ω)) is onto. The space Y p is endowed with the norm || f || Y p = inf div x g = f ||g || (L p (Ω T )) N .
Given F ∈ Y p , there is a unique solution u ∈ X p of the following initial boundary value problem

∂ t u -∆u = F in Ω T , u = 0 on (0, T ) × Γ, u(0, x) = 0 in Ω. We set Λ -1 = ∂ t -∆, so that u = Λ(F ). Let g be defined by g (p) := ||Λ|| L (Y p ;X p ) . It is well-known that g (2) = 1. Now, let A ∈ (L ∞ (Ω)) N ×N be such that there exists α > 0 satisfying N i , j =1
A i , j (x)ξ i ξ j ≥ α|ξ| 2 for a.e. x ∈ Ω and for all ξ ∈ R N .

We set β := max 1≤i , j ≤n ||A i , j || L ∞ (Ω) and

A u = - N i , j =1 ∂ x i A i , j ∂ x j u .
We state the following Lemma (cf [7]).

Lemma 3.2. Assume that

A is symmetric. Let f ∈ L 2 (0, T ;V ), u 0 ∈ H and u ∈ L 2 (0, T ;V ) be the solution of ∂ t u + A u = f in Ω T , u(0) = u 0 in Ω.
There exists r > 2, depending on α, β and Ω, such that if u 0 ∈ W 1,r 0 (Ω) and f ∈ L r (0, T ;W -1,r (Ω)), then u ∈ L r (0, T ;W 1,r 0 (Ω)). Furthermore, the following estimate holds true

||u|| L r (0,T,W 1,r 0 (Ω)) ≤ C (α, β, r )(|| f || L r (0,T ;W -1,r (Ω)) + ||u 0 || W 1,r 0 (Ω) ), (3.2)
where the constant C (α, β, r ) > 0 depends on Ω, α, β and r (but not on T ) as follows: Then, the smaller (1µ), the larger r .

C (α, β, r ) ≤ g (r ) (1 -k(r )) β , k(r ) = g (r )(1 -µ) ( 3

MAIN RESULTS

We aim giving an existence result of physically admissible weak solutions for model (M ) completed by initial and boundary conditions.

We introduce function T l defined by

T l (u) = (u -h bot ) ∀u ∈ (h bot , h soi l ),
which is extended continuously and constantly outside (h bot , h soi l ). Function T l (h) represents the thickness of the saturated freshwater zone in the reservoir while l = h soi l -h bot corresponds to the thickness of the subsoil. We also emphasize that the function T l also acts on the source term Q for avoiding the pumping when the thickness of freshwater zone is smaller than 0.

Setting Ω + h = Ω x × (h, h soi l ) and Ω - h = Ω x × (h bot , h), we thus consider the following system : As we shall see later, it is possible to establish an existence result for h in the space W (0, T, Ω x ) in the degenerate case, but if we want to extend this result to the space X (0, T, Ω x ), we must introduce a regularization of the diffusive term. So, let δ > 0, we introduce the regularization T δ = T l + δ and the following regularized equation

∂ t θ(P ) + θ(P )α P ∂ t P + ∇ • q = Q, q = -K (θ(P ))∇( P ρ o g + z), in (0, T ) × Ω + h , (3.5) {equation1} {equation1} P (t , x, z) = ρ 0 g P s ρ 0 g + h -z in (0, T ) × Ω - h , (3.6) {equation2} {equation2} φ∂ t h -∇ • T l (h) K ∇ h = T l (h) Q in (0, T ) × Ω x . ( 3
φ∂ t h δ -∇ • T δ (h δ ) K ∇ h δ = T l (h δ ) Q in (0, T ) × Ω x . ( 3
P | Γ h = P s in (0, T ), ∇P • ν = 0 on (0, T ) × (Γ soi l ∪ Γ ver ), P (0, x, z) = P 0 (x, z) in Ω + h 0 . (3.9) ∇h • ν = 0 on (0, T ) × ∂Ω x , h(0, x) = h 0 (x) in Ω T . (3.10)
Function P s is constant with respect to the time and the space. Function P 0 ∈ H 2 (Ω) satisfies the compatibility condition P 0 (x, h 0 ) = P s in Ω + h 0 . We also assume that h

0 ∈ L ∞ (Ω x ) is such that h 0 ≥ h bot a.e. in Ω x . Source term Q is given function of L 2 (0, T ; H ).
Let us now detail the mathematical assumptions. We begin with the characteristics of the porous structure. We limit our study to the isotropic case so K 0 is assumed to be a scalar. In the saturated part, the averaged hydraulic conductivity K is thus equal to the constant K 0 ρ 0 g µ . From now, the density ρ 0 will be denoted by ρ. Functions θ and κ are pressure-dependent and we assume

θ ∈ C 1 (R), 0 ≤ θ(x) ≤ θ + , θ (x) ≥ 0 ∀x ∈ R, (3.11) {theta} {theta} κ ∈ C (R), 0 ≤ κ(x) ≤ κ + ∀x ∈ R + . ( 3

.12) {kappa} {kappa}

For the previous parabolic system, we state and prove the following existence result.

Theorem 3.5. Assume that there exist two real numbers θ -and κ -such that Remark 3.6. As mentioned above, the equation (3.7) becomes degenerate when the free interface touch the bottom of the aquifer. For the first step of the proof, it is possible to overcome this degenerescence thanks to an entropy functional. But to prove more regularity for h, we have to assume that the thickness of fresh water inside the aquifer is always positive i.e. (hh bot ) ≥ δ > 0. This provides an interpretation of the diffusion coefficient δ. Another interpretation of δ is to see it as the thickness of the interface between the saturated zone and the unsaturated zone (see [16]) which thus would no longer be assumed sharp.

θ(x) ≥ θ -> 0 ∀x ∈ R, κ(x) ≥ κ -> 0 ∀x ∈ R + . ( 3
Functions θ and κ characterize the mathematical type of the system. More precisely, problem (3.5)-(3.7) is of parabolic type if θ and κ are positive functions and of degenerate parabolic type if θ and κ are non-negative functions. The second result of this article is devoted to the degenerate parabolic framework of the system. In this case, we prove the following result Theorem 3.7. Assume that functions θ and κ are non-negative and verify (3.11)-(3.12). Assume, moreover, there exists ε 0 > 0 such that κ is increasing in (0, ε 0 ). Next section is devoted to the proof of Theorem 3.5.

PROOF OF THEOREM 3.5

Let us sketch our strategy. The problem is characterized by the presence of a free interface between the two domains, by the difficulties inherent in the Richards equations and by the coupling between the two equations. Furthermore we must face the gradual disappearance of water in the desaturation zone and thus the disappearance of a main unknown of the problem. Since the system is strongly coupled, we apply a fixed point approach to solve it. The key is to first solve equation in h. It is possible thanks to the continuity of the normal component of the Darcy flux through the interface which allows to express the flux with respect to the time derivative of the interface depth h. We need to get sufficiently regularity result for h and its time derivative in order to "linearized' (in some sens) the pressure equation. This regularity can be obtained thanks to a regularization which enables to ensure a thickness of freshwater zone always greater than δ > 0 inside the aquifer. To overcome the difficulty related to the strong nonlinearities in pressure equation, we perform a change of variable in pressure equation. We thus use the fundamental Kirchoff's transform to linearize the divergence part of Eq. (3.5) on the variable domain depending on h computed at the previous step. We finally establish sufficient uniform estimates for the pressure on the whole domain. By using a Schauder fixed point theorem, we prove an existence result for the full problem.

Without lost of generality, we can simplify the equations by taking null source term Q = 0 in the desaturation zone for the existence proof.

We claim the following result. Let us sketch our strategy. First step consists in using a Schauder fixed point theorem for proving an existence result for an auxiliary regularized problem. More precisely we regularize the function T l with the parameter δ > 0. Subsequent difficulty is that the mapping used for the fixed point approach has to be continuous in L 2 (0, T ; H 1 (Ω x )). We show that the regularized solution satisfies the maximum principle announced in Lemma 1. We finally establish sufficient uniform estimate (thanks to an entropy functional) to let the regularization δ tends to zero.

Proof.

• Step 1 : Existence for the regularized system : Let δ > 0 we introduce the regularization T δ = T l + δ and the regularized equation

φ∂ t h δ -∇ • T δ (h δ ) K ∇ h δ = T l (h δ ) Q in (0, T ) × Ω x . Let h ∈ W (0, T, Ω x ), we introduce the linear problem : Find h ∈ W (0, T, Ω x ) verifying φ∂ t h δ -∇ • T δ ( h) K ∇ h δ = T l ( h) Q in (0, T ) × Ω x , (4.1)
∇h δ • ν = 0 on (0, T ) × ∂Ω x , h δ (0, x) = h 0 in Ω x .
From now, we omit the subscript δ in h δ . It follows from ( [25,27]) that for every h ∈ W (0, T, Ω x ), there exists a unique solution h := F ( h) ∈ W (0, T, Ω x ) of problem (4.1).

Sequential continuity of F in L 2 (0, T ; H ) when restricted to any bounded subset of W (0, T, Ω x ).

Assume given a bounded sequence ( hn ) in W (0, T, Ω x ) and a function h ∈ W (0, T, Ω x ) such that

h n → h in L 2 (0, T ; H ).
We thus have hn h weakly in W (0, T, Ω x ); that is, hn h weakly in L 2 (0, T,V ) and d hn d t d h d t weakly in L 2 (0, T,V ). Set h n = F ( hn ) and h = F ( h). We intend to show that h n → h weakly in W (0, T, Ω x )and thus strongly in L 2 (0, T ; H ) thanks to a classical result of Aubin.

Pick a constant M > 0, that we will precise later on, such that

||∇ h|| L 2 (0,T ;H ) ≤ M . ( 4

.2) {estiM} {estiM}

For all n ∈ N, h n satisfies (4.1). Pick any τ ∈ [0, T ] and take w = h n χ (0,τ) (t ) in (4.1). It leads to

φ τ 0 〈∂ t h n , h n 〉 V ,V + Ω x,τ K T δ ( hn )∇h n • ∇h n = Ω x,τ T l ( hn ) Q h n d xd t .
The functions h n belong to W (0, T ) and hence to C ([0, T ]; L 2 (Ω)). Thanks to Lemma 3.1, we can write

τ 0 〈∂ t h n , h n 〉 V ,V d t = 1 2 ||h n (•, τ)|| 2 H - 1 2 ||h 0 (•, 0)|| 2 H .

GLOBAL EXISTENCE RESULT FOR A MODEL OF SHALLOW AQUIFER IN INTERACTION WITH OVERLAND WATER. 14

On the other hand, we have that

Ω x,τ K T s ( hn )∇h n • ∇h n d xd t ≥ K δ||∇h n || 2 L 2 (0,τ;H ) .
Using Cauchy-Schwarz and Young inequalities, we obtain

Ω x,τ T δ ( hn ) Q h n d xd t ≤ Q 2 L 2 (0,T ;H ) φ l 2 + φ 4 ||h n || 2 L 2 (0,T ;H ) .
Using the above estimates, we obtain for all τ ∈ [0, T ]

φ 4 ||h n (•, τ)|| 2 H + K δ||∇h n || 2 L 2 (0,τ;H ) ≤ φ 2 ||h 0 (•, 0)|| 2 H + Q 2 L 2 (0,T ;H ) φ l 2
We infer from the above inequality that there exist real numbers A M = A M ( K , h 0 , l , φ) and B M = B M ( K , h 0 , l , φ) depending only on the data of the problem such that

||h n || L ∞ (0,T ;H ) ≤ A M , ||h n || L 2 (0,T ;V ) ≤ B M . (4.3) {pbh_reg} {pbh_reg} Thus the sequence (h n ) n is uniformly bounded in L ∞ (0, T ; H ) ∩ L 2 (0, T ; H 1 (Ω x )). Set C M = max(A M , B M ), V 1 = H 1 (Ω x ) and V 1 = (H 1 (Ω x )) .
We now prove that (∂ t h n ) n is bounded in L 2 (0, T ;V 1 ). We have

||∂ t h n || L 2 (0,T ;V 1 ) = sup ||w|| L 2 (0,T ;V 1 ) ≤1 T 0 〈∂ t h n , w〉 V 1 ,V 1 = sup ||w|| L 2 (0,T ;V 1 ) ≤1 T 0 - 1 φ Ω T K T δ ( hn )∇h n • ∇w - Ω T Q T l ( hn ) w .
Besides (since δ is expected to tend to 0, hence we may assume that δ < 1)

Ω x,T
K T δ ( hn )∇h n • ∇w ≤ K (l + 1) ||h n || L 2 (0,T ;V 1 ) ||w|| L 2 (0,T ;V 1 ) , and since h n is uniformly bounded in L 2 (0, T ; H 1 (Ω x )), we write

Ω x,T K T δ ( hn ) ∇h n • ∇w d xd t ≤ K l C M ||w|| L 2 (0,T ;V 1 ) .
Furthermore we have

Ω T
QT l ( hn )w d xd t ≤ Q L 2 (0,T ;H ) l ||w|| L 2 (0,T ;V 1 ) .

Summing up these estimates, we conclude that

||∂ t h n || L 2 (0,T ;V 1 ) ≤ M := l φ ( K C M + Q L 2 (0,T ;H ) ).
We have proved that the sequence h n n is uniformly bounded in the space W (0, T, Ω x ). Using Aubin-Lions' lemma, we can extract a subsequence (h n k ) k , converging strongly in L 2 (Ω x,T ), almost everywhere in (0, T ) × Ω x , and weakly in W (0, T, Ω x ) to some limit denoted by v. From the convergence a.e. in Ω x,T , we see that for all w ∈ W (0, T, Ω x ), T δ ( hn )∇w → T δ ( h)∇w strongly in L 2 (Ω T ) by dominated convergence. It follows that v solves (4.1) and (3.10). By uniqueness of the solution of that system, we conclude that v = h and that the whole sequence h n → h weakly in W (0, T, Ω x ) and strongly in L 2 (0, T ; H ). The sequential continuity of F in L 2 (0, T ; H ) is established.

Existence of C ⊂ W (0, T, Ω x ) such that F (C ) ⊂ C .

We aim now to prove that there exists a nonempty bounded closed convex set of W (0, T, Ω x ) × L 2 (0, T ; (H 1 (Ω x )), C , such that F (C ) ⊂ C . We notice that this result will imply that there exists a real number M > 0, depending on initial data, such that for h = F ( h) ∈ W , we have ||∇h|| L 2 (0,T ;H ) ≤ M .

Taking w = h ∈ L 2 (0, T ;V 1 ) in (4.1) yields

φ T 0 〈∂ t h, h〉 V 1 ,V 1 d t + Ω x,T K T δ ( h)∇h • ∇h d xd t = Ω x,T QT l ( h) h d xd t .
We apply Lemma 3.1 to the function f defined by f (u) = u for u ∈ R in order to compute the first terms of previous equality. We get for all τ ∈ (0, T ) Applying Gronwall's lemma, we deduce that there exists a constant M := M (φ, l , δ, K , h 0 , Q) > 0 such that h L 2 (0,T ;V 1 ) ≤ M .

τ 0 〈∂ t h, h〉 V 1 ,V 1 d t = 1 2 Ω x h 2 (
Furthermore we deduce the estimate in L 2 (0, T,V 1 ) of ∂ t h as previously. Introduce the set C := {h ∈ W (0, T, Ω x ); h(0, . = h 0 , h L 2 (0,T ;V 1 ) ≤ M , ∂ t h L 2 (0,T,V 1 ) ≤ M }. (4.4)

Then C is a nonempty, closed, convex, bounded set in L 2 (0, T ; H ), defined such that F (C ) ⊂ C . Since C is also a bounded set in W (0, T, Ω x ), we also proved that F is sequentially continuous in L 2 (0, T ; H ). For the fixed point strategy, it remains to show the compactness of F (C ). Since we work in metric spaces, proving its sequential compactness is sufficient. The compactness of F (C ) is straightforward due to the Aubin's theorem. We now have the tools for using the Schauder's fixed point theorem [ We set h m = hh bot -= inf(0, hh bot ) ∈ L 2 (0, T ;V ).

Let τ ∈ (0, T ). We recall that h m (0, •) = 0 a.e. in Ω thanks to the maximum principle satisfied by the initial data h 0 . We now claim that there exists a constant C > 0 such that ∀δ > 0, h δ L 2 (0,T ;V 1 ) ≤ C .

We fix δ > 0 and we write h instead of h δ . We define the nonnegative entropy functional as it is done in [3] S (a) -1 =    a ln aa for a > 0, 0 for a = 0, +∞ for a < 0. Formally, the variation of the entropy functional

Ω x
S ((h-h bot )d x may be evaluated by letting w = S (hh bot ) in (3.7). However, S (hh bot ) is unbounded as h → h bot , so we are led to substitute S ((h -h bot )+δ) to w = S (h -h bot ) (see [33]). Letting w = S ((h -h bot )+δ) in (3.8) leads to We now proceed to the last step in the proof of Lemma 4.1.

φ T 0 〈∂ t h, S ((h -h bot ) + δ)〉 V 1 ,V
• Step 4 : δ → 0

We infer from (4.7) that h δ is uniformly bounded in W (0, T, Ω x ). Extracting a subsequence if needed, we may assume that for some h ∈ W (0, T, Ω x ) h δ → h in W (0, T, Ω x ) as δ → 0.

The maximum principle readily follows from (4.5). Letting δ → 0. in (3.8), we get at once (3.7) using the convergence a. e. in Ω x,T (and thus T δ (h δ ) → T l (h) a. e. in Ω x,T ). Finally, Eq. (3.10) holds true since the map h ∈ W (0, T, Ω x ) → h(t = 0, .) ∈ H 1 (Ω x ) is continuous. The proof of Lemma 4.1 is complete.

In order to process the coupling with the unsaturated zone, we need more regularity for the interface depth h more specially for the partial time derivative of h. We can not obtain this regularity in the case where the saturation interface h touches the bottom of the aquifer, which corresponds to the previous degenerate case. We thus consider the regularized problem (3. Proof. Let δ > 0 and h ∈ W (0, T, Ω x ) a solution of (3.8)-(3.10). We multiply Eq. (3.8) Assume that h 0 ∈ W 1,r (Ω x ) and Q ∈ L r (0, T ;W -1,r (Ω)). Then ∇h belongs to (L r (Ω x,T )) 2 , and is bounded as follows, uniformly with regard to T : δ + < 1 which could be very restrictive.

||∇h|| (L r (Ω x,T )) 2 ≤ C r (l , δ, K , h 0 ). ( 4
Proof. We adapt the proof of Lemma 4.1. We turn back to the construction of the intermediate solution which appears as the fixed point of an application in Step 1 of the proof of Theorem 3.5. We recall its outline. If F is the application defined in (4.1) and if C is the nonempty (strongly) closed convex bounded subset of the space (L 2 (0, T ; H )) 2 defined in (4.4), we have shown that F (C ) ⊂ C and that there exists h ∈ C such that F (h) = h. This fixed point for F is a weak solution of problem (3.8)-(3.10) in (L 2 (0, T ; H 1 (Ω x )) 2 . Now, we prove that, if the assumptions of Proposition 4.2 are fulfilled, this solution is actually in L r (Ω x,T ), r > 2. To this aim, we modify the definition of the convex bounded subset C by including an estimate in the norm L r (0, T ;W 1,r (Ω x )) of its elements.

Let M be a strictly positive real number that we will define later on. We set D := {h ∈ L r (0, T ;W 1,r (Ω x )), h(0) = h 0 , h L 2 (0,T ;V ) ≤ M , ∇h (L r (Ω x,T )) N ≤ M }. We emphasize that the real M does not depend on the real number M . We have the tools to perform a fixed point analysis similar to the one in the proof of Theorem 3.5. We have already chosen M so that the bounded convex D defined by (4.13) satisfies F (D) ⊂ D. Let us show that D is closed in L 2 (0, T ; H ). In fact, it is sufficient to check that, if h n denotes a sequence of function of D such that h n → h in L 2 (0, T ; H ), then ∇h ∈ (L r (Ω x,T )) 2 with ∇h (L r (Ω x,T )) N ≤ M . Due to the definition of D, the sequence (∇h n ) n is uniformly bounded in the space (L r (Ω x,T )) 2 . Thus, there exits v ∈ (L r (Ω x,T )) 2 such that, for an appropriate subsequence here characterized by an increasing function ϕ, the convergence ∇h ϕ(n) v holds true weakly in (L r (Ω x,T )) 2 . It means Since r > 2, we have L 2 (Ω x,T ) ⊂ L r (Ω x,T ) and then we infer from the latter convergence together with (4.15) that ∇h = v in L r (Ω x,T ). We conclude the proof thanks to (4.14). In brief, D is a nonempty convex, bounded closed set in (L 2 (0, T ; H )), satisfying F (D) ⊂ D.

Ω x,T ∇h ϕ(n) • Φ d xd t → Ω x,T v • Φ d xd t , ∀Φ ∈ (L r (Ω x,T )) 2 , 1 r + 1 r = 1, ( 4
The remainder of the proof follow the lines of the one of Lemma 4.1. It follows from Schauder fixed point theorem that there exist h ∈ D such that F ( h) = h. This fixed point is a weak solution of problem (3.8)-(3.10) and its gradient is uniformly bounded in the space (L r (Ω x,T )) 2 . The proof of Proposition 4.3 is complete. Remark 4.5. If we take into account the contribution of water coming from upper zone of the aquifer, the source term can be expressed by the flux q |z=h + .∇(zh) where q |z=h + is the Darcy flux in the capillary fringe. We thus have to estimate the L 2 -norm of this flux, which represents the main difficulty of the mathematical analysis. It implies the L 4 -norm of q (and thus the H 2norm of the pressure) which could be estimated by applying Lemma 3.2 to a linearization of Richards equation. We now construct the framework to apply the Schauder fixed point theorem (see [20,39]). For the fixed point strategy, we introduce the convex subset K p of W (0, T, Ω) . We set

q |z=h + • ∇(z -h) = q |z=h -• ∇(z -h) = -K ∇h • ∇(z -h) = K ∇h • ∇h, (4.23)
= K |∇(z -h)| 2 -K = K ρ 2 g 2 |∇P (z = h -)| 2 -K .
We state the following result. Proof. Sequential continuity of F in L 2 (0, T ; H ) when restricted to any bounded subset of W (0, T Ω).

Assume given a bounded sequence pn in (W (0, T, Ω) and a function p ∈ W (0, T Ω) such that pn → p in L 2 (0, T ; H ). We thus have pn p weakly in W (0, T, Ω); that is, pn p weakly in L 2 (0, T,V ) ( and d pn d t d p d t weakly in L 2 (0, T,V )). Set p n = F ( pn ) and p = F ( p). We intend to show that p n → p weakly in W (0, T, Ω)and thus strongly in L 2 (0, T ; H ) thanks to a classical result of Aubin. Lemma 4.7 implies that the sequence p n n is uniformly bounded in the space W (0, T, Ω). Using Aubin-Lions' lemma, we can extract a subsequence (p n k ) k , converging strongly in L 2 (Ω T ), almost everywhere in (0, T ) × Ω, and weakly in W (0, T, Ω) to some limit denoted by v. From the convergence a.e. in Ω T , we see that for all w ∈ W (0, T, Ω), κoθoB -1 ( pn )w → κoθoB -1 ( p)w (and also τ( pn )w → τ( pn )w ) strongly in L 2 (Ω T ) by dominated convergence. It follows that v solves (4.21)- (4.22). By uniqueness of the solution of that system, we conclude that v = p and that the whole sequence p n → p weakly in W (0, T, Ω) and strongly in L 2 (0, T ; H ). The sequential continuity of F in L 2 (0, T ; H ) is established. Existence of C ⊂ W (0, T, Ω) such that F (C ) ⊂ C .

We aim now to prove that there exists a nonempty bounded closed convex set of W (0, T, Ω), C , such that F (C ) ⊂ C . We notice that this result will imply that there exists a real number K (M ) > 0, depending on initial data, such that for p = F ( p) ∈ W (0, T, Ω), we have ||∇p|| L 2 (0,T ;H ) ≤ K (M ).

We take the set C = K p and the constant K (M ) = T M p . Then C is a nonempty, closed, convex, bounded set in L 2 (0, T ; H ), defined such that F (C ) ⊂ C (thanks to Lemma 4.7). Since C is also a bounded set in W (0, T, Ω), we also proved that F is sequentially continuous in L 2 (0, T ; H ). For the fixed point strategy, it remains to show the compactness of F (C ). Since we work in metric spaces, proving its sequential compactness is sufficient. The compactness of F (C ) is straightforward due to the Aubin's theorem.

GLOBAL EXISTENCE RESULT FOR A MODEL OF SHALLOW AQUIFER IN INTERACTION WITH OVERLAND WATER. 22

We now have the tools for using the Schauder's fixed point theorem [39,Corollary 9.7]. There exists p ∈ C such that F (p) = p. Then p is a weak solution of problem of (4.21)- (4.22). We end the proof by considering the inverse Kirchoff's transform to turn back to the original problem (3.5)-(3.7), (3.10)-(3.9). The proof of Theorem (3.5) is achieved.

PROOF OF THEOREM 3.7

We now aim to establish the existence result for the degenerate problem of Theorem 3.7. Degenerate means that we no longer assume a residual positive saturation in the desaturation zone. So saturation may be zero in some areas of the aquifer. From Theorem 3.5, we can assert that there exists a weak solution (p ε , h ε ) ∈ Z (0, T ) 2 of the following parabolic problem, for any ε > 0 ∂ t θ ε (P ε ) + θ ε (P ε )α P ∂ t P ε + ∇ • q ε = 0, q ε = -κ ε (θ ε (P ε ))∇( P ε ρg + z), in (0, T ) × Ω + h , (5.1) {eq1} {eq1} P ε (t , x, z) = ρ 0 g P s By following the original idea presented in [15], we prove that there exists an extracted subsequence of solution of regularized problem (5.1)-(5.5) that weakly converges to a solution of original problem. We first recall some uniform estimates. Due to (3.14), the conductivity κ is an increasing function in (0, ε) as soon as ε < ε 0 . The estimates stated for p ε in section 4 written for P ε become θ(P ε ) + θ (P ε ) + ε + κ θ ε (P ε ) + ε ∂ t P ε L 2 2(Ω T ) ≤ C (5.7)

ρ 0 g + h ε -z in (0, T ) × Ω - h , (5.2) {eq2} {eq2} φ∂ t h ε -∇ • (T l (h ε ) + δ) K ∇ h ε = T l (h ε ) Q in (0, T ) × Ω x . ( 5
κ θ ε (P ε ) + ε ∇P ε (L ∞ (0,T ;L 2 2(Ω))) 3 ≤ C .

(5.8)

These estimates are completely useless in potentially zero saturation areas. Indeed functions θ, θ and κoθ are null in (-∞, P d ) so we cannot ensure that P ε (x, t ) > P d . The idea in [15] is to introduce a convenient truncature function. First let H be a primitive of the function θ(κoθ) H (q) = q θ(s)(κoθ)(s) d s.

(5.9) {primi_H} {primi_H} From (5.7)-(5.8), function H (P ε ) is uniformly bounded in H 1 (Ω T ). We then define the limit function H such that

H (P ε ) → H in L 2 (

.2. 6 )

 6 Équation d'état pour la compréssibilité du fluide On considère que le fluide est compressible ainsi la pression P est liée à la masse volumique ρ par la relation d ρ ρ = α P d P , c'est à dire :

1 FIGURE 1 . 2 -

 112 FIGURE 1.2 -Évolution de la saturation du sol obtenue par le modèle de Richards 3d : cas test 1.

1 FIGURE 1 . 3 -

 113 FIGURE 1.3 -Évolution de la saturation obtenue par le modèle de Richards 3d : cas test 2.

10 zFIGURE 1 . 4 -

 1014 FIGURE 1.4 -Évolution de la surface de la nappe phréatique h 2d sat et h R sat pour R ∈ {0, 7, ∞}. Les quatre premiers tracés correspondant au cas test 1 et les quatre suivants au cas test 2.

.5. 11 )

 11 Alors le système (1.5.3)-(1.5.7), (1.5.4), (1.5.6)-(1.5.8) admet une solution faible (P, h) satisfaisant (a) la fonction P

Figure 1 :Figure 2 :

 12 Figure 1: Bidimensional representation of the cylindrical geometry of the problem: Ω x ⊂ R is an interval.

  20) This setting corresponds to the simplest form of the model (3.5)-(3.9) since (3.20) is a classical boundary value problem. Nevertheless the simulations below illustrate that it is not the better form of approximation for the 3D-Richards equation.Velocity of the flow. The velocity v of the flow turns out to be the superposition of the two velocities u and w which respectively describe the fast and slow components of the flow. Actually u (resp. w) is the dominant component of the flow in the short time scale (resp. large time scale).

3. 3 . 12 )

 312 Comments on the model in the cases (3.11) and (3.Now we come back to the model (3.5)-(3.9) in which we set the virtual interface h by

  soil level h soil . In this sense Ω - h sat does not represent the physical water table which possibly touches the soil level. We only have done this choice for the definition of h sat to recover the unknown h in the maximal case(3.11) and thus to facilitate the visualisation. Numerical scheme. For the numerical approximation of the problem (3.5)-(3.9) we use mass-conservative fully implicit time schemes associated with finite elements methods in space for both horizontal and vertical directions. The schemes for (3.10) and (3.12) differ: • In the case (3.10), we solve directly equation (3.7) in which the right-hand side (u •e 3 )| Γ + h is seen as a Dirichlet to Neumann operator depending on H and obtained by solving the 1D-vertical Richards equations. This non-linear term is treated with a Newton method.

  sat the ones coming from the model (3.5)-(3.9) with the function h given respectively by (3.10), (3.11) and (3.12).

1 Figure 3 :Figure 4 :

 134 Figure 3: Solution of the classical Richards problem in the reference test case.

Figure 5 :

 5 Figure 5: Cumulative error in space and time h 2d sat -h κ sat L 1 ([0,T ]×Ω) versus the ratio length/deepness= 1/ε of the aquifer (κ ∈ {a, b,c}). Function h 2dsat is the iso-pressure P = P s in the original 2d-Richards problem and h κ sat is the one associated with the coupled problems for three different choices of h satisfying (3.8). The test case is the one of figure3.

•

  The saturation and relative conductivity satisfy s(P ) = s(P ), k r (P ) = k r (P ). (4.1)

  Richards model. By introducing the same parameter ε 1, the rescaled coupled problem (3.5)-(3.9) reads:

) 19 70 CHAPITRE 2 .

 192 DÉRIVATION DU MODÈLE COUPLÉand for the flux at the soil level

. 19 )

 19 Moreover, since s and k r are C ∞ by part functions, we write 

26 ) 4 . 1 .

 2641 Proposition Let (P γ ε , v γ ε ) be the solution of the rescaled 3D-Richards problem (4.9)-(4.10) or of the rescaled coupled model (4.12)-(4.16) for T = ε -γ T and γ ∈ {0, 1, 2}. Assume that (4.17)-(4.20) hold true, then • satisfies (4.22) and (4.23) if α = 0, or (4.22) and (4.24) with the compatibility condition F 0 = 0 if α = 0. (4.22) and (4.23) if α = 0, or (4.22) and (4.25) with the compatibility condition F 0

(4. 30 )

 30 Short time case. We prove the first claim of Proposition 4.1 which is associated with the short characteristic time scale T = ε -γ T for γ = 0. The equation (4.9) here reads φ ∂s(P ) ∂t + ε div x (v) + ∂v • e 3 ∂z = 0. (4.31) Some computations show that the main order terms in the latter equation give φ ∂s(P 0 ) ∂t + ∂v 0 • e 3 ∂z = 0 in ]0, T [×Ω. (4.32)

47 )

 47 Taking into account(4.22),(4.46) and the fact that α = F 0 = 0, the third equation of (4.28) and the second equations of (4.29) for k = 2 become

( 4 .

 4 55) By (4.53) for k = 1, ∂ z P 1 = 0 on Ω - h 0 (t ). Then by (4.50) and the first equation of (4.53)

  )-(4.20) are substituted in the previous equations: for t ∈]0, T [ and (x, z) ∈ Ω + • e 3 )| Γ h 0 = 0 for (t , x) ∈]0, T [×Ω x . (4.59)

  intermediate time case, we substitute asymptotics (4.17)-(4.20) in the previous equations. Identifying the coefficients associated with ε n for n ∈ {1, 2}, we get ∂ z (u n • e 3 ) = 0 in Ω + h (t ) and u n • e 3 = 0 on Γ soil . This leads to u 0 • e 3 = u 1 • e 3 = 0 on Ω same arguments we obtain P 0 = ρ g (H 0 -z) and v 0 = 0 in whole Ω. System (4.22) is satisfied. The characterization of H 0 depends on the values of α. Similar arguments to those employed in the intermediate time case when α = 0 lead to (4.23).

  A CLASS OF MODELS COUPLING SLOW AND FAST FLOW IN SHALLOW AQUIFERS: CONSERVATIVE SCHEME AND VALIDITY EXPERIMENTSCHRISTOPHE BOUREL, MUNKHGEREL TSEGMID ABSTRACT. In this work we study numerically a class of models describing the flow in shallow aquifers. Those models approximate the classical 3d-Richards one by coupling the two dominant behaviors of the flow that hold in shallow geometry. The first one corresponds to a fast 1d-flow in the vertical direction. The second one is a slow 2d-horizontal one associated to an hydrostatic pressure profile (at the bottom of the aquifer). This model has been introduced in[2].

FIGURE 1 .

 1 FIGURE 1. Bidimensional representation of the cylindrical geometry of the problem: Ω x ⊂ R is an interval.

Proposition 1 . 1 .

 11 Let P , u, h and H solution of problem (1.14)-(1.18). It holds in ]0, T [×Ω x :

1 FIGURE 2 .

 12 FIGURE 2. Solution of the classical 2d-Richards problem in the first reference test case.

40 FIGURE 3 .

 403 FIGURE 3. On the left: Evolution of the ratio |t r 2d -t r c |/t r 2d with respect to the thickness of the aquifer |L x /h bot | (in per cent). On the right: saturation profile obtained with the coupled model at the time where the infiltrating water reached the water table

FIGURE 4 .

 4 FIGURE 4. Evolution of the iso-pressure P = P s obtained from the classical Richards equation (h 2d sat ) and from the coupled model for three choices of h characterized by R > 0 (h R sat for R ∈ {0, 7, ∞} respectively). The test case is the one of Figure2.

FIGURE 5 .

 5 FIGURE 5. Cumulative error in space and time1 L x h 2d sat -h R sat L 1 ([0,T ]×Ω) versus the ratio length/deepness of the aquifer (R ∈ {0, 4, 8, ∞}). Function h 2d sat is the iso-pressure P = P s in the original 2d-Richards problem and h κ sat is the one associated with the coupled problems for different choices of R characterizing h by (1.13). The test case is the one of Figure2with h 0 sat = -2.2.

1 FIGURE 6 .

 16 FIGURE 6. Solution of the 3d-Richards problem in the reference case of infiltrating overland water.

40 FIGURE 7 .

 407 FIGURE 7. On the left: Evolution of the ratio |T r 2d -T r c |/T r 2d with respect to the thickness of the aquifer |L x /h bot |. On the right:

  FIGURE 8. Evolution of the iso-pressure P = P s obtained from the classical Richards equation (h 2d sat ) and from the coupled model for three choices of h given by R = 0, R = 3 and h = h bot (h κ sat for κ ∈ {a, b, c} respectively). The test case is the one of Figure6.

  FOR A MODEL OF SHALLOW AQUIFER IN INTERACTION WITH OVERLAND WATER. CHRISTOPHE BOUREL, CAROLE ROSIER, MUNKHGEREL TSEGMID ABSTRACT. In this work, we analyze a model which describes the water flow in shallow aquifers. It is an alternative to the 3D-Richards model which is classically used in this kind of porous media. The model couples the two dominant kind of flows holding in the aquifer. The first one is described by the classical 3d-Richards problem in the upper capillary fringe. The second results from Dupuit approximation by vertically integrating the conservation laws between the bottom of the aquifer and the saturation interface.

Remark 3 . 3 .

 33 .3) {estim_c} {estim_c} where µ = α/β. In view of (3.3), the value of r depends on the characteristics (α, β) of the elliptic operator A , roughly on the ratio α/β. Actually the real number r may be chosen in the range 2 < r ≤ max{r 0 ∈ R; k(r 0 ) < 1}.(3.4) {Choixr} {Choixr}

Remark 3 . 4 .

 34 Taking into account the relation (2.20), we replaced the expression of the flux q |z=h + • ∇(z -h) with that of the time derivative of h in (3.7). Besides, since the compressibility coefficient is very small, we also neglected the term with the storage coefficient in Eq. (3.7).

  .8) {equation3bis} {equation3bis} System (3.5)-(3.7) is completed by the following boundary and initial conditions:

  .13) {theta-kappa} {theta-kappa} Then system (3.5)-(3.9), (3.6), (3.8)-(3.10) admits a weak solution (P, h) satisfying (a) the function P ∈ L ∞ (0, T ; H 1 (Ω)) ∩ H 1 (0, T ; L 2 (Ω)) is solution of (3.5)-(3.9) and (3.6); (b) the function h ∈ L ∞ (0, T ; H 1 (Ω)) ∩ H 1 (0, T ; L 2 (Ω)) is solution of (3.8)-(3.10).

( 3 .

 3 14) {kappa_reg} {kappa_reg} Then system (3.5)-(3.9), (3.6), (3.8)-(3.10) admits a weak solution (P, h) satisfying (a) the function P ∈ W (0, T, Ω) ∩ L 2 (Ω T ) is solution of (3.5)-(3.9) in L 2 (0, T ; H -1 (Ω)) and the Darcy velocity q ∈ L 2 (Ω T ) 3 ; (b) the function h ∈ L 2 (0, T ; (H 1 (Ω)) ) ∩ L 2 (0, T ; H 1 (Ω)) is solution of (3.7)-(3.10).

( 4 .

 4 6) {entropy} {entropy}

1 . 4 . 2 .

 142 8)-(3.10) instead of the original problem (3.7)-(3.10). It is obvious that the result of existence of the regularized problem results from the previous Lemma 4.Proposition Let δ > 0 and h ∈ W (0, T, Ω x ) a solution of (3.8)-(3.10), it satisfies u L ∞ (0,T ;H 1 (Ω)) + u H 1 (0,T ;L 2 (Ω)) ≤ M u ,(4.8){estimation_u_p} {estimation_u_p} where u = ((hh bot ) + δ) 2 and constant M u only depends on the data of the problem.

. 12 ) 4 . 4 .

 1244 RemarkThe characterization (4.11) of r depends on the function g which is the norm of the inverse of the Laplacian, g (p) = ||Λ|| L (L p (W -1,p );L p (W 1,p 0 )) and which could appear hard to compute explicitely. Nevertheless (see the proof of Lemma 3.2), mentioning g (r ) is only necessary when g (r ) > g (2) = 1. Thus (4.11) requires g (r )(1 -

(4. 13 )

 13 Our aim is to check that F (D) ⊂ D for some appropriate choice of M . Let h ∈ D and h = F ( h). Applying Lemma 3.2 to (4.1), we deduce that, with the notations of (3.3) and (4.11),||∇h|| (L r (Ω T )) N ≤ g (r ) φ ||∇ Q|| (L r (Ω x,T )) 2 + ||h 0 || W 1,r (Ω x ) (1 -k(r ))β ≤ g (r ) φ ||∇ Q|| (L r (Ω x,T )) 2 + ||h 0 || W 1,r (Ω x )(1 -g (r )(1constant M such that the initial condition and source term satisfyg (r ) φ ||∇ Q|| (L r (Ω x,T )) 2 + ||h 0 || W 1,r (Ω x ) (1 -g (r )(1 -δ δ + ))β ≤ M .Then, we obtain ||∇h|| (L r (Ω x,T )) 2 ≤ M .(4.14)

  .15) {50_1} {50_1} besides h (L r (Ω x,T )) 2 ≤ lim inf n→∞ ∇h ϕ(n) (L r (Ω 2,T )) 2 ≤ M .But we know (see the proof of the closeness of C in L 2 (0, T ; H )) thath ϕ(n) h weakly in L 2 (0, T,V ) thus in particular Ω x,T ∇h ϕ(n) • Φ d xd t → Ω x,T ∇h • Φ d xd t , ∀Φ ∈ (L 2 (Ω x,T )) 2 .

Assumptions ( 3 .

 3 11)-(3.12) are sufficient to define the primitive function B such thatB (P ) = P κ(θ(P )) (4.16) {kirchoff} {kirchoff}The application B is bijective by (3.13) and so the existence of p such that p = B (P ) is equivalent to the existence of P solution of (3.5). Applying Kirchoff's transform to (3.5), we now consider the transformed problem in the upper capillary fringeτ(p)∂ t p -K ∆p -∇ • ρg K κ(θ(B -1 (p)) e 3 = 0 in (0, T ) × Ω + h , (4.17) {equation1bis} {equation1bis} p| Γ h = B (P s ) in (0, T ), ∇p • ν = 0 on (0, T ) × (Γ soi l ∪ Γ ver ), p(0, x, z) = B (τ(p) = θ + α P θ B -1 (p) B -1 (p). Note that there exists a nonnegative real τ -such that 0 < τ -:= α P θ - κ + ≤ τ(p) ≤ τ + := α P θ + κ - . (4.19) {theta-mino} {theta-mino}

20 Remark 4 . 6 .

 2046 K p = {v ∈ X (0, T, Ω); v X (0,T,Ω) ≤ M p }, the constant M p being defined thereafter. Let p ∈ K p and h ∈ W (0, T, Ω) a solution (3.8)-(3.10), we solve the following problem:p(t , x, z) = B ρ g P s ρg + hz in (0, T ) × Ω - h , (4.20) {equation2bis} {equation2bis} τ( p)∂ t p -K ∆pρg K κoθoB -1 ( p)∂ z p = 0 in (0, T ) × Ω + h , (4.21) {equation1ter} {equation1ter} p| Γ h = B (P s ) in (0, T ), ∇p • ν = 0 on (0, T ) × (Γ soi l ∪ Γ ver ), p(0, x, z) = B (P 0 ) in Ω + h(t =0) , (4.22) 126 CHAPITRE 4. ASPECTS THÉORIQUES GLOBAL EXISTENCE RESULT FOR A MODEL OF SHALLOW AQUIFER IN INTERACTION WITH OVERLAND WATER. Taking into account the definition of the pressure P in the saturated part (2.22), we can give another expression for the continuity of the flux across Γ h . More precisely we have

Lemma 4 . 7 . 1 .

 471 Let δ > 0 such that g (4) + δ < Let p ∈ K p and h ∈ W (0, T, Ω x ) be a solution of (3.8)-(3.10) satisfying (4.12) for r = 4. There exists a unique function p ∈ W (0, T, Ω) solution of (4.21)-(4.22) that satisfies p L ∞ (0,T ;H 1 (Ω)) + p H 1 (0,T ;L 2 (Ω)) ≤ M p , (4.24) {estimation_p_u} {estimation_p_u} where constant M p only depends on the data of the problem. Proof. Let δ > 0, h ∈ W (0, T, Ω x ) a solution of (3.8)-(3.10

Lemma 4 . 8 .

 48 )) ≤ C (ρ, θ, κ, τ -, δ, M u , M , |Ω x |, T, p 0 ). (4.25) {estimation_p} {estimation_p} Finally we infer from (4.20) and (4.10) the estimate p L ∞ (0,T ;H 1 (Ω - h )) + p H 1 (0,T ;L 2 (Ω - h )) ≤ C . The system (4.21)-(4.22) is a linear parabolic problem with smooth coefficient on variable domain. By adapting the proofs of [26]-[30], the estimate (4.25) is sufficient to assert the existence of a solution p of problem (4.21)-(4.22). The uniqueness of p results from the linearity of the problem. The Lemma 4.7 is proved. Let p ∈ K p and h ∈ W (0, T, Ω) a solution of (3.8)-(3.10) satisfying (4.10), Lemma 4.7 enables to define an application F : W (0, T, Ω) → W (0, T, Ω) by F ( p) = p, (4.26) {def_F} {def_F} where the function p is the solution of (4.21)-(4.22). The end of the present subsection is devoted to the proof of the existence of a fixed point of F in some appropriate subset. Let h ∈ W (0, T, Ω x ) a solution of (3.8)-(3.10), thus• The application F defined by (4.26) is sequential continuity in L 2 (0, T ; H ) when restricted to any bounded subset of W (0, T, Ω). • There exists C a nonempty, closed, convex, bounded set in L 2 (0, T ; H ), defined such that F (C ) ⊂ C .

  .3) {eq3} {eq3} System (3.5)-(3.7) is completed by the following boundary and initial conditions:∇h ε • ν = 0 on (0, T ) × ∂Ω x , h ε (0, x) = h 0 (x) in Ω T ,(5.4)P ε | Γ h = P s in (0, T ), ∇P ε • ν = 0 on (0, T ) × (Γ soi l ∪ Γ ver ), P ε (0, x, z) = P 0 (x, z) in Ω + h 0 ,(5.5)whereθ ε = θ + ε and κ ε = κ + ε, (5.6) {theta_reg} {theta_reg} functions θ and κ satisfying (3.11) and (3.12).

  Nous utilisons maintenant les approximations introduites dans 1.2.4 pour intégrer verticalement l'équation (1.2.11), réduisant ainsi le problème 3d à un problème 2d . Nous effectuons l'intégration verticale entre les profondeurs h bot et h. Puisque θ(P ) = φ dans la zone saturée, la moyenne verticale de (1.2.11) conduit à

	h
	h bot

de continuité à travers l'interface Notre

  .2.18) A ce stade, nous avons obtenu un système indéterminé de deux équations aux dérivées partielles ((1.2.10)-(1.2.18)) avec trois inconnues P , H et h. objectif est maintenant d'inclure dans le modèle les propriétés de continuité et de transfert à l'interface. En conséquence, nous allons exprimer les deux termes de flux apparaissant dans l'équation (1.2.18) pour réduire le nombre d'inconnues.

	Équations

  17) Il s'agit de la forme la plus simple du modèle (1.3.7)-(1.3.11) puisque (1.3.16) est un problème aux limite classique du fait que h soit fixé. On remarque que la condition au bord sur Γ soil est identique à celle donnée dans le problème de Richards 3d original. Par ailleurs, la condition sur le bord Γ bot est à présent de type Dirichlet. Celle-ci permettra le bon couplage entre les composantes lente et rapide de l'écoulement. En fait, même si cette condition est une condition de Dirichlet, il n'est pas permit à l'eau de quitter l'aquifère par le fond. En effet, le flux (u • e 3 )| Γ bot (pas forcément nul) apparaît comme terme source dans la première équation de (1.3.17). Ce terme source est calculé par l'intermédiaire du problème de Richards 1d vertical et transfère l'eau depuis sa description verticale vers sa description horizontale. Ces dernières correspondent respectivement aux composantes rapide et lente de l'écoulement. Regardons plus précisément ces vitesses. H joue le rôle d'une charge hydraulique et ne dépend pas de z. Dans cette situation, la vitesse associée à P a est caractérisée par une loi de Darcy associée à la conductivité k r (P a ) M 0 . La pression auxiliaire P a satisfait

					∂ ∂z	k r (P a ) K zz	∂P a ∂z	+ 1 = 0	dans ]0, T [×Ω.
	Il s'agit de la version stationnaire de l'équation (1.3.18). L'évolution de P a est donnée par la première équation de (1.3.17) par le biais de H . Ainsi, le couple (P a , w) repré-
	sente grossièrement la composante lente de l'écoulement qui nécessite une longue
	durée d'expérience pour apparaître non négligeable. À l'inverse avec ce type de du-
	rée, l'écoulement vertical apparaît instantané.
	Composante rapide de l'écoulement : L'écoulement rapide est globalement verti-
	cal. La vitesse de l'eau et caractérisée par u dans (1.3.15). On peut voir cette vitesse
	comme provenant d'une loi de Darcy 3d (comme dans (1.2.25)) dont la partie hori-
	zontale serait négligée.
	φ	∂s(P ) ∂t	-	∂ ∂z	k r (P ) K zz	∂P ∂z	+ 1 = 0	dans ]0, T [×Ω.	(1.3.18)

La vitesse v de l'écoulement s'avère être la superposition des deux vitesses auxiliaires u et w. D'après

(1.3.16)

, on a que la pression P satisfait l'équation de Richards 1d verticale suivante dans laquelle la variable x ∈ Ω x n'apparaît que comme un paramètre.

Finalement, on peut remarquer que le problème de Richards 3d se réduit à l'équation précédente lorsque les termes de diffusion horizontale sont négligés. En effet, en temps court (pour lequel seulement la composante rapide de l'écoulement est visible) cette diffusion horizontale s'avère être non-dominante.

Composante lente de l'écoulement : La composante lente de l'écoulement est globalement horizontale et admet pour vitesse w. Cette inconnue dépend de la pression auxiliaire P a (t , x, z) := H (t , x) -z par l'intermédiaire

de (voir (1.3.15)

)

w = -k r (P a ) M 0 ∇ x H .

où 1.3.3. Commentaires sur le modèle dans les cas (1.3.13) et (1.3.14) Dans cette partie, on revient au modèle général (1.3.7)-(1.3.11) dans lequel l'interface virtuelle h(t , x) dépend de t et vérifie

Composante rapide de l'écoulement : globalement vertical, une partie est instan- tanée

  . L'écoulement est décrit comme précédemment dans Ω + h . Il y est donc globalement vertical et vérifie la même équation de Richards 1d. La vitesse associée est encore donnée par une loi de Darcy 1d. Contrairement au cas précédent, la partie Ω - h est à présent non-vide et décrit une partie saturée de l'aquifère. De plus, la pression y suit le profil affine donné par (1.3.9) dans tout Ω - h . L'écoulement vertical apparaît comme étant instantané. Il est associé à une vitesse u nulle. Cette hypothèse de vitesse verticale nulle dans la partie saturée est classique dans le cadre des écoulements dans des aquifères peu profond. Celle-ci est l'Hypothèse de Dupuit (voir partie 1.2.4). Le problème (1.3.7)-(1.3.11) dans les cas (1.3.14) peut donc être vu comme le couplage d'un écoulement horizontal de Dupuit avec de nombreux problèmes de Richards 1d. Ces derniers permettant de décrire plus précisément la recharge de la nappe phréatique par de l'eau provenant de la surface.

	h peut être vu h n'est pas vide par comme la nappe phréatique. D'autre part, la partie supérieure Ω + construction puisque h(t , x) ≤ h max . L'intérêt est donc qu'il n'est jamais nécessaire de décrire un couplage entre Ω -h et l'écoulement de surface, celui couplant Ω -h et Ω + h sera suffisant.

3.21)

On en déduit en particulier, puisque l'on a R ≥ 0 et d'après (1.2.15), que le sol est complètement saturé dans Ω - h (t ) pour tout t ∈]0, T [. Ainsi lorsque R = 0, Ω -

.4. ASPECTS NUMÉRIQUES DU PROBLÈME COUPLÉ

	1L'objectif de cette partie est double.
	-Le premier est de proposer un schéma numérique adapté à l'approximation
	du problème couplé (1.3.7)-(1.3.11). Celui-ci sera basé sur une reformulation
	du problème couplé dans laquelle ne subsistera qu'une unique équation 2d.
	Le couplage initial ce réduira à un terme faisant apparaître un opérateur "Di-
	richlet to Neumann".
	-Le second objectif se base sur l'utilisation de ce schéma afin de comparer nu-
	mériquement les problèmes couplés, en fonction du choix de Q dans (1.3.5),
	avec le problème original de Richards 3d. On montrera en pratique que dans
	toutes les situations testées, il n'est pas nécessaire que l'aquifère admette un
	ratio ε=profondeur/largeur très petit pour que l'approximation donnée par le
	problème couplé soit satisfaisante. De plus on montre qu'un bon choix de la
	valeur R dans (1.3.5) permet d'améliorer les approximations.
	1[×Ω x	(1.3.25)
	Ces problèmes effectifs déterminent les comportements dominants de l'écoulement
	dans un aquifère peu profond, et ce, en fonction de l'échelle de temps considérée.
	Ces comportements sont la base du modèle couplé (1.3.7)-(1.3.11).
	Le résultat principal de cette partie est que les problèmes effectifs issus du pro-
	blème couplé (1.3.7)-(1.3.11) sont exactement les mêmes que ceux ci-dessus. Ce ré-
	sultat est montré dans le chapitre 2. Il justifie que le modèle couplé approche bien
	(asymptotiquement) le problème original de Richards 3d dans des aquifères peu pro-
	fond quelque soit l'échelle de temps considérée.
	Pour quantifier cette approximation, nous allons comparer numériquement ces
	deux modèles. C'est l'objet du chapitre 3 qui est résumé en français dans la section
	suivante.

1.4.1. Description du schéma implicite en temps Dans

  l'approximation numérique du problème couplé, la tache non-triviale est la linéarisation du couplage non-linéaire en temps. Ainsi, on se concentrera uniquement dans la suite sur le schéma en temps et l'on ne spécifiera pas celui en espace.On introduit la discrétisation suivante. Pour M ∈ N * , les temps discrets sont donnés par t n = n ∆ t pour n = {0, ...., M } et ∆ t = T /M . Les inconnues discrètes au temps t n pour n ∈ {0, ...., M } fixé sont

	P n (x, z) P (t n , x, z),	H n (x) H (t n , x),
	h n (x) h(t n , x)	et	u n (x, z) u(t n , x, z).
	On utilise un schéma d'Euler implicite pour approcher les dérivées en temps dans
	(1.3.7)-(1.3.10) and (1.3.14). On obtient le problème discret suivant qui caractérise (P n , u n 3 , H n , h n ) pour n ∈ {1, ...., M } et P n-1 donnés.

Présentation du premier cas test. Dans ce test on considère le niveau du sol comme étant horizontal donné par le graphe de h soil (t , x) := 0. Au dessus de la nappe phréa- tique (dans l'aquifère) on considère un réservoir d'eau, de forme rectangulaire, dans lequel le sol est complètement saturé. Cet ensemble est noté ω est donné par

	Impermeable rock	h soil , h bot
	ω =]L x /10, 3L x /10[× ]-3.5, -1.7[ .	(1.4.11)
	Dans ce premier test, on s'intéresse seulement à l'infiltration de l'eau provenant de
	ω. Pour cette raison, nous considérons une condition de Neumann homogène au
	niveau du sol (α = F = 0 dans (1.3.8)) comme sur les bords verticaux. On représente la situation initiale dans le tracé de gauche de la Figure 1.2.
	Présentation de la seconde expérience. Dans ce cas test, la recharge en eau de la
	nappe phréatique vient de la surface par l'intermédiaire d'une condition au bord de
	Dirichlet sur Γ soil . Cette situation représente physiquement l'infiltration de l'eau sous

  14) with respect to H holds only through the boundary value at level h cons . In particular, the level h cons itself do not depends on H .A natural choice for the stabilization term Λ in (2.10) is the derivative of Θ with respect to H (if Θ enough regular). The procedure (2.10) becoming indeed the classical Newton method. With this in mind, we recalling the definition of Θ in (2.6) and make formally the derivative of Θ with respect to H . It follows since h cons does not depend on H that

	that ( ∂P ∂H , ∂u ∂H ) is solution of the following linear problem of unknown (P , u )
	∂Θ ∂H	(x, H , P ) =	∂u 3 ∂H h bot	+	1 ∆t	h soil (x) h bot (x)	φs (P )	∂P ∂H	d z,	(2.15)
	where s is the derivative of s with respect to P . It remains to describe this unknown ( ∂P ∂H , ∂u ∂H ). By computing formally the derivative of each equation of (2.14) with respect to H , we obtain

C (x, H , P ) :

  .26) {mass_DR} {mass_DR} since s(P ) = 1 in ]h bot (x), h(t , x)]. Next we set,

					M tot (t ) := M + h (t ) + M -h (t ).			(2.27) {mass_conservation} {mass_conservation}
	Proposition 2.1. Assuming an homogeneous density ρ, the total mass satisfies for all t ∈ (0, T ) : ∂ ∂t Ω M tot = ρ Q.
	Proof. By using relation (2.26) and (2.27) it comes
	∂ ∂t	M tot =	Ω x	h soil (x) h(t ,x)	ρ ∂ t θ(P ) + θα P ∂ t P d z d x +	Ω x	h(t ,x) h bot (x)	ρ S 0 B f ∂ t h d z d x	(2.28) {proof_mass_1} {proof_mass_1}

  Lemma 4.1. Let h 0 ∈ H 1 (Ω x ) and Q ∈ L 2 (0, T ; Ω x ), there exists a function h ∈ W (0, T, Ω x ) solu-

	tion of (3.7)-(3.10) that satisfies
	h L 2 (0,T ;H 1 (Ω x )) ≤ M and h L 2 (0,T ;(H 1 (Ω x )) ) ≤ M ,
	where M and M only depend on the data of the problem.
	Futhermore the following maximum principle holds true:
	h bot ≤ h(t , x) for a.e. x ∈ Ω x and for any t ∈ (0, T ).

  τ, x) d x -1 2 Ω x h 2 (0, x) d x.

	So we obtain
	φ 2 Ω x	h 2 (τ, x) d x +	Ω x,τ	T δ ( h) K ∇h • ∇h d xd t =	φ 2 Ω x	h 2 (0, x) d x +	Ω x,τ	QT

l ( h)h d xd t . Besides Ω x,τ K T δ ( h)∇h • ∇h d x ≥ K δ||∇h|| 2 L 2 (0,τ;H ) ,

and

| Ω x,τ QT l ( h)h d xd t | ≤ Ω x,τ |h| 2 d xd t + l 2 Ω x,τ Q2 d xd t .

  39, Corollary 9.7]. There exists h ∈ C such that F (h) = h. Then h is a weak solution of problem (3.8)-(3.10).

	• Step 2 : Maximum Principle We claim that
	h bot ≤ h(t , x) a.e in Ω x,T .	(4.5) {ppemax} {ppemax}

  Moreover, ∇h • ∇h m = χ {h bot -h>0} |∇(h -h bot )| 2 . Thus, taking w(t , x) = h m (x, t )χ (0,τ) (t ) in (3.8) GLOBAL EXISTENCE RESULT FOR A MODEL OF SHALLOW AQUIFER IN INTERACTION WITH OVERLAND WATER. 16 τ 0 φ〈∂ t h, h m (x, t )〉 V 1 ,V 1 d t + QT l (h) h m d xd t .By definition of T l (h), T l (h)χ {h<h bot } = 0, we can simplify the above equation as follows {h<h bot } δ K |∇(h -h bot )| 2 d xd t . = 0

		Ω x,τ	(δ + T l (h)) K ∇h • ∇h m d xd t =	Ω x,τ
	φ 2 Ω x	h 2 m (τ, x)d x +	Ω x,τ
	The previous equation leads to		1 2 Ω x,τ	h 2 m (τ, x)d x ≤ 0
	and then h m = 0 a.e. in Ω x,T .
	• Step 3 : Uniform estimate

χ

  h bot ) + δ) d x +C Q L 2 (Ω x,T ) , |∇h| 2 d xd t ≤ |Ω x | S ∞,(0,h soi l -h bot +1) +C Q L 2 (Ω x,T ) .(4.7)

	Hence
	Ω x S ((h 0 thus (since δ < 1) Ω x,T |∇h| 2 d xd t ≤ Ω x
	K
	Ω x,T

1 d t + Ω x,T K ∇h • ∇h d xd t = Ω x,T QT δ (h) S ((hh bot ) + δ) d xd t .

We infer from Lemma 3.1 that

T 0 〈∂ t h, S ((hh bot ) + δ)〉 V 1 ,V 1 d t = Ω x S ((h(T, .) -h bot ) + δ) d x -Ω x S ((h 0 -h bot ) + δ) d x,

and (since the function T

δ (h) ln T δ (h) is bounded over (h bot , h soi l )) | Ω x,T QT δ (h) S ((h-h bot )+δ) d xd t | ≤ Q L 2 (Ω x,T ) T δ (h) S ((h-h bot )+δ) L 2 (Ω x,T ) < C Q L 2 (Ω x,T ) .

S ((h(T, .) -h bot ) + δ) d x + K

  by ∂ t u and integrate by parts over Ω x . We thus obtain

	φ	Ω x	((h -h bot ) + δ)∂ t h∂ t u d x +	K 2	d d t Ω x	|∇u| 2 d x =	Ω x	T l (h) Q∂ t u d x	(4.9)
	Besides
					|	Ω x	T ≤	l 2 φ	Q 2 L 2 (Ω x ) +	φ 4	∂ t u 2 L 2 (Ω) ,
	hence
	φ 4	0	T	Ω x	|∂ t u| 2 d x +	K 2 Ω x	|∇u| 2 d x ≤	l 2 φ	Q 2 L 2 (Ω x ) +	K 2 Ω x	|∇u 0 | 2 d x.
	We get inequality (4.8) by taking M 2 u := sup( We directly infer from (4.8) that	2 φ	,	1 K )	2 l 2 φ	Q 2 L 2 (Ω x,T ) + K	Ω x	|∇u 0 | 2 d x .
							h L ∞ (0,T ;H 1 (Ω)) + h H 1 (0,T ;L 2 (Ω)) ≤	M u δ	.	(4.10) {estimation_h_p}

l (h) Q∂ t u d x| ≤ Ω x |T l (h) Q| 2 d x × ∂ t u 2 L 2 (

Ω) , {estimation_h_p} Proposition 4.3. Let h ∈ W (0, T, Ω x ) be a solution of Problem (3.8)-(3.10), and let r > 2 be the greatest real number such that k(r ) = g (r )(1µ) < 1. (4.11) {choix_r} {choix_r}

  ). Assume p is a solution of (4.21)-(4.22). In the following, we denote by C a generic constant. We multiply Eq. (4.21) by ∂ t p and integrate by parts over Ω + h . We getSince P s is constant Γ bot ∂ t B (P s )∇p. νd σ = 0.Using Cauchy-Schwarz and Young inequalities and the regularity of functions θ and κ, we get ,x, h)| 2 ∂ t hd x| ≤ ρ 2 g 2 κ + M 2 + |Ω x | 1/2 M u δ .Choosing = τ -and applying Gronwall lemma, we deduce

		Ω + h	τ( p)|∂ t p| 2 d z d x +	1 2	d d t Ω +
	|	Ω +							Ω + h	|∂ z p| 2 d z d x	1/2	Ω + h	|∂ t p| 2 d z d x	1/2
					≤ C (ρ, θ, κ)	Ω + h	|∇p| 2 d z d x	1/2	Ω + h	|∂ t p| 2 d z d x	1/2
					≤ 2 Ω +
							Ω + h	|∇h| 4 d z d x) 1/2 + |Ω x | 1/2	Ω x	|∂ t h| 2 d x	1/2 .
		So, taking into account inequalities (4.10) and (4.14) written for r = 4, we get
		| |∇p(t τ -Ω x 2 T 0 Ω + h |∂ t p| 2 d x d t + 1 2 Ω + h	|∇p| 2 d x ≤	C 2 Ω + h	|∇p 0 | 2 d x e	T C (ρ, θ, κ, δ, M u , M , |Ω x |) τ -	,
	therefore
		p L ∞ (0,T ;H 1 (Ω + h )) + p H 1 (0,T ;L 2 (Ω + h

h |∇p| 2 d z d x -Γ bot ∂ t B (P s )∇p. νd σ -Ω x |∇p(t , x, h)| 2 ∂ t hd x = Ω + h ρg κoθoB -1 ( p)∂ z p∂ t pd z d x. h ρg κoθoB -1 ( p)∂ z p ∂ t pd z d x| ≤ C (ρ, θ, κ) h |∂ t p| 2 d z d x + C (ρ, θ, κ) 2 Ω |∇p| 2 d z d x.

Besides, using (4.23), we can write

| Ω x |∇p(t , x, h)| 2 ∂ t hd x| = Ω x ρ 2 g 2 B (P )|∇h| 2 + 1 ∂ t hd x ≤ ρ 2 g 2 κ + (

  Ω T ) and a.e. in Ω T .

	We thus introduce the truncature function T P d such that
	Let	T P d (x) =	x P d if x < P d . if x ≥ P d ,	(5.10)
		t ε = T P d (P ε ).	(5.11) {def_t} {def_t}

however the numerical simulations below show good results even for a ratio of order 0.1, which is not exceeded by the large majority of the unconfined aquifers.7

This profile has to be compared with the corresponding one in Figure2which is obtained with the same physical parameters. The idea is that at almost the same amount of water has infiltrates in the soil (5% of difference) but some of this water has flee in the horizontal direction. In conclusion, the hypothesis of only vertical flow in Ω + h which is done in the coupled model seems to be less valid in this situation. A larger ratio L x |h sat | is needed to minimized the

Remerciements

We end the proof by noting that, thanks to the equality α = 0 and the first equation of (4.55) for k = 2, we have

The paper is organized as follows: In Section 2, we present the model coupling 3D-Richards and Dupuit horizontal flow. The main results regarding existence theorems are given in Section 3. Finally, the proofs of the theorems are performed in Sections 4 and 5.

NUMERICAL ASPECTS OF A CLASS OF MODELS COUPLING SLOW AND FAST FLOW IN SHALLOW AQUIFERS

21 choices for the function h satisfying (1.18): the minimal one for R = +∞ (and so h = h bot ), the maximal one for R = 0 and the intermediate one for R = 7.

We denote by h 2d sat the level obtained from the reference 2d-Richards model and we denote by h R sat the one coming from the model (1.14)- (1.18) with the function h characterized by R ≥ 0. In this test, we fix the largeness of the aquifer as L x = 40. Moreover we choose the initial level of the water table to be h 0 sat = -2.2. The functions h 2d sat and h R sat for R ∈ {0, 7, ∞} are plotted in Figure 4 at time t ∈ {0.15, 1.5, 5, 10} (in days). We of course do not plot the initial situation which is the same for each model and is the one of the reference test case described in the previous paragraph. The curve h 2d sat is the reference one and is plotted with a black solid line in Figure 4.

Bear in mind that the function h characterizes the level below which the vertical flow is assumed to be instantaneous (instead of being described by the 1D-Richards equation). In every case, the horizontal flow is ruled by equation (1.17).

• In the case R = 0 we have h = h 0 sat . In this case the vertical flow in the whole saturated zone Ω - h sat is considered to be instantaneous. Then, when the water coming from the supply ω reaches the water table, the flux (u •e 3 )| Γ h increases very quickly. So does the corresponding hydraulic head H and the horizontal flow is very and even too fast (see the black dotted line compared to the black solid line in Figure 4). • In the case R = +∞, we have h = h bot (see (1.13)). The vertical flow is described by the 1D-Richards model in the whole domain, even in the saturated part below the level z = h ∞ sat . The horizontal flow in this case seems to be slower than the one given by the Richards model (compare the gray dot-dashed line with the black solid one in Figure 4). Roughly the idea is that in this case the water have to travel along the whole vertical direction before reaching the level z = h = h bot . Then the flux (u •e 3 )| Γ bot at the bottom of the aquifer takes a lot of time to increase when the water coming from rectangle ω reaches the water table. This flux being the source term in equation (1.17), the function H increases with some delay and the corresponding horizontal flow is slower. • In the case R = 7, h bot ≤ h ≤ h 7 sat . The corresponding flow should exhibit an intermediate behavior between the too previous ones. Here, the value R = 7 was chosen so that h 7 sat is very close to the reference one h 2d sat (see the gray dashed line).

Error v.s. ratio: whole space-time. We represent in Figure 5 the error 1

for R ∈ {0, 3, 7, ∞} and for L x ∈ [20,160]. The initial water table is given by h 0 sat = -3.6. As expected, this error decrese when the ration L x /|h bot | increase. Of course a good choice of the parameter R provide a reduction of the error. Moreover, here the situation is simple since the value R = 7 remains better independently on the ratio.

SECOND EXPERIMENT: CONFINED AQUIFER AND RIVER

Geometry. In this second experiment, the geometry of the aquifer is again given by (1.2). As previously the horizontal component is the interval ]0, L x [ for L x > 0. Moreover, the bottom of the aquifer is characterized by an horizontal boundary h bot (x) = h bot < 0. In the other hand we consider an almost horizontal soil with a hole in the middle. The corresponding function

versus the ratio length/deepness of the aquifer (R ∈ {0, 4, 8, ∞}). Function h 2d sat is the iso-pressure P = P s in the original 2d-Richards problem and h κ sat is the one associated with the coupled problems for different choices of R characterizing h by (1.13). The test case is the one of Figure 6. error. To complete the investigation of this case, we quatify as before the total error between evolutions of the water table and also the influance of parameter R.

Dependence with respect to R. As for the first experiment, we compare in this part the solution of the classical Richards model with the one obtained by using the coupled model (1.14)- (1.18). We test the same three particular functions h satisfying (1.18): R = +∞, R = 0 and R = 7.

We use the same notation h 2d sat for the level obtained from the reference 2d-Richards model and h R sat for the one coming from the model (1.14)-(1.18) with the function h characterized by R > 0. In this test, we fix the largeness of the aquifer as L x = 40 as in the Figure 6. The initial level of the water table is also the same as in the expertiment of Figure 6 and is h 0 sat = -3.6.

The functions h 2d

sat and h R sat for R ∈ {0, 7, ∞} are plotted in Figure 8 at time t ∈ {1.5, 3, 5, 10} (in days). The curve h 2d sat is the reference one and is again plotted with a black solid line in Figure 8.

The same kind of behavior is obtain in this second experiment. We notice nevertheless that the influance of the value R is here more significant.

• In the case R = 0 we have h = h 0 sat . In this case the vertical flow in the whole saturated zone Ω - h sat is considered to be instantaneous. The resulting horizontal flow is again too fast. This continue during the whole flow experiment until reaching at the end a solution quite far from the general one (see the black dotted line compared to the black solid line in Figure 8). • In the case R = +∞, we have h = h bot (see (1.13)). The horizontal flow is again slower than the one given by the Richards model. Here the difference is significant and increase with the time, the solution of the becoming bader (compare the gray dotdashed line with the black solid one in Figure 8).

As before this intermediate choice gives an intermediate behavior which is closer to the original one. Here, the value R = 3 was chosen to minimise the error. It turns out as we will see in Figure 9 that this choice depends a priori on the value ratio L x |h bot | .

DERIVATION OF THE MODEL

The basis of the modeling is the mass conservation law written for fresh water coupled with the classical Darcy law for porous media. Fluid and soil are considered to be weakly compressible.

For the three-dimensional description, we denote by x := (x, z), x = (x 1 , x 2 ) ∈ R 2 , z ∈ R, the usual coordinates.

CONSERVATION LAWS

We begin with the conservation of momentum. In view of the (large) dimensions of an aquifer (related to the characteristic size of the porous structure of the underground), we consider a continuous description of the porous medium.

The effective velocity q of the flow is thus related to the pressure P through the Darcy law associated with a non-linear anisotropic conductivity

where ρ and µ are respectively the density and the viscosity of the fluid, K 0 is the permeability of the soil, κ(P ) is the relative conductivity and g the gravitational acceleration constant.

Introducing the hydraulic head H defined by

we write the previous equation as follows:

.2) {darcy0} {darcy0}

In this relation, the matrix K is the hydraulic conductivity which expresses the ability of the underground to conduct the fluid. We have denoted by ρ 0 the reference density of the fluid.

Next, the conservation of mass during displacement is given by the following equation

where Q denotes a generic source term (for production and replenishment). The function θ is the volumetric moisture content defined by

where φ is the porosity of the medium and s is the saturation. If we assume that the air present in the unsaturated zone has infinite mobility, the saturation s and then the function θ are thus considered as monotone functions depending on the pressure as we will detail latter.

STATE EQUATION FOR THE FLUID COMPRESSIBILITY

We consider that the fluid are compressible by assuming that pressure P is related to the density ρ as follows:

.4) {alpha} {alpha}

The real number α P ≥ 0 is the fluid compressibility coefficient and P 0 is the pressure of reference. Further assuming α P = 0 we would recover the incompressible case.

By definition of θ, H and T P d , we remark that

and the last result of convergence can be written as follows

Function H is of course not bijective in (-∞, P d ). We thus define a bijective continuous extension H of

we have H Estimates (4.8) being uniform with respect to ε, we directly pass to the limit in Eq. ( 5.3) and Eq.

(5.2) when ε → 0. Theorem 3.7 is proved.