Yves Khaled

Thierry Guy

Chiara Simeoni

Christian Bourdarias

Stéphane Gerbi

Enrique Zuazua

Merci Pour

Mehmet Ersoy

Courant, Friedrichs Warnecke Muscl

Monotone Upwind

Serre Sgn

Naghdi Green

Xviii

Discontinuous Galerkin PFS Pressurised and Free Surface AMR Adaptive Mesh Refinement AMM Adaptive Moving Mesh BB-AMR Block-Based Adaptive Mesh Refinement CPU Central Processing Unit CFL

Garde en mémoire que l'avenir te sera plus doux si tu emploies bien le temps présent.

-Les proverbes et pensées grecques (1812)

R É S U M É

Ce mémoire est le fruit de mes dix dernières années de recherche à l'IMATH. Il concerne l'analyse mathématique, numérique des équations compressibles et incompressibles dans lequel j'aborde les équations hydrostatiques et les équations non-hydrostatiques.

Je présente et je développe des outils numériques de résolution adaptative sur des maillages non structurés en dimension quelconque pour des méthodes volumes finis et Galerkin discontinus.

Je propose également la réduction de modèle, i. e. une approximation asymptotique, dans le cadre des ondes longues pour dériver ou revisiter des modèles non-linéaires hydrostatiques et non-hydrostatiques. En particulier, je propose la première généralisation des équations de Serre-Green-Naghdi monodimensionnel en géométrie quelconque pour l'hydraulique fluviale. Celle-ci est obtenue par une réduction 3D vers 1D.

A B S T R A C T

This dissertation is the fruit of my last ten years of research at IMATH. It concerns the mathematical, numerical analysis of compressible and incompressible equations in which I tackle hydrostatic and nonhydrostatic equations. I present and develop numerical tools for adaptive mesh refinement on unstructured meshes in any dimension for finite volume and discontinuous Galerkin methods.

I also propose model reduction, i.e. asymptotic approximation, in the long wave framework to derive or revisit hydrostatic and nonhydrostatic nonlinear models. In particular, I propose the first generalisation of the one-dimensional Serre-Green-Naghdi equations for arbitrary geometry for river hydraulics. This is achieved by a 3D to 1D reduction. contents iv cfd: dg method, hp-amr and applications to wavedriven groundwater flows [START_REF] Benjamin | The stability of solitary waves[END_REF] Exemples de résultats de simulation (1). . . . . . Figure 3. [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] Exemples de résultats de simulation (2). . . . . . Figure 3. [START_REF] Baron | Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation[END_REF] Cas test 1d d'infiltration. . . . . . . . . . . . . . . Figure 5. [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF] Numerical illustration of the sensitivity of (v ε , δϕ ε ). 46 Figure 6. [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF] Multiple steady shocks (dashed line: source term, solid line: numerical solution). . . . . . . . . . . . Figure 6. [START_REF] Altmann | A local time-stepping discontinuous Galerkin algorithm for the MHD system[END_REF] A steady solution emerging in finite time with g c (dashed line: source term, solid line: numerical solution). . . . . . . . . . . . . . . . . . . . . . Figure 6. [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF] A local finite time emerging steady solution with g dc (dashed line: source term, solid line: numerical solution). . . . . . . . . . . . . . . . . . Figure 7. [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF] Example of hierarchical dyadic tree. . . . . . . . Figure 7. [START_REF] Altmann | A local time-stepping discontinuous Galerkin algorithm for the MHD system[END_REF] Mesh refinement and mesh coarsening: construction of the fluxes. . . . . . . . . . . . . . . . Figure 7. [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF] Illustration of a global time advancement with a first order scheme at time t n + δt n . . . . . . . . Figure 7. [START_REF] Arpaia | r-adaptation for Shallow Water flows: conservation, well balancedness, efficiency[END_REF] Illustration of a global time advancement with a first order scheme at time t n + 2δt n . . . . . . . . Figure 7. [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] Illustration of a global time advancement with a first order scheme at time t n + 3δt n . . . . . . . . Figure 7. [START_REF] Baron | Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation[END_REF] Illustration of a global time advancement with a first order scheme at time t n + 4δt n . . . . . . . . The function f for the mesh refinement criterion S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25 and x B = 3.75. . . . . . . . . . . Figure 8. [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF] The mesh refinement criterion S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25 and x B = 3.75. . . . . . . . . . . . . . . . . . . . . . . Figure 8. [START_REF] Arpaia | r-adaptation for Shallow Water flows: conservation, well balancedness, efficiency[END_REF] Numerical results for the water height at time t = 2 s using the automatic threshold α PE . N m is the mean number of cells used during the simulation. . . . . . . . . . . . . . . . . . . . . . . Figure 8. [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] Numerical results for the water height at time t = 2 s using the criterion 1 with the mean method. N m is the mean number of cells used during the simulation. . . . . . . . . . . . . . . . Figure 8 Convergence rate of the L 2 error for ε = 0.4. The reference solution is computed with 10 000 cells. The order is computed through the first order interpolation polynomial . . -Marie-Madeleine Davy

1 P R É A M B U L E
Le rêve de tout scientifique n'est-il pas de comprendre et de prévoir les réactions de la nature qui l'entoure? Depuis toujours, connus sous différents noms : savant, "cheikh" ou encore chercheur, ces hommes ont établi des lois reliant des grandeurs observables à des grandeurs non ou moins observables, en se basant simplement à leur intuition.

La seule vraie chose précieuse est l'intuition.

-Albert Einstein L'idéologie de ces savants était basée sur les quatre éléments : la terre, l'eau, l'air et le feu. Cette idéologie, bien que simpliste, avait une signification profonde, qui dans les faits n'est pas tellement éloignée de nos connaissances à ce jour. Cette tétralogie se traduit par :

• la terre désignant l'état solide,

• l'eau désignant l'état liquide,

• l'air désignant l'état gazeux,

• le feu désignant l'énergie.

Entre terre et ciel, ces éléments sont omniprésents et en particulier ces "fluides" sont en nous, et, leurs rôles est essentiel puisqu'ils sont les supports de toutes vies.

Contrairement à l'eau, l'air est omniprésent et fort heureusement. On la trouve en abondance dans certaines régions alors que dans d'autres elle est une denrée rare. C'est ainsi que la recherche de ce liquide met l'hydraulique au premier rang des sciences les plus anciennes. Il est facile de s'en convaincre. De la plus lointaine civilisation à celle de nos jours, l'homme a laissé une trace de son action sur l'eau. Les canaux, les aqueducs des mésopotamiens aux romains en sont des exemples, mais pas seulement, l'homme antique fit entrer ce "fluide" dans le panthéon de ces divinités. Outre le caractère vital de l'eau, ce "fluide" a toujours été un élément de fascination parfois religieuse et parfois scientifique.

préambule La notion de "fluide" est bien-sûr bien plus générale que celui de l'eau, bien que celui-ci reste dans tous les esprits comme étant l'exemple typique. L'air et plus généralement tous les gaz sont des exemples plus subtils. On peut étendre à d'autres corps dont les métaux en fusion, la lave des volcans, les roches, la matière stellaire constituant une nébuleuse, ou encore une foule d'individus s'écoulant dans une rue ou encore le trafic routier etc. permettant de classer les fluides en fluide newtonien, un fluide qui garde sa fluidité même si on la remue contrairement au fluide non-newtonien. Le fluide non-newtonien permet de marcher dessus et sa viscosité augmente lorsque la contrainte mécanique augmente.

Tout au long des siècles, presque exclusivement entre les mains des artisans, comme l'est le travail du fer au forgeron, la mécanique des fluides d'Archimède à Vinci, a été et est une contemplation esthétique et un sujet de recherche vaste et fascinant dans tous les domaines scientifiques. De l'invention rudimentaire à l'invention moderne, de la roue hydraulique à l'éolienne en passant par les moulins à vent, à eau, et les premières "machines volantes", de nombreuses créations ont permis d'exploiter le fluide non pas comme une ressource vitale, mais comme un élément utile à rendre notre quotidien meilleur.

Par cette intuition encore et toujours, l'artisan a créé et innové. Il a fallu attendre plusieurs siècles pour voir le mariage entre le savant et l'artisan. On peut citer à titre d'exemple, les travaux de Belidor, ingénieur militaire et professeur de mathématiques, qui en 1737 publia les quatre tomes de son "Architecture hydraulique" ou "L'art de conduire, d'élever, et de ménager les eaux, pour les différents besoins de la vie", Bernouilli, membre d'une célèbre famille de mathématiciens de Bâle, qui formule l'expression de la conservation de l'énergie où figurent les trois termes de position, de pression et de vitesse (Traité d'hydrodynamique de 1738), Toricelli, en 1644 avait précédé et limité son travail aux jets d'eau des villas italiennes, Euler en 1741 donne l'expression dite "des moments cinétiques" qui est à la base des calculs de toutes les turbomachines, ou encore Lagrange qui éclairera la seconde moitié de ce siècle.

Si Vinci fut un exemple, Euler en est un autre. On oublie généralement qu'il est parti de la théologie, pour le conduire ensuite aux mathématiques, à la topologie, à la mécanique des fluides etc. tout en restant le philosophe de son siècle. Cette polyvalence s'oppose à la spécialisation des chercheurs modernes. Mais cette universalité a une limite. Si Euler était un éminent philosophe et savant, il n'était pas ingénieur pour autant. Il faudra attendre quelques années avant que les théories de génie voient leurs réalisations se concrétiser par Fourneyron, Girard, Kaplan etc. Les réflexions de ces savants ne sont plus seulement des formules aujourd'hui, elles sont chaque jour utilisées et concrétisées, c'est la magie de la théorie que je décrirais sous la forme suivante :

La théorie est folle ! Elle n'est pas assez précise pour être juste, Mais assez pour la mettre en pratique.

-

Mehmet Ersoy

Loin d'un fluide réel, les fluides théoriques dits "fluides parfaits" ignorent toutes les difficultés que présente l'autre. Le phénomène essentiel qui caractérise le fluide réel est la combinaison des effets cinétiques d'inertie et de la viscosité génératrice de forces tangentes aux lignes de courant, forces qui vont engendrer des gradients de vitesse perpendiculairement aux parois limitant l'écoulement. Les premières études de cette la viscosité sont dues à un médecin et physicien, Poiseuille qui portera le nom de "Loi de Poiseuille" (1840). Cette loi communément appelée loi de Hagen-Poiseuille a été découverte en parallèle par Hagen. Ce dernier était un scientifique de renom connu pour ses travaux dans de multiples domaines comme le génie civile et maritime, les sciences hydrauliques et la mécanique des fluides.

Deux notions très importantes des fluides réels se déduisent de l'influence de la viscosité. C'est ici que se situent les travaux de Reynolds et en particulier son célèbre nombre sans dimension qui sépare les écoulements dits laminaires où la viscosité joue un rôle stabilisant, des écoulements dits turbulents où les forces d'inertie l'emportent et créent une instabilité interne. L'influence prépondérante de la viscosité au voisinage des parois est à l'origine de la couche limite à l'étude de laquelle s'illustrèrent entre autres Prandtl et Von Karman. Toutes les questions de sillage et décollement se déduisent de ces considérations. La mécanique des fluides "modernes" initiée par Euler et d'Alembert, le stade des fluides parfaits, fait des progrès avec Henri Navier et George Gabriel Stokes pour s'attaquer aux fluides incompressibles puis compressibles.

De la représentation concrète de l'intuition, aux pensées des grands philosophes de l'époque, tout mène indéniablement à la description par des équations, à savoir, à la modélisation mathématique, les mathématiques appliquées et le calcul scientifique dont Borel, Lagrange, Lions, etc. furent les précurseurs indéniables de la modélisation mathématique moderne.

Borel a consacré une cinquantaine d'articles à la philosophie et à l'histoire des sciences, à la psychologie, à la pédagogie et à l'économie politique. Fervent défenseur de la conception intuitionniste des mathématiques, il a toujours insisté sur la nécessité de ne jamais perdre de vue le "réel" : pour lui, Les mathématiques ne sont pas un jeu purement abstrait de l'esprit, mais sont, Au contraire, en étroite connexion avec la réalité concrète.

-Émile Borel

Il tenait à rappeler, dans les raisonnements les plus abstraits de l'analyse, comment ces raisonnements expriment des réalités, et par suite nous aident à comprendre la réalité physique. C'est pourquoi Borel a toujours préféré aux définitions descriptives les définitions constructives ; ces dernières sont indispensables pour passer de la théorie aux applications et permettent auparavant de s'assurer de l'existence et éventuellement de l'unicité de l'élément dont on parle. Jacques-Louis Lions qui était sans aucun doute le mathématicien d'exception, précurseur de biens de choses en analyse mathématique appliquée moderne et notamment en analyse numérique. Sa philosophie laisse rêveur, et comme le citait Antoine de Saint-Exupéry, "Les rêves fournissent la nourriture de l'âme. Celui qui ne cherche pas à réaliser son rêve suicide son âme." : Ce que j'aime dans les mathématiques appliquées, C'est qu'elles ont pour ambition de donner du monde des systèmes une représentation qui permette de comprendre et d'agir. Et, de toutes les représentations, la représentation mathématique, Lorsqu'elle est possible, est celle qui est la plus souple et la meilleure.

Du coup, ce qui m'intéresse, c'est de savoir jusqu'où On peut aller dans ce domaine de la modélisation des systèmes, c'est d'atteindre les limites.

-

Jacques-Louis Lions

On ne peut pas parler de mathématiques appliquées sans parler d'algorithme et sans faire allusion au père de l'algèbre. Fondateur du premier observatoire permanent au monde, pionniers de l'algèbre avec son traité, Kitab al jabr w'al muqabalah, fondateur de la maison de la sagesse, Al-Khwarizmi traduit des textes scientifiques et philosophiques du monde entier. Al-Khwarizmi était philosophe et un excellent passeur de connaissance. Il introduit dans son aire culturelle les connaissances mathématiques indiennes, notamment le système décimal de numération. La traduction latine de son ouvrage Algorithmi de numero indorum permit la transmission de ces connaissances jusque dans l'Occident du XIIe siècle. D'ailleurs, le mot Algorithmi, traduction latine d'Al-Khwarizmi, a donné naissance au mot "algorithme" et fait de lui l'irréfutable "père fondateur" de nos ordinateurs. À la fin du Xème siècle et au début du XIème siècle, le mathématicien persan Al-Karaji franchit une nouvelle étape dans l'histoire de l'algèbre : il la détache de la géométrie -ce que ni Al-Khwârismî ni Abu Kamil n'avaient pu faire -dans ses ouvrages intitulés Gloire de l'algèbre et de la muqabala et Merveilles du calcul. Ensuite, la grande étape franchie par l'algèbre musulmane est la résolution de l'équation cubique dont le célèbre poète et mathématicien persan Omar Khayyâm (1048-1131), bercé et inspiré par la poésie Ceux qui ont embrassé science et littérature ont récité leur fable et se sont endormis.

-Omar Khayyâm

Les bases de l'analyse numérique et du calcul numérique "moderne" furent, sans aucun doute, initiées par Newton (1669) sur la recherche de zéro de fonction. Dix ans plus tard, c'est à Raphson, qu'on doit l'algorithme de Newton-Raphson qui met en évidence une formule de récurrence. Un siècle plus tard, Mouraille et Lagrange étudient la convergence des approximations successives en fonction des conditions initiales par une approche géométrique. Cinquante ans plus tard, Fourier et Cauchy s'occupe de la rapidité de la convergence etc.

Du traité d'Al-Khwarizmi à la naissance de la première machine de calcul, le calcul scientifique puis le calcul haute performance a vu le jour. Frontière entre les mathématiques appliquées et la mécanique, la Mécanique des Fluides Numériques (MFN) ou par le terme anglaissaxon Computational Fluid Dynamics (CFD), a grandi d'une curiosité mathématique pour devenir aujourd'hui un outil essentiel voire incontournable dans pratiquement toutes les branches scientifiques dont la dynamique des fluides et notamment en hydrodynamique côtière et celle des vagues.

De l'université de Savoie à l'Université de Toulon, sous la direction de Christian Bourdarias et Stéphane Gerbi, puis sous celle d'Enrique Zuazua lors de mon post-doc au Basque Center for Applied Mathematics, ainsi, j'ai l'immense plaisir de vous présenter sans aucune prétention mes travaux de recherche, fruit de nombreuses années de passion.

Ce mémoire est dédié à l'étude, l'approximation et la simulation des équations aux dérivées partielles issues de la mécanique des fluides compressibles ou incompressibles -hydrostatique et non hydrostatique dans les domaines de l'hydrauliques et de l'hydrodynamique.

En amont, mes travaux comportent une part importante de modélisation -approximation de systèmes hyperboliques de systèmes de type "hydrostatique", modèles primitives compressibles, modèles de type Saint-Venant non homogène, modèles multi-fluide, modèles à faible Mach -dont les applications sont par exemples, les écoulements à surface libre et/ou en charge, les problèmes d'entrainement d'air et plus généralement les écoulements multi-fluides, la sédimentation, la dynamique de l'atmosphère, les problèmes de déferlements, les problèmes d'inondations, les écoulements non-hydrostatique . . . En aval, mes travaux abordent un certain nombre de problèmes théoriques -existence et unicité des solutions faibles, stabilité des solutions faibles, contrôle optimal, problèmes de Riemann (systèmes autonomes et non autonomes), analyse asymptotique en temps long (propriétés des solutions stationnaires vs solutions non stationnaires) -dont le but principal est d'aboutir à certaines caractérisations nécessaires lors de l'approximation numérique de ces équations aux dérivées partielles. Mes travaux explorent donc l'analyse numérique, la construction de méthodes numériques Volumes Finis en multi-dimension et le calcul scientifique -solveur de Godunov, solveur cinétique, schémas équilibres, schémas anti-diffusif, techniques de raffinement de maillage, critère de raffinement -pour la simulation numérique des problèmes de modélisation mentionnés précédemment. Cette thèse est consacrée à la mise en oeuvre de nouveaux modèles et méthodes numériques pour des écoulements compressibles ou incompressibles en domaine "couche mince". S'il est possible de calculer la dynamique d'un écoulement (par exemple en conduite fermée) à partir de la résolution des équations de Navier-Stokes, il est nécessaire d'effectuer la réduction de modèle afin de s'intéresser à la modélisation hydrodynamique aux échelles qui nous intéressent. Les modèles obtenus sont souvent plus exploitables d'un point de vue numérique en termes de gain de temps et introduit des nouveaux problèmes. En l'état actuel, les modèles concernés ne sont pas entièrement satisfaisants tant d'un point de vue théorique que numérique notamment en domaine "couches minces".

Ce travail apporte un nouveau regard à des problèmes existants par des dérivations de modèles originales, des preuves mathématiques astucieuses ou des schémas numériques alliant précision et simplicité.

motivations.

Le couplage entre un écoulement en charge et un écoulement à surface libre (écoulements mixtes transitoires en charge/à surface libre) est, comme le montre l'expérience, un problème délicat compte tenu de l'instabilité du phénomène. Celui-ci est en outre peu accessible aux mesures directes. Le suivi en temps réel, à l'aide d'une simulation numérique, est d'un intérêt majeur pour les ingénieurs en tant qu'outil de validation et de prévision. Il existe en effet des risques très importants liés à des surpressions, compromettant l'intégrité de la structure. Ce type de problème se rencontre, par exemple, dans la modélisation des écoulements dans des réseaux d'assainissement lors d'orages violents, dans des canaux d'amenée en amont des usines hydroélectriques, etc. l'état de l'art.

Nous avons proposé un nouveau modèle pour les écoulements mixtes à géométrie variable, nommé PFS (voir [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF] pour une dérivation complète ou [54, Chapitre 1]). Ce système est donné par les équations suivantes :

∂ t A + ∂ x Q = 0, ∂ t Q + ∂ x Q 2 A + p(x, A, E) = -gA d Z dx +Pr(x, A, E) -G(x, A, E) -gK(x, A, E) Q|Q| A
où t est la variable de temps, x l'abscisse curviligne (suivant l'axe d'écoulement),

d Z(x) dx
le terme de pente, Pr(x, A, E) le terme source de pression, G(x, A, E) le terme de courbure et K(x, A, E) le terme de frottement donné par la formule de Strickler (voir [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF] ou [54, Chapitre 1] pour le détail des termes). L'inconnue A(t, x), homogène à une aire, est la variable "mixte" car elle permet de calculer tantôt le niveau d'eau à surface libre, tantôt la pression en charge. L'inconnue Q(t, x) correspond au débit. La donnée de l'aire mixte A n'est en effet suffisante ni pour déterminer la loi de pression ni les termes sources. Il faut suivre l'historique du point courant en termes d'état E (surface libre si E = 0 ou charge si E = 1). Le terme de pression p(x, A, E), conforme à la physique du problème, est égal à la pression à surface libre lorsque la variable d'état E est égale à 0 (i. e. p(x, A, 0) = gI 1 (x, A) cos θ la pression hydrostatique, c. f. Figure 2.1), ou bien de type mixte hydrostatique-acoustique en cas de charge si E = 1 (i. e. p(x, A, 1) = gI 1 (x, S) cos θ + c 2 (A -S), c. f. originalité : le modèle.

Ce nouveau modèle permet de prendre en compte la compressibilité de l'eau lors d'un écoulement en charge et les points de transition (points de passage d'un écoulement en charge vers un écoulement à surface libre et vice et versa) en géométrie quelconque. Contrairement aux méthodes basées sur l'artifice de la fente de Preissmann, il permet de simuler correctement les écoulements sub-atmosphériques (considérés jusqu'alors comme un passage d'un écoulement en charge vers un écoulement à surface libre).

Ce modèle est construit à partir des équations d'Euler compressibles tridimensionnelles pour la partie en charge avec une loi de pression linéarisée et les équations d'Euler incompressibles tridimensionnelles pour la partie à surface libre. Les effets de courbure locale sont pris en compte en écrivant les équations dans un repère local mobile (suivant l'axe d'écoulement moyen). Les deux modèles obtenus sont de type Saint-Venant et sont formellement proche dans l'écriture. Ils sont couplés pour obtenir un système d'équations hyperboliques : le modèle PFS.

originalité : aspects numériques.

L'approximation numérique naturelle utilisée est la méthode des Volumes finis de par la nature hyperbolique du système PFS. La présence des termes sources "complexes" rendent la construction d'un schéma Volumes Finis consistant, stable et entropique difficile. C'est d'ailleurs la difficulté majeure de la discrétisation de ces équations. En combinant les techniques de décentrement classique pour les termes conservatifs, en utilisant la théorie des produits non conservatifs et en introduisant la notion de pente dynamique (la friction est intégrée au terme Z), nous avons obtenu des résultats très satisfaisants avec des solveurs VFRoe et cinétique (voir par exemple [START_REF] Bourdarias | A kinetic scheme for pressurized flows in non uniform pipes[END_REF][START_REF] Bourdarias | A model for unsteady mixed flows in non uniform closed water pipes and a well-balanced finite volume scheme[END_REF][START_REF] Bourdarias | A kinetic scheme for transient mixed flows in non uniform closed pipes: a global manner to upwind all the source terms[END_REF] ou [54, Chapitre 2] et Figure 2.2a, par exemple).

originalité : traitement des points de transition et conditions aux limites.

Les points de transition entre les deux types d'écoulement sont traités via la technique des ondes fantômes, i.e l'interface est considérée comme une frontière libre, correspondant à une discontinuité du gradient de pression. Les états de part et d'autre d'une telle interface sont obtenus via la résolution d'un problème de Riemann linéaire à matrice discontinue, la ligne de discontinuité coïncidant avec la trajectoire du point de transition. La généralisation de ces travaux en présence de termes sources complexes est présentée dans [START_REF] Bourdarias | A kinetic scheme for pressurized flows in non uniform pipes[END_REF] (voir aussi [54, Chapitre 2.2] et Figure 2.2b). Nous avons proposé une nouvelle méthode pour le traitement des points de transition dans [START_REF] Bourdarias | Unsteady mixed flows in non uniform closed water pipes: a Full Kinetic Appraoch[END_REF], alternative à la méthode des ondes fantômes jusqu'alors utilisée (c. f. [START_REF] Bourdarias | A kinetic scheme for pressurized flows in non uniform pipes[END_REF]), en adéquation avec le schéma cinétique (voir Figure 2.2c). Enfin, ces études théoriques et numériques ont contribué aux développements des codes EDF-CIH Roemix (solveur VFRoe) et Flowmix (solveur cinétique).

2.2 la dynamique de l'atmosphère.

motivations.

Les équations primitives sont en pratique les équations les plus utilisées pour l'étude de la circulation atmosphérique et océanique à grande échelle et en particulier par les modèles de prévision météorologique. Elles sont issues de la dynamique des fluides géophysiques en domaine "couche mince" . Les équations primitives compressibles dérivent de l'approximation hydrostatique des équations de Navier-Stokes compressibles.

l'état de l'art.

Les premiers résultats mathématiques concernant les équations primitives à viscosité constante (celles utilisées en météorologie par exemple) ont été établis par Lions et al. [START_REF] Lions | New formulations for the primitive equations for the atmosphere and applications[END_REF] (formulation mathématique des équations et existence de solutions faibles globales en temps).

Dans la littérature, le modèle que nous appellerons équations primitives compressibles simplifiées (car les termes de diffusion de chaleur émanant du soleil et la quantité d'eau dans l'air sont négligés), est introduit la première fois par Kochin [START_REF] Kochin | On simplification of the equations of hydromechanics in the case of the general circulation of the atmosphere[END_REF] en 1936. Curieusement, les premiers résultats mathématiques sont de Gatapov et Kazhikhov [START_REF] Gatapov | Existence of a global solution to one model problem of atmosphere dynamics[END_REF] en 2005 et qui plus est, pour un modèle similaire à celui de Kochin [START_REF] Kochin | On simplification of the equations of hydromechanics in the case of the general circulation of the atmosphere[END_REF] que nous appellerons modèle intermédiaire.

originalité.

Dans un premier temps, nous nous sommes intéressés à la dérivation des équations primitives compressibles introduites par Kochin [START_REF] Kochin | On simplification of the equations of hydromechanics in the case of the general circulation of the atmosphere[END_REF]. Pour ce faire, nous avons proposé une approximation hydrostatique des équations bidimensionnelles et tridimensionnelles de Navier-Stokes compressibles en utilisant un tenseur de viscosité anisotrope non constant (voir [START_REF] Ersoy | Compressible primitive equation: formal derivation and stability of weak solutions[END_REF] et [START_REF] Ersoy | A Construction of Biorthogonal Wavelets With a Compact Operator[END_REF]Chapitre 3.2]). Cette analyse asymptotique permet de justifier les modèles en dimensions deux et trois. Le modèle en dimension deux s'écrit : Les résultats dans [START_REF] Ersoy | Existence of a global weak solution to one model of Compressible Primitive Equations[END_REF][START_REF] Ersoy | Compressible primitive equation: formal derivation and stability of weak solutions[END_REF] sont obtenus en exploitant les informations données par l'équation hydrostatique ∂ y p(ρ) = -ρg. En effet, la solution de cette équation s'écrit sous la forme (stratifiée) ρ = ξ(t, x)e -g/c 2 y où on montre que ξ est solution du modèle intermédiaire via le changement de variable z = 1e -g/c 2 y . Le modèle intermédiaire est beaucoup simple à étudier et, en dimension deux, il s'écrit :

     ∂ t ρ + ∂ x (ρu) + ∂ y (ρv) = 0, ∂ t (ρu) + ∂ x (ρu 2 ) + ∂ y (ρuv) + ∂ x p = ∂ x (ν 1 ∂ x u) + ∂ y (ν 2 ∂ y u), ∂ y p = -ρg
     ∂ t ξ + ∂ x (ξu) + ∂ z (ξw) = 0 ∂ t (ξu) + ∂ x (ξu 2 ) + ∂ z (ξuw) + ∂ x ξ = ∂ x (∂ x u) + ∂ z (∂ z u) ∂ z ξ = 0 où w(t, x, z) = e -y v
∈ L ∞ (0, T; W 1,2 (Ω)), ∂ t ρ ∈ L 2 (0, T; L 2 (Ω)) et u ∈ L 2 (0, T; W 2,2 (Ω)) ∩ W 1,2 (0, T; L 2 (Ω)), v ∈ L 2 (0, T; L 2 (Ω)).
       ∂ t h + div(q) = 0, ∂ t q + div q⊗q h + ∇ g h 2 2
= -gh∇b l'état de l'art et originalité. Nous avons proposé une nouvelle manière d'aborder ce problème en décrivant la sédimentation par une description microscopique couplée à des équations de la mécanique des fluides.

∂ t b + ξdiv(q b (h, q)) = 0 où h représente
Plus précisément, la dynamique des sédiments est décrite par l'équation de Vlasov :

∂ t f + div x (v f ) + div v ((F + g) f ) = r∆ v f
où f est la densité de particules, v ∈ R 3 la vitesse cinétique d'une particule soumise à la force de gravité g = (0, 0, -g) t et F la force de friction de Stokes. Le terme de diffusion r∆ v f décrit le mouvement Brownien des particules où r > 0 est la vitesse de diffusion donnée par la formule d'Einstein. Cette équation est couplée au système de Navier-Stokes compressible :

∂ t ρ w + div(ρ w u) = 0, ∂ t (ρ w u) + div(ρ w u ⊗ u) + ∇p(t, x) = divσ(ρ w , u) + F, où u(t, x) = (u 1 , u 2 , u 3 )(t
= k(t, x 1 , x 2 )ρ(t, x) 2 où k(t, x 1 , x 2 ) = gh(t,x 1 ,x 2 ) 4ρ f
, ρ := ρ w + ρ s est une "densité de mélange" avec ρ s la densité macroscopique de sédiments, ρ s = R 3 f dv et ρ f une densité caractéristique du fluide. Le second couplage est décrit par le terme -

R 3 F f dv présent dans F -R 3 F f dv + ρ w g.
Outre l'originalité d'introduire une loi de pression barotrope nonhomogène (voir aussi [START_REF] Bresch | Shallow-water equations and related topics[END_REF][START_REF] Bresch | Existence of global weak solutions for a 2D viscous shallow water equations and convergence to the quasi-geostrophic model[END_REF]), la réduction de ces équations vers un modèle de Saint-Venant-Exner a été obtenue par deux analyses asymptotiques successives. La première analyse asymptotique est la limite hydrodynamique de l'équation de Vlasov. En intégrant cette équation contre 1 et v, nous obtenons un système macroscopique que nous couplons avec les équations de Navier-Stokes compressibles pour aboutir à un modèle de type mélange (similaire à un modèle à deux phases). Nous procédons ensuite à la réduction de ce modèle par une analyse asymptotique en couche mince et nous obtenons les équations de Saint-Venant-Exner. Les détails de ces travaux sont disponibles dans [START_REF] Ersoy | A Construction of Biorthogonal Wavelets With a Compact Operator[END_REF]Chapitre 5].

T R AVA U X D E R E C H E R C H E : P É R I O D E P O S T D O C T O R A L E
Vous trouverez ci-dessous le résumé de mes travaux de recherche : contrôle, analyse numérique, mécanique numérique et applications, et modèles asymptotiques hydrostatique et non-hydrostatique.

analyse de sensibilité motivations.

Les applications du contrôle optimal sont nombreuses et variées, par exemple, pour le contrôle des flux routiers, ferroviaires, aériens, fluviaux, barrages EDF, etc. Dans certains de ces exemples, les équations sont des lois de conservation hyperbolique (non-linéaires) pour lesquelles l'utilisation des techniques standard basées sur la linéarisation est délicate, en raison de la présence de discontinuités. Le problème qui a été abordé ici, à notre connaissance n'a pas encore été traité. Il s'agit d'étudier l'analyse de sensibilité pour des lois de conservation scalaire (stationnaire) avec terme source. Cette analyse est pertinente dans les applications aux problèmes de contrôle optimal et d'identification des paramètres (voir Part ii). l'état de l'art et originalité.

Dans ce travail [START_REF] Ersoy | Sensitivity analysis of 1-d steady forced scalar conservation laws[END_REF], on considère un modèle simplifié de loi de conservation scalaire

∂ x f (v(x)) + v(x) = g(x), x ∈ R, (3.1)
avec des "conditions aux limites" lim x→±∞ v(x) = 0. L'objectif principal est d'effectuer une analyse de sensibilité rigoureuse des solutions v du problème (3.1) par rapport aux perturbations du terme source (forcage) g. En notant v ε = v + ε(δv) ε une solution du problème perturbé (3.1) avec g = g + εδg nous identifions la limite de la quantité

v ε -v ε = (δv) ε → δv faiblement-* dans M(R) as ε → 0, où (
δv) est l'unique solution de dualité (c. f. [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF]) du problème linéaire

∂ x f (v)δv + δv = δg.
Dans le cadre fonctionnel de Bressan and Marson [START_REF] Bressan | A variational calculus for discontinuous solutions of systems of conservation laws[END_REF], on peut aussi établir ce résultat à partir du problème non-stationnaire ∂ t u(t, x) + ∂ x f (u(t, x)) + u(t, x) = g(x), u(0, x) = u 0 (voir [START_REF] Ersoy | Sensitivity analysis of 1-d steady forced scalar conservation laws[END_REF] et Figure 3.1). 3.2 mécanique des fluides numériques et applications motivations.

La mécanique des fluides numérique (MFN) ou par le terme anglais computational fluid dynamics (CFD), a grandi d'une curiosité mathématique pour devenir un outil essentiel voire incontournable dans pratiquement toutes les branches de la dynamique des fluides et notamment en hydrodynamique côtière dont celle des vagues. Dans ce contexte, nous nous intéressons aux défis de la modélisation numérique (voir Part iii) qui consiste à simuler avec précision ces processus sur de très grandes échelles spatiales à l'aide de raffinement adaptatif de maillage (Adaptive Mesh Refinement (AMR)) basé sur des blocs (Block-Based-AMR ou BB-AMR).

l'état de l'art et originalité.

De manière générale, nous nous intéressons à l'approximation numérique par Volumes Finis de systèmes hyperboliques non linéaires multidimensionnels, écrit ici sous la forme 1d,

∂w ∂t + ∂ f (w) ∂x = 0, (t, x) ∈ R + × R w(0, x) = w 0 (x), x ∈ R. où w : R + × R → R d représente le vecteur d'état inconnu et f : R d → R d désigne le flux.
Nous avons développé un code numérique multidimensionnel pour l'approximation numérique de ces équations par Volumes Finis qui intègre un algorithme de pas de temps local, un schéma d'ordre 1 et 2 en espace-temps (MUSCL, AB2, RK2) basé sur un solveur cinétique et le solveur de Godunov dans un cadre de raffinement de maillage adaptatif [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF][START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF][START_REF] Ersoy | Adaptive scheme based on entropy production: robustness through severe test cases for hyperbolic conservation laws[END_REF][START_REF] Golay | Entropy production as mesh refinement criterion -Application to wave breaking[END_REF][START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 2: Application to tsunamis propagation[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -application to wave breaking[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -Application to wave breaking[END_REF] type décomposition de domaine (BB-AMR).

On rappelle succinctement le principe des méthodes AMR dans le cadre de la dimension 1. En notant S n k b un critère de raffinement positif sur une cellule C k b à l'instant t n du maillage, on compare ce dernier à un seuil α, par exemple la moyenne,

α = S m = 1 |Ω| ∑ k b S n k b .
On définit ensuite deux coefficients 0 < β min ≤ β max , qui déterminent le pourcentage des cellules à raffiner ou à déraffiner. Ainsi, pour chaque cellule C k b : Dans le cadre multidimensionnel, nous avons développé une méthode de type décomposition de domaine Block-Based AMR (BB-AMR) qui permet un contrôle efficace de la mémoire partagée, ce qui permet d'obtenir un temps de calcul bien mieux équilibré entre les coeurs de l'ordinateur (voir Chapter 7). Enfin, afin d'améliorer la méthode AMR, nous avons proposé une méthode permettant de calculer le seuil α = α PE automatiquement de manière à "optimiser" les zones à raffiner (jusqu'alors le seuil α était calibré manuellement). Ce dernier est construit à partir du réarrangement décroissant du critère de raffinement (voir Chapter 8). La compréhension de la dynamique des écoulements souterrains dans les zones côtières est d'un intérêt majeur dans les domaines de l'ingénierie côtière et du développement durable, notamment dans la région sud-est de la France. Une meilleure compréhension de la circulation de ces écoulements est d'une importance primordiale pour analyser et prévoir un ensemble de processus physiques et biogéochimiques des zones littorales, tels que le transport des sédiments (stabilité du lit), la diffusion des matières dissoutes tels que les polluants ou les éléments nutritifs ou encore le mélange entre les eaux continentales (douces) et marines (salées). Chacune de ces questions requiert une attention particulière dans le contexte du changement global (élévation du niveau de la mer, événements de submersion, salinisation des terres arables) et de la pression anthropique croissante (urbanisation côtière). l'état de l'art et originalité.

• si S n k b > α max = β max α, alors la cellule est raffinée et divisée en deux cellules C k b0 et C k b1 , • si S n k b0 < α min = β min α et S n k b1 < α min ,
Ces travaux s'inscrivent dans un projet de recherche (voir Part iv) visant à simuler les interactions entre les vagues et les eaux souterraines sur les plages de sable. Des études récentes en laboratoire et sur le terrain [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] ont permis de bien comprendre expérimentalement la dynamique des eaux souterraines et ont montré un schéma de circulation sous la zone de swash. Jusqu'à présent, peu de modèles ont été proposés et leur portée reste limitée. De nouveaux développements numériques sont nécessaires pour capturer un maximum de processus physiques, en particulier pour résoudre tout au long de la phase d'onde la dynamique complexe des champs de saturation et de pression. La modélisation est assez difficile à réaliser car elle nécessite de résoudre à la fois les équations non linéaires pour l'écoulement des eaux souterraines et les équations des eaux peu profondes pour l'écoulement de surface. En outre, ces modèles doivent être couplés de manière appropriée afin de respecter la physique du problème. De ce fait, les infiltrations/exfiltrations de la zone de déferlement, les fluctuations de la saturation en sable, les multiples échelles de temps et les diverses échelles spatiales nécessitent une méthodologie robuste et précise tant pour les schémas numériques que pour les algorithmes, ainsi qu'un apercu minutieux de la modélisation.

L'équation de Richards est une équation parabolique non linéaire classique pour décrire ces écoulements : 3.4 modèles asymptotiques pour les ondes longues motivations.

∂ t (θ(ψ)) -∇ • (K(ψ)∇(ψ + z)) = 0
Par opposition à l'esprit "mécanique numérique", l'approche asymptotique permet de réduire la dimension du problème mathématique pour aboutir à des équations plus simples à résoudre dont le gain en temps de calcul est significatif. C'est le cas par exemple, lorsque les écoulements sont en couche mince. On peut réduire ainsi le modèle mathématique 3d de départ en un modèle 2d, voire 1d (voir Part v).

Dans la continuité de ma thèse, j'ai justifié et développé quelques modèles réduits en couche mince (ou approximation onde longue) et notamment en introduisant le premier modèle réduit 3d vers un modèle 1d non-linéaire et faiblement dispersif pour l'hydraulique fluviale/urbaine. Les applications de ces modèles sont nombreuses et variées.

l'état de l'art et originalité : modèle hydrostatique Dans le cadre des modèles à surface libre dans une conduite fermée ou canal, par le biais de calcul asymptotique, j'ai déterminé l'expression du terme de friction en fonction du périmètre mouillé P m pour une géométrie quelconque par la formule

K(x, u) = K 0 (u) Γ b (t,x) ds A ds où u est la vitesse moyenne suivant x, K 0 (u) = C l u + C t |u| u est la friction de coefficient laminaire C l et turbulent C t , et Γ b (t,
x) désigne la frontière mouillée (i. e. la partie de la paroi en contact avec le fluide). Ainsi, la quantité A Γ b (t,x) ds représente le rayon hydraulique qui tient compte de la géométrie de la section. Un résultat similaire a été obtenu pour le modèle en charge. Je vous invite à consulter ces articles [START_REF] Ersoy | Dimension reduction for incompressible pipe and open channel flow including friction[END_REF][START_REF] Ersoy | Dimension reduction for compressible pipe flows including friction[END_REF] qui ne sont pas détaillés dans ce mémoire. l'état de l'art et originalité : modèle hydrostatique intégrant des termes de pluie et d'infiltration L'un des modèles les plus utilisés pour décrire l'hydraulique pluviale est le modèle de Saint-Venant. Nous proposons un nouveau modèle permettant non seulement de modéliser ce phénomène mais aussi de prendre en considération l'influence de la trainée générée par l'apport de la pluie ou par l'infiltration. Le modèle obtenu s'écrit comme suit

∂ t h + ∂ x q = S := R -I, ∂ t q + ∂ x q 2 h + g h 2 2 = -gh∂ x Z + S q h -k + (R) + k -I + k 0 q h q h
où les inconnues h(t, x) and q(t, x) = h(t, x)u(t, x) désignent respectivement la hauteur d'eau et la vitesse moyenne de l'écoulement, g étant la gravité et Z représente la topographie. Le terme source S représente l'apport par les eaux pluviales R ou la perte d'eau par l'infiltration I. Enfin les termes de la forme k p représentent la trainée générée par l'apport de la pluie p = + et l'infiltration p =et la topographie p = 0. Je vous invite à consulter cet article [START_REF] Ersoy | A Saint-Venant shallow water model for overland flows with precipitation and recharge[END_REF] qui n'est pas détaillé dans ce mémoire.

l'état de l'art et originalité : extension au cas non-hydrostatique La modélisation de l'hydrologie des bassins versants et des rivières occupe une place centrale dans les sciences de l'environnement, notamment en ce qui concerne la disponibilité de l'eau, les réseaux d'égouts urbains, les risques d'inondation, etc.

L'un des modèles les plus utilisés pour décrire le mouvement des cours d'eau est le modèle à surface libre moyenné par section (voir [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF] ou [54, Chapitre 1]), à savoir, les modèles de type Free Surface model (FS-model). Cependant, pour certains régimes, ces modèles sont incapables de reproduire les trains d'ondes ou "ondes de choc dispersives" qui sont induits par une distribution de pression non hydrostatique. En procédant à un développement asymptotique à l'ordre deux, on peut ainsi construire un nouveau modèle non-linéaire et faiblement dispersif [START_REF] Debyaoui | A Generalised Serre-Green-Naghdi equations for variable rectangular open channel hydraulics and its finite volume approximation[END_REF][START_REF] Debyaoui | Generalised Serre-Green-Naghdi equations for open channel and for natural river hydraulics[END_REF] (voir Part v) qui généralise le modèle FS-model et les équations de Serre-Green-Naghdi :

         ∂ t A + ∂ x Q = 0 ∂ t Q + ∂ x Q 2 A + I 1 (x, A) + µ 2 ∂ x (G(x, A)D(u)) = I 2 (x, A) + µ 2 G(x, A, Q) + O(µ 2 2 )
où A est l'aire mouillée, Q = Au est le débit moyen, u est la vitesse moyenne par section, I 1 (resp. I 2 ) est la pression hydrostatique (resp. terme source hydrostatique) , G(x, A) généralise le coefficient h 3 3 dans les équations de Serre-Green-Naghdi,

D(v) = (∂ x v) 2 -∂ x ∂ t v -v∂ xx v et G(x, A, Q) un

Part II S E N S I T I V I T Y A N A LY S I S

Dans la nature, tout a toujours une raison.

Si tu comprends cette raison, tu n'as plus besoin de l'expérience.

-

Léonard de Vinci

In this part, in collaboration with E. Feireisl and E. Zuazua, we analyze a 1d forced steady state scalar conservation laws [START_REF] Ersoy | Sensitivity analysis of 1-d steady forced scalar conservation laws[END_REF].

In Chapter 4, we first show the existence and uniqueness of entropy solutions as limits as t → ∞ of the corresponding solutions of the scalar evolutionary hyperbolic conservation law.

In Chapter 5, we then linearise the steady-state equation with respect to perturbations of the forcing term. This leads to a linear first order differential equation with, possibly, discontinuous coefficients. We show the existence and uniqueness of solutions in the context of duality solutions. We also show that this system corresponds to the steady state version of the linearised evolutionary hyperbolic conservation law. This analysis leads us to the study of the sensitivity of the shock location with respect to variations of the forcing term, an issue that is relevant in applications to optimal control and parameter identification problems.

Possible extensions and further discussions can be found in Chapter 6.

Optimal control of solutions to (non-linear) hyperbolic conservation laws is hampered by the presence of discontinuities (shock waves) making the use of the standard techniques based on linearization rather delicate. The same can be said about optimal design and parameter identification problems. For an analysis of these issues in the context of evolutionary hyperbolic conservation laws in 1d we refer to [START_REF] Castro | An alternating descent method for the optimal control of the inviscid Burgers equation in the presence of shocks[END_REF] and [START_REF] Castro | Flux identification for 1d scalar conservation laws in the presence of shocks[END_REF] and the references therein. The issue of the convergence of timeevolution controls towards the steady state ones as the time horizon tends to infinity is a subject that recently has attracted attention (see, for instance, [START_REF] Cardaliaguet | Long time average of mean field games[END_REF] and [START_REF] Porretta | Long time versus steady state optimal control[END_REF]). But, as far as we know, this has not yet been addressed in the context of the scalar conservation laws considered here.

In this work, we consider a simplified model problem based on the 1d scalar steady driven conservation law

∂ x f (v(x)) + v(x) = g(x), x ∈ R, (4.1)
supplemented with a "far field" boundary condition

lim x→±∞ v(x) = 0.
Note that, as we shall see in Section 4.1, v can be viewed as the asymptotic limit for t → ∞ of solutions to the associated evolutionary scalar hyperbolic conservation law

∂ t u(t, x) + ∂ x f (u(t, x)) + u(t, x) = g(x), u(0, x) = u 0 . (4.2)
Our main goal is to perform a rigorous sensitivity analysis of solutions v of the steady problem (4.1) with respect to perturbations of the forcing term g. In particular, denoting v ε the solution of (4.1) with g = g + εδg we identify the limit of the quantity

v ε -v ε = h ε → h as ε → 0,
where v is the solution associated to g. It turns out that h is in general a measure on R, with a singular part sitting on the set of discontinuities of v, taking account of the sensitivity of the shock location of v, see Chapter 5. Our approach is based on an adaptation of the concept of duality solutions introduced by Bouchut and James [START_REF] Bouchut | One-dimensional transport equations with discontinuous coefficients[END_REF][START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF] and Strömberg [START_REF] Strömberg | On viscosity solutions of irregular Hamilton-Jacobi equations[END_REF]. introduction 4.1 the asymptotic limit as t → ∞ This chapter is devoted to analysis of the models under consideration from close, using the notion of entropy solutions. Of course, our analysis can be viewed as a particular instance of the general theory of nonlinear semigroups in L 1 (R) (see [START_REF] Crandall | Generation of semi-groups of nonlinear transformations on general Banach spaces[END_REF]). The steady problem and the evolution one are linked in the sense that the later can be solved in an unique manner by the semigroup of L 1 (R)-contractions generated by the accretive operator associated to the steady state problem.

Throughout this work we suppose that the flux function f is locally Lipschitz,

f ∈ W 1,∞ loc (R).
We start with the nowadays standard definition of entropy solution to the evolutionary problem (4.2) introduced by Kružkov [START_REF] Kruzhkov | Generalized solutions of the Cauchy problem in the large for first order nonlinear equations[END_REF], [START_REF] Kruzhkov | Results on the nature of the continuity of solutions of parabolic equations, and certain applications thereof[END_REF] and [START_REF] Kružkov | First order quasilinear equations in several independent variables[END_REF].

Definition 4.1.1. Let g ∈ L ∞ (R). A function u = u(t, x), u ∈ L ∞ ((0, T) × R) ∩ C([0, T]; L 1 loc (R)),
is an entropy solution to problem (4.2) if it solves the equation in the sense of distributions, and, moreover, fulfils the integral inequality

R E (u(τ, •)) ϕ(τ, •)dx - R E u 0 ϕ(0, •) dx ≤ τ 0 R E(u)∂ t ϕ + F(u)∂ x ϕ dx dt + τ 0 R -E (u)uϕ + E (u)gϕ dx dt, (4.3
)

for any τ ∈ [0, T], any test function ϕ ∈ C ∞ c ([0, T] × R), ϕ ≥ 0, and any convex entropy E : R → R with F ≡ E f .
Note that the existence of global-in-time entropy solutions for u 0 ∈ L ∞ (R) can be established by means of artificial viscosity approximations (see Kružkov [START_REF] Kružkov | First order quasilinear equations in several independent variables[END_REF]).

It can also be shown that two entropy solutions u 1 , u 2 emanating from the initial data u 1 0 , u 2 0 satisfy

|x|≤M u 1 (τ, •) -u 2 (τ, •) dx + τ 0 |x|≤M+λ(τ-t) u 1 -u 2 dx dt ≤ |x|≤M+λτ u 1 0 -u 2 0 dx
for any τ ≥ 0, M > 0, where

λ = sup | f (z)| |z| ≤ (1 + τ) max u 1 0 L ∞ (R) , u 2 0 L ∞ (R) , g L ∞ (R) . Consequently, R u 1 (τ, •) -u 2 (τ, •) dx + τ 0 R u 1 -u 2 dx dt ≤ R u 1 0 -u 2 0 dx.
In particular, the entropy solutions are uniquely determined by the initial data.

As an immediate consequence of the previous estimate we deduce the exponential stability of entropy steady state solutions.

, •) = u 0 ∈ L ∞ (R), such that u 0 -v L 1 (R) < ∞. Then u(t, •) -v L 1 (R) ≤ exp(-t) u 0 -v L 1 (R) for any t ≥ 0.
Remark 4.1.1. In Corollary 4.1.1 the stationary solution v is arbitrary, in particular, it does not need to satisfy the far field condition provided u 0v belongs to L 1 (R).

The conclusion of Corollary 4.1.1 remains valid for solutions of general scalar conservation laws in R N provided the flux field is locally Lipschitz.

stationary solutions

For the sake of simplicity, we focus on the stationary solutions satisfying the far field condition. To this end, we restrict ourselves to the class of forcing terms

g ∈ L 1 ∩ L ∞ (R).
Under these circumstances, in view of (4.1), it is easy to see that the stationary problem (4.1) admits a unique entropy solution v in the same class, namely,

v ∈ L 1 ∩ L ∞ (R).
In addition, by the conservation of mass we have

R v 1 -v 2 dx = R g 1 -g 2 dx, the comparison principle R v 1 -v 2 + dx ≤ R g 1 -g 2 + dx and the L 1 -contraction property R v 1 -v 2 dx ≤ R g 1 -g 2 dx. (4.4)
where v 1 , v 2 are the entropy solutions of (4.1) with g = g 1 , g = g 2 , respectively. The existence of stationary solutions out of (4.1) can be proved in several different ways. For instance, using the L 1 -distance as Lyapunov function and LaSalle's invariance principle (see [START_REF] Zuazua | A dynamical system approach to the self-similar large time behavior in scalar convection-diffusion equations[END_REF]) one can show that there is an unique stationary solution and that all other solutions converge exponentially as t → ∞ to it. One can also construct the stationary solutions as limits of time periodic solutions of time-period τ > 0 with τ → 0. These periodic solutions can be built as fixed points (using Banach contraction principle) of the semigroup map associating the value of the solution at time t = τ to the initial datum. Furthermore, one can also build the stationary solutions as vanishing viscosity limits of elliptic equations. In fact, as mentioned above, these results are also a particular example of the classical theory of nonlinear contraction semigroups in L 1 (R).

The main result of this chapter is as follows:

Proposition 4.2.1. Let the flux function f be continuously differentiable on R and non-degenerate in the sense that the critical points f (y

) = 0 are isolated in R. Let g ∈ L 1 ∩ L ∞ ∩ BV(R).
Then the stationary problem (4.1) possesses an entropy solution v determined uniquely in the class

v ∈ L 1 ∩ L ∞ (R)
such that v is continuous in R with a possible exception of a countable set of points {s i }, with

lim x→s i - v(x) = v i -= v i + = lim x→s i + v(x), i = 1, . . . , N.
Each open interval with end points v i -, v i + contains at least one critical point y i of f such that either

f (y i ) < f (v i -) = f (v i + ) yielding v i -> v + i or f (y i ) > f (v i -) = f (v i + ) yielding v i -< v + i .
Proof: Assuming, in addition to (4.2), that

g ∈ BV(R), (4.4) yields immediately v ∈ BV(R).
Identifying v with its Lebesgue means,

v(x) = lim h→0 1 2h x+h x-h v(z) dz,
we deduce that v is continuous with a possible exception of countably many points at which the left and right limits exist. Assume that

x 0 ∈ R is a point of discontinuity of v, specifically, lim x→x 0 - v(x) = v -= v + = lim x→x 0 + v(x),
where, as

f (v) is Lipschitz continuous (since ∂ x ( f (v)) = g -v), f (v-) = f (v + ).
Note that this null jump condition is the natural limit of the Rankine-Hugoniot condition for the damped evolutionary hyperbolic conservation law (4.2) or, in other words, the condition characterizing stationary shocks.

Since v is an entropy solution of (4.1) we deduce that

∂ x [( f (v) -f (k))sgn(v -k)] ≤ c in the sense of distributions,
for some finite c. This can be easily seen by the vanishing viscosity argument, or directly from (4.3) in view of the fact that both g and v are bounded, for instance. In particular,

lim x→x 0 - ( f (v) -f (k)) sgn(v -k) = ( f (v -) -f (k)) sgn(v --k) ≥ ( f (v + ) -f (k)) sgn(v + -k) = lim x→x 0 + ( f (v) -f (k)) sgn(v -k)
provided k belongs to the open interval I with the end points v -, v + . We distinguish three complementary cases:

• f (k) < f (v -) = f (v + )
for some k in the interval linking v -and v + . In this case, relation (4.2) implies

v -> v + . • f (k) > f (v -) = f (v + )
for some k in the interval linking v -and v + . Similarly, we deduce from (4.2) that

v + > v -. • The degenerate case in which f (v-) = f (k) = f (v +
) can be excluded by the main assumptions of the Proposition in which we impose the set of critical points of f to be isolated.

S E N S I T I V I T Y W I T H R E S P E C T T O T H E F O R C I N G T E R M
In this chapter, we study the sensitivity of solutions to the stationary problem (4.1), supplemented with the far field condition, with respect to perturbations of the right-hand side g. To this end, we introduce the concept of duality solutions in the spirit of Bouchut and James [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF].

duality solutions for the stationary problem

We start with a prototype example of a scalar conservation law, where the flux function f is strictly convex with a (global) minimum attained in R. By virtue of Proposition 4.2.1, solutions of problem (4.1) admit a countable (possibly empty) set of singularities where the solutions "jump down" across the shock, meaning they satisfy the so-called Oleinik condition. Note that, in particular, the degenerate case can be excluded because of the strict convexity assumption on f . Furthermore, by comparison (the maximum principle), it can be shown that

∂ x v ≤ g in D (R) (5.1)
provided that

∂ x g ≤ g in D (R). (5.2)
We now study the effect of a perturbation

g ε = g + ε(δg), with δg ∈ L 1 ∩ L ∞ ∩ BV(R) δg of compact support in R, ∂ x (δg) ≤ δg in D (R).
This leads to an entropy solution v ε ∈ L 1 ∩ L ∞ (R) of (4.1) that can be written in the form

v ε = v + ε(δv) ε ,
sensitivity with respect to the forcing term more specifically,

∂ x f ([v + ε(δv) ε ]) + [v + ε(δv) ε ] = g + ε(δg).
Note that, since we deal with entropy solutions,

(δv) ε L 1 (R) ≤ δg L 1 (R) .
Furthermore, we have

∂ x 1 ε ε 0 f (v + z(δv) ε ) dz (δv) ε + (δv) ε = δg;
whence the perturbation h = (δv) ε solves the linear problem:

∂ x (A ε h ε ) + h ε = δg, where A = 1 ε ε 0 f (v + z(δv) ε ) dz , A ε L ∞ (R) ≤ c g L ∞ (R) , (δg) L ∞ (R) , A ε → A weakly-* in L ∞ (R).
Furthermore, in accordance with (5.1), (5.2),

∂ x A ε = 1 ε ε 0 f (v + z(δv) ε ) ∂ x (v + z(δv) ε ) dz ≤ c( f )(g + δg), ( 5.3)
with

∂ x g ≤ g, ∂ x (δg) ≤ δg.
This yields local uniform BV-bounds on A . These arguments, together with the continuity of f allow us to deduce that

A = f (v),
As a consequence we expect the limit measure h to be a solution of

∂ x f (v)h + h = δg. (5.4)
Note however that f (v) may be discontinuous at the shock discontinuities of v. Thus, we need to introduce a suitable concept of weak solution in the class of measures in order to determine the solution h in a unique manner and to justify the limit process above. For doing that, the fact that the coefficients involved in the linearised equation fulfill a one-sided Lipschitz condition of the form (5.3) will play a key role.

Definition 5.1.1. We say that h is a duality solution of

∂ x (Ah) + h = δg (5.5) if d dt R ph dx + R ph dx = - R (δg)p dx (5.6)
for any reversible solution p of the evolutionary problem

∂ t p + A∂ x p = 0 (5.7)
in the sense of Bouchut and James [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF].

Remark 5.1.1.

• In [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF] reversible solutions of the adjoint equation (5.7) are built as limits of approximating sequences of solutions of a final value problem at t = T for regularized potentials A. In [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF] it is shown that, under the one-sided Lipschitz condition on A, this leads to a well identified unique solution, the so-called reversible one.

• Typically, as shown in [START_REF] Castro | An alternating descent method for the optimal control of the inviscid Burgers equation in the presence of shocks[END_REF], when linearizing scalar conservation laws along a solution with a shock discontinuity, the adjoint system of the linearised one has the structure (5.7) with the potential A being discontinuous along the shock. The corresponding solution p can be defined by characteristics away from the zone of influence of the shock. The reversible solution is that taking over this set the corresponding value at the shock location at the final time t = T. In this way, for locally Lipschitz continuous data at time t = T reversible solutions are locally Lipschitz.

Accordingly, the duality formula (5.6) makes sense for h a measure.

• Note that duality solutions of (5.5) are uniquely determined by the right-hand side δg in the class of measures h ∈ M(R).

Seeing that, by virtue of (5.3), the partial derivatives ∂ x A ε are bounded from above so that they fulfill uniformly the one-sided Lipschitz condition. Then we may use the abstract convergence result of Bouchut and James [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF] to conclude that

h ε ≡ (δv) ε → δv weakly-* in M(R),
where δv ∈ M(R) is the unique duality solution of the linearised problem

∂ x f (v)δv + δv = δg.
Note that the compactness results in [START_REF] Bouchut | Differentiability with respect to initial data for a scalar conservation law[END_REF] refer to the evolution equation but, accordingly, they can also be applied to steady state solutions. This procedure yields a duality solution to the limit linearised steady-state problem (5.4).

The uniqueness of the duality solution for the steady problem (5.4) is easy to prove. In case there were two distinct solutions, h 1 and h 2 ,

then h = h 1 -h 2 would satisfy d dt R ph dx + R ph dx = 0
for all reversible solution of (5.7). Considering all possible reversible solutions of (5.7) associated to all Lipschitz continuous data p T = p T (x) at t = T we would conclude that R p T h dx = 0, and this would yield h ≡ 0.

We have proved the following result.

Theorem 5.1.1. Let f : R → R be a strictly convex function in C 2 (R). Let g, δg ∈ L 1 ∩ L ∞ ∩ BV(R)
be given such that

∂ x g ≤ g, ∂ x (δg) ≤ δg in D (R). Let v ε = v + ε(δv) ε ∈ L 1 ∩ L ∞ (R)
be the solution of the perturbed problem (5.1), where v

∈ L 1 ∩ L ∞ (R) is the unique solution of (4.1). Then (δv) ε → δv weakly-* in M(R) as ε → 0,
where (δv) is the unique duality solution of the linear problem

∂ x f (v)δv + δv = δg (5.8)
in the sense specified in Definition 5.1.1.

Remark 5.1.2. In accordance with [16, Theorem 2.2], the function f (v) in (5.8) may be redefined on the set of zero measure in such a way that the resulting function is Borel and (5.8) holds in the sense of distributions. In particular, the singular part of the measure δv is supported by shocks of v.

sensitivity of shocks

As in the previous chapter, we assume that the flux function f is strictly convex with a (global) minimum attained in R. Actually, to simplify the presentation we assume that f is even, as it is for instance the case for the Burgers equation where f (u) = u 2 /2. Assume that an entropy solution v exhibits a shock discontinuity at x = 0 ∈ supp g. Our goal is to justify the following relation:

δϕ = lim ε→0 ϕ ε ε = [ f (v)δv] 0 [v] 0 (5.9)
where [u] x 0 := u(x + 0 )u(x - 0 ) stands for the jump of a function u across x = x 0 , ϕ ε stands for the shock location of the perturbed solution v ε = v + εδv ε under forcing g ≡ g + εδg ε . In particular, we have lim

x→0 - v(x) = v(0 -) > v(0 + ) = lim x→0 + v(x) = -v(0 -) . Note that v(0 + ) = lim x→0 + v(x) = -v(0 -) since f is even and f (v(0 + )) = f (v(0 -)).
The function x → f (v(x)) is Lipschitz continuous (note that the Rankine-Hugoniot condition ensures in this case the continuity of f (v) across the shock), in particular, there exists a < 0 and b > 0 such that the solution has the following structure

v(a) = 0, v(x) > 0 in (a, 0], v ∈ C ([a, 0]) , and v(b) = 0, v(x) < 0 in [0, b), v ∈ C ([0, b]) .
Accordingly, the solution may not have an infinite number of shocks accumulating at a point and shocks are isolated. Moreover, if g is as in Theorem 5.1.1, the function f (v) satisfies the one-sided Lipschitz condition (5.3) and the linearised equation

∂ x ( f (v)δv) + δv = δg, δg ∈ L 1 ∩ L ∞ , admits a unique solution δv ∈ L 1 (a, 0) and δv ∈ L 1 (0, b).
Consider now the solutions v ε = v + εδv ε of the perturbed problem

∂ x ( f (v ε )) + v ε = g + εδg .
Since v ε , v are entropy solutions, we get

||v ε -v|| L 1 (R) ε ||δg|| L 1 (R)
and, in particular,

δv ε → δv weakly-* in M(R) .
In accordance with Theorem 5.1.1, the limit measure δv is unique, the so-called duality solution of the limiting stationary problem. We have

f (v ε ) → f (v) in C ([a, b]) ,
and we deduce that v ε possesses a unique shock discontinuity at the

point x = ϕ ε in any compact interval [a + δ, b -δ] for any δ > 0 provided ε = ε(δ) is small enough.
Moreover, ϕ ε → 0 as ε → 0 and one can quantify the variation of the shock location with respect to the perturbation g → g + εδg. Indeed, I ε being an open interval with the end points 0, ϕ ε . We have

||δv ε -δv|| C([a+δ,b-δ]\I ε ) → 0 as ε → 0.
On the other hand,

I ε δv ε dx = 1 ε I ε (v ε -v) dx = ϕ ε ε 1 |I ε | ε |v ε -v| dx,
where

I ε δv ε dx → f (v)δv (0+) -f (v)δv (0-) while 1 |I ε | ε |v ε -v| dx → [v] (0-) -[v] (0+).
Thus we have proved (5.9).

stationary problem seen as an evolutionary one

In this chapter, we still assume that the flux function f is strictly convex with a (global) minimum attained in R.

The sensitivity analysis of the steady equation (4.1) with respect to forcing can be also obtained as a formal limit as t → ∞ of the unsteady one (4.2) in view of the exponential stability of those solutions (see Corollary 4.1.1).

First of all we consider the evolution problem (4.2) and, for the sake of simplicity, we assume that the solution develops at most one shock. Let us assume that u(t, x) is a Lipschitz continuous solution of equation (4.2) on Ω ± separated by a regular curve

Σ = {(t, ϕ(t)), t > 0},
where it satisfies the Rankine-Hugoniot condition

ϕ (t) = [ f (u)] ϕ(t) [u] ϕ(t)
and the entropy condition [u] ϕ(t) 0.

In the presence of shocks, the state of equation (4.2) needs to be viewed as a pair (u(t, •), ϕ(•)) and problem (4.2) can be written as follows:

               ∂ t u + ∂ x ( f (u)) + u = g(x), t 0, x ∈ Ω(t), ϕ (t) = [ f (u)] ϕ(t) [u] ϕ(t) , t ∈ (0, T), u(0, x) = u 0 (x), x ∈ Ω(0), ϕ(0) = ϕ 0
where Ω(t) stands for the union of the two regions to both sides of the shock:

Ω(t) = {x < ϕ(t)} ∪ {x > ϕ(t)}.
We then analyse the sensitivity of the pair (u, ϕ) with respect to perturbations δg of the source term g. The functional framework adopted here is based on the generalized tangent vectors introduced by Bressan and Marson [START_REF] Bressan | A variational calculus for discontinuous solutions of systems of conservation laws[END_REF] (see also [START_REF] Castro | An alternating descent method for the optimal control of the inviscid Burgers equation in the presence of shocks[END_REF]).

Let

g ∈ L 1 ∩ L ∞ ∩ BV(R) with supp g compact in R
be the source term. Let

u 0 ∈ L 1 ∩ L ∞ (R) with supp u 0 compact in R
be the initial datum with a single discontinuity at

x = ϕ 0 ∈ supp g
and consider a generalized tangent vector

(δg, δu 0 , δϕ 0 ) ∈ (L 1 ∩ L ∞ ∩ BV(R)) × (L 1 ∩ L ∞ (R)) × R.
Then, for ε small enough, the perturbed solution

u ε (t, •) = u(t, •) + εδu(t, •) of (4.2) with g = g ε (≡ g + εδg) is Lipschitz continuous devel- oping a single discontinuity at x = ϕ ε (t) for t > 0.
Thus, u ε (t, •) generates a generalized tangent vector

(δu(t, •), δϕ(t)) ∈ (L 1 ∩ L ∞ (R)) × R.
One can prove (see for instance [START_REF] Bressan | A maximum principle for optimally controlled systems of conservation laws[END_REF]) that it solves the following linearised system:

                 ∂ t δu + ∂ x ( f (u)δu) + δu = δg, (t, x) ∈ (0, T) × Ω(t), δϕ (t) [u] ϕ(t) + δϕ(t) ϕ (t) [u x ] ϕ(t) -[ f (u)u x ] ϕ(t) +ϕ (t) [δu] ϕ(t) -[ f (u)δu] ϕ(t) = 0, t ∈ (0, T), δu 0 (x) = δu(0, x), x ∈ Ω(0), δϕ(0) = δϕ 0 .
(5.10)

This system has a unique solution which can be computed in two steps. Characteristics of (5.10) and (4.2) being the same, providing u solution of (4.2), the method of characteristics determines δu outside the shock curve. Thus the values of u and u x to both sides of the shock curve allow to determine the jump relation through x = ϕ(t). As a consequence all coefficients of the ODE are known and δϕ is then obtained by solving the second equation in (5.10).

Taking into account the exponential stability of u (see Corollary 4.1.1), the linear structure of equation (5.10), and passing formally to the limit as t → ∞ we get the following linearised steady equation

∂ x ( f (v)δv) + δv = δg, x ∈ R/{ϕ}, δϕ [ f (v)v x ] ϕ + [ f (v)δv] ϕ = 0.
Note that this is a formal argument. To make it rigorous one should carefully prove that the limit as t → ∞ of all terms actually exists what requires uniform regularity estimates on solutions.

In particular, formally, the second equation is obtained as limit as t → ∞ of the linearised evolution shock condition in (5.10). Indeed, using, formally, that both ϕ (t) and δϕ (t) tend to zero as t → ∞ we get the linearised shock condition in (5.3).

Moreover it provides the sensitivity of the shock location (which can be also obtained easily by linearizing the static shock condition [ f (v)] ϕ = 0) and can be written as

δϕ = [ f (v)δv] ϕ [v] ϕ
(see also equation (5.9)).

Let us illustrate the sensitivity of the shock location with respect to small perturbations of the source term for T large enough. We consider in Figure 5.1 g ε (x) = g(x) + εδg(x) where

δg(x) = 0 if 0 x π g(x) if π < x 2π
as a perturbation of g(x) = sin(x). In Figure 5.1a we show the behavior of the solution v ε under the perturbation g ε and in Figure 5.1b the sensitivity of the shock location. In particular, according to numerical simulations the shock seems to disappear for ε 8.6. In this chapter, throughout several examples, we describe some qualitative properties of steady solutions such as their support, shock location, etc. We also discuss a property that numerical simulations seem to indicate, according to which, solutions, at least in part, stabilize in finite time. This constitutes an interesting open problem.

For the sake of simplicity, let us again consider the Burgers equation

(4.1) with f (v) = v 2 /2 and g ∈ L 1 (R) ∩ L ∞ (R) of compact support.

structure of steady solutions

In what follows, given a finite number of points {a 0 , . . . , a N } we consider g as the sum of

g(x) = N-1 ∑ i=0 g i (x)1 (a i ,a i+1 ) (x) (6.1)
where g i (x) = g(x) |(a i ,a i+1 ) is either positive or negative and so that sign changes from a subinterval to the other. Following Mascia and Terracina [START_REF] Mascia | Large-time behavior for conservation laws with source in a bounded domain[END_REF], under such circumstances, one can easily prove the following: Properties 6.1.1. For a given 1 i N, the stationary solution has at most one shock s ∈ [a i-1 , a i ].

Proof. Assume by contradiction that there exist two shocks, say, (s 1 , s 2 ) ∈ [a i-1 , a i ] 2 with s 1 < s 2 . Then there exist three C 1 solutions of problem (4.1) v -, v m and v + such that the solution on the whole interval

[a i-1 , a i ] is: v(x) = v -(x)1 (a i-1 x<s 1 ) (x) +v m (x)1 (s 1 <x<s 2 ) (x) + v + (x)1 (s 2 <x a i ) (x)
where

v -(s - 1 ) > 0 > v m (s + 1 ) and v m (s - 2 ) > 0 > v + (s + 2 ). On one hand, since we have v m (s + 1 ) < 0 < v m (s - 2 ) then there exists s 1 < s < s 2 such that v m (s) = 0. On the other, v m (x) satisfies ∀x ∈ (s 1 , s 2 ), f (v m (x))v m (x) = g(x) -v m (x) with f (v m (s)) = f (0) = 0.
As a consequence, we have g(s) = 0 which is a contradiction with the decomposition (6.1) of g.

In what follows, we illustrate examples of steady solutions of equation (4.1) for f (v) = v 2 /2 when considering continuous source terms (see Figure 6.1a and Figure 6.2-6.3),

g c (x) = sin(4x)1 (0,2π) (x)
and discontinuous source terms (see Figure 6.1b),

g dc (x) = 1 (0,1) (x) -0.8 1 (1,2) (x) + 0.4 1 (2,3) (x) -0.5 1 (3,4) (x) + 0.7 1 (4,5) (x) -2 1 (5,6) (x) + 2 1 (6,7) (x) -1 (7,8) (x) + 3 1 (8,9) (x).
The results are displayed on Figure 6.1 and, for both cases, there is only one shock per interval where g keeps the same sign.

Note also that, while in the first example, the support of the steady solution v is included in the one of the forcing term g c , that is not the case in the second case. Indeed, in that one the support of v goes beyond that of g dc to the right. Over there the solution exhibits a simple linear behavior in agreement with the equation fulfilled beyond the support of the forcing term (∂ x (v 2 /2) + v = 0).

evolution problem

Let us come back to the two previous examples. The inspection of the time evolution of several numerical simulations1 seems to indicate that there exists one part of the solution (the one localized on the support of g) that stabilizes in finite time (this can be observed, for instance in Figure 6.2), while the other one decays exponentially to zero (as displayed on Figure 6.3) beyond the support of g where the equation reads

∂ t u + ∂ x f (u) + u = 0.
More precisely, it seems that the following assertion holds: Conjecture: For any u 0 and g with compact support there exists a finite time

T * > 0 such that ∀t T * , ∀x ∈ supp g, u(t, x) = v(x),
where v is the unique entropy steady state solution.

Note that rigorous analysis leads to the exponential convergence of the total mass but that the behavior inside the support of the forcing term g requires further work. indeed, let m(t) = R u(t, x) dx be the mass of u over R. It satisfies

m + m = R gdx.
Thus, one has

m(t) = R g(x) dx + m(0) -R g(x) dx e -t = R v(x) dx + m(0) -R g(x) dx e -t since R v(x) dx = R g(x) dx.
For a better understanding of the large time behavior and the possible stabilization in finite time a finer analysis of the characteristic of Equation (4.2). The numerical simulations seem to indicate the following:

• In Figure 6.2, one has m(0) = R g(x) dx = 0 and thus m(t) = 0, for all t > 0. Moreover it seems that supp u ⊂ supp g, for all t > 0. Such a solution should stabilize in finite time.

• In Figure 6.3, one has m(0) = 0, R g(x) dx = 0 and thus m(t) = (1e -t ) R g(x) dx, for all t > 0. Part of u(t) escapes supp g, then it seems to converge to zero exponentially but not in finite time. It seems that, nevertheless, the part of the solution restricted to supp g reaches the steady state in finite time.

As we said above these global or local finite time stabilization properties constitute interesting open problems to be rigorously analyzed. -Isaac Newton

In this part, I present the work in collaboration with T. In Chapter 7, we present the finite volume method developed for non-linear hyperbolic systems in adaptive mesh refinement framework. The mesh refinement algorithm is driven by

• a mesh refinement indicator and,

• a mesh refinement threshold.

A mesh refinement indicator is based on an error estimation procedure or a feature detection technique. In practice, we use the well-known numerical density of entropy production as a mesh refinement indicator and we refer to [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF][START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF][START_REF] Ersoy | Adaptive scheme based on entropy production: robustness through severe test cases for hyperbolic conservation laws[END_REF][START_REF] Golay | Entropy production as mesh refinement criterion -Application to wave breaking[END_REF][START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 2: Application to tsunamis propagation[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -application to wave breaking[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -Application to wave breaking[END_REF] and the reference therein, where we show its efficiency through several numerical test cases in 1-d, 2-d and 3d cases.

As a well-known fact, the mesh refinement threshold is a problem-dependent parameter whose choice might be the critical weakness of mesh adaptation methods. This issue is addressed and an automatic selection of the threshold is proposed in Chapter 8 for which one can show that whatever the mesh refinement indicator is, the obtained results is almost of the same order of accuracy.

CFD has grown from a mathematical curiosity to become an essential tool in all branches of fluid dynamics, particularly in coastal hydrodynamics, including wave hydrodynamics. Boosted by the continuous improvement of computer technology, CFD simulations are expected to provide a unique insight, in terms of spatio-temporal resolution and controlled conditions, into such wave dynamics. Therefore, the numerical simulation of waves is both an attractive research topic and a challenging task for coastal, environmental and hydraulic engineering.

The involved physical processes are quite complex and the scales are often large in general. A major issue of many modelling challenges is thus to accurately simulate those processes over very large ranges in spatial scale inexorably leading to heavy computational time on uniform meshes.

In principle, AMR [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] or Adaptive Moving Mesh (AMM) [START_REF] Arpaia | r-adaptation for Shallow Water flows: conservation, well balancedness, efficiency[END_REF][START_REF] Khakimzyanov | Numerical simulation of conservation laws with moving grid nodes: Application to tsunami wave modelling[END_REF] methods allow to solve in a reasonable Central Processing Unit (CPU) time these equations adjusting the computational effort locally to maintain a quite uniform level of accuracy. The first one relies on macro cells which can be refined (and then possibly coarsened) while the other method allows all mesh nodes to move to generate big or small cells (with a fixed nodes number for all the simulation). The zones where the mesh needs to be moved or refined is determined thanks to a mesh refinement criterion. Depending on the application, equations and numerical methods, a variety of different criteria might be used based on an error estimation procedure or a feature detection technique. Basically, the zones where a given threshold is exceeded are moved or refined. In particular, for AMR methods, the more the threshold is small, the more accurate is the result at the expense of the computational time. In this work, we focus on h-AMR method.

governing equations and finite volume approximation

In the most of the applications concerned in this work, wave dynamics are governed by non-linear hyperbolic systems:

∂w(t) ∂t + ∇ • f (t, w) = G, (t, x) ∈ R + × R d w(0, x) = w 0 (x), x ∈ R d
where w, f , G stands respectively for conservative variables, flux and source.

In what follows, for the sake of simplicity, we present the numerical algorithms in 1d. We refer to the corresponding papers for the multidimensional case (see for instance [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF]).

We recall here the well-known construction of a numerical approximation of the following general non-linear hyperbolic conservation laws

∂w ∂t + ∂ f (w) ∂x = 0, (t, x) ∈ R + × Ω w(0, x) = w 0 (x), x ∈ R. (7.1)
where Ω ⊂ R, w ∈ R d denotes the vector state (referred as the conservative variables) and f denotes the flux governing the physical description of the flow.

The computational domain is split into control volumes

C k = ]x k-1/2 , x k+1/2 [ of mesh size h k with x k±1/2 = x k ± h k /2.
The unknowns w(t, x) are approximated by their mean values on the cell C k at time t:

w k (t) 1 h k C k w (t, x) dx .
Integrating (7.1) over each cell and applying the Green's formula we obtain:

∂w k ∂t + F k+1/2 (t) -F k-1/2 (t) h k = 0 , (7.2)
where F k±1/2 stands for the numerical flux at the interface x k±1/2 (see for instance [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Godlewski | Numerical approximation of hyperbolic systems of conservation laws[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]).

In our case, F = f (R(0 ± , w L , w R )) is determined from the exact solution R(0 ± , w L , w R ) of the Riemann problem associated with the left w L and the right w R state:

w 0 (x) = w L if x < 0 , w R if x > 0 .
Equation (7.1) is completed with an entropy inequality of the form:

S(w) = ∂s(w) ∂t + ∂ψ(w) ∂x ≤ 0 , (7.3)
where (s, ψ), with ψ (w) = s (w) f (w), stands for a convex entropyentropy flux pair. S is called the density of entropy production. 3) leads to the so-called numerical density of entropy production which is a measure of the amount of violation of the entropy equation (as a measure of the local residual as in [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF][START_REF] Houston | A posteriori error analysis for numerical approximations of Friedrichs systems[END_REF][START_REF] Karni | Local error analysis for approximate solutions of hyperbolic conservation laws[END_REF][START_REF] Karni | A smoothness indicator for adaptive algorithms for hyperbolic systems[END_REF]). As a consequence, the numerical density of entropy production provides information on the need to locally refine the mesh (e. g. if the solution develops discontinuities) or to coarsen the mesh (e. g. if the solution is smooth and well-approximated) and can be useful for AMR methods (as already done in [START_REF] Puppo | Numerical entropy production on shocks and smooth transitions[END_REF][START_REF] Puppo | Numerical entropy production for central schemes[END_REF][START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF] and [1, 59, 60, 72-74, 121, 122, 166, 167]).

The Equation (7.3) is discretized using the same scheme as Equation (7.1). Thus, the semi discrete formulation for the entropy inequality reads:

S k (t) = ∂s(w k ) ∂t + δψ k (t) h k (7.4)
where

δψ k (t) = ψ R(0 -, w k (t), w k+1 (t)) -ψ R(0 + , w k-1 (t), w k (t)) .
To also achieve a second order approximation in space, the Monotone Upwind Scheme for Conservation Laws (MUSCL) reconstruction can be used as done in [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF].

The numerical time integration of Equations (7.2) and (7.4) can be performed through Runge-Kutta or Adams-Bashforth (we refer to [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF] for details). Thus, for instance, the first order space-time scheme is :

w n+1 k -w n k δt n + F n k+1/2 -F n k-1/2 h k = 0 ,
and

S n k := s(w n+1 k )) -s(w n k ) δt n + δψ k n h k (7.5)
where δt n is the time step,

F n k+1/2 = f (R(0 ± , w n k , w n k+1 ))
and,

δψ n k = ψ R(0 -, w n k , w n k+1 ) -ψ R(0 + , w n k-1 , w n k ) .

adaptive mesh refinement algorithm

We recall here the adaptive mesh refinement algorithm using the numerical density of entropy production, for instance Equation (7.5), as a mesh refinement indicator.

Mesh refinement criterion based on the numerical density of entropy production

The construction of an efficient criterion remains a difficult subject, since it depends on physical phenomena and numerical methods.

Multi-resolution analysis (see for instance, [START_REF] Harten | Multiresolution algorithms for the numerical solution of hyperbolic conservation laws[END_REF][START_REF] Müller | Adaptive multiscale schemes for conservation laws[END_REF]) provides some robust tools to estimate and control rigorously the error at each refinement level. The a posteriori error estimates can be also based on the measure of the local residual as done in [START_REF] Verfürth | A posteriori error estimation and adaptive mesh-refinement techniques[END_REF] and in [START_REF] Karni | A smoothness indicator for adaptive algorithms for hyperbolic systems[END_REF][START_REF] Zhang | A posteriori error estimation for finite-volume solutions of hyperbolic conservation laws[END_REF] for instance. Based upon Richardson-type estimates of the local truncation error [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF], one can also refine or coarsen to reach a given accuracy for a minimum amount of work. However, near discontinuous solutions, Richardson extrapolation is invalid [START_REF] Powell | Adaptive-mesh algorithms for computational fluid dynamics[END_REF]. It will work but it is no more reliable than heuristic methods and the extra cost cannot be justified. The mesh refinement criterion can be also constructed, as simple as possible, as the combination of gradients based on the underlying physical quantities to detect pertinent regions to refine or to coarsen [START_REF] Quirk | A parallel adaptive grid algorithm for computational shock hydrodynamics[END_REF] under the condition that the mesh refinement threshold is wellchosen.

But paradoxically, to our knowledge, very few works use a physical criterion as [START_REF]Contribution à l'étude théorique et à l'approximation par éléments finis du système hyperbolique de la dynamique des gaz multidimensionnelle et multi-espèce[END_REF]: the numerical density of entropy production. More recently, Puppo [START_REF] Puppo | Numerical entropy production on shocks and smooth transitions[END_REF][START_REF] Puppo | Numerical entropy production for central schemes[END_REF] has proved that the numerical density of entropy production can be used as a discontinuity indicator and a local error indicator. Indeed, the entropy inequality being related mathematically and physically to the system (7.1), the numerical density of entropy production (7.5) should be (and actually is, as we will see in Section 8.3, see also [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF][START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF] for instance) an accurate and a useful mesh refinement criterion whenever the mesh refinement threshold is well-calibrated (we refer to Chapter 8 for a possible construction of an automatic "well-calibrated" mesh refinement threshold).

Mesh refinement process

We present here a local mesh refinement procedure driven by the numerical density of entropy production, for instance Equation (7.5). In order to reduce the time necessary to manage the refinement, we use "macro cells" which could be refined by generating hierarchical grids. Each cell can be split in two sub-cells. We thus produce a dyadic cells graph, whose numbering (in basis 2) allows a quick computing scan to determine the adjacent cells. For the convenience of the reader, we make use of the following notations: let k b be the index which makes reference to the macro cell numbered k and b a binary number which contains the hierarchical information of a sub-cell. In particular, the level of a sub-cell C k b is defined as the length(b) -1. For instance, a macro cell C k 0 of level 0 will be split into two sub-cells C k 00 and C k 01 of level 1. A mesh refinement example is proposed in Figure 7.1. The mesh refinement procedure can be simply expressed as follows. We introduce a non negative mesh refinement criterion S n k b , for instance Equation (7.5), computed on each cell C k b at time t n and compared to a threshold α, for instance,

α = S m = 1 |Ω| ∑ k b S n k b
In this case, this approach is the so-called mean method (see [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] for other methods).

We then define two coefficients 0 < β min ≤ β max , which determine the ratio of the cells to be refined or coarsened. Thus, for each cell C k b :

AL 1. • if S n k b > α max = β max α, the cell is refined and split into two sub-cells C k b0 and C k b1 . • if S n k b0 < α min = β min α and S n k b1 < α min , the cell is coarsened into a cell C k b .
The adaptive algorithm is stopped when the level l of a cell C k b reaches the maximum level l max .

Remark 7.2.1.

• The threshold parameters β min and β max allow to set a percentage of mesh refinement and mesh coarsening with respect to the quantity S m . The grid is not refined or coarsened for S n k b ∈ [α min , α max ].

• It is not surprising that these settings can deteriorate or improve the accuracy of the numerical solution. For instance, the more β min and β max are small, the more accurate are the results at the expense of the computational time. Thus β min and β max are a problem-dependent parameters whose choice might be the critical weakness of the mesh adaptation method. This issue is addressed and an automatic selection of the threshold is proposed in Chapter 8.

On one hand, for the first order scheme, if a cell C k b is split into two sub-cells C k b0 and C k b1 , averaged values at time t n are projected on each sub-cell:

w n k b0 = w n k b1 = w n k b
and the numerical fluxes between the two sub-cells (see Figure 7.2a) are defined as:

                 F n k b0 -1/2 = F n k b -1/2 , F n k b0 +1/2 = F n k b1 -1/2 = f (w n k b ) , F n k b1 +1/2 = F n k b +1/2 .
In the case of second order method, as the gradient is available, we also use:

         w n k b0 = w n k b - h k b 4 ∂w n k b ∂x , w n k b1 = w n k b + h k b 4 ∂w n k b
∂x . On the other hand, if two sub-cells C k b0 and C k b1 are coarsened, we initialize the new cell C k b as (see Figure 7.2b):

w n k b = 1 2 w n k b0 + w n k b1
and

       F n k b -1/2 = F n k b0 -1/2 , F n k b +1/2 = F n k b1 +1/2 .
(a) Mesh refinement.

(b) Mesh coarsening.

the local time stepping method

Unlike the Runge-Kutta method, the Adams-Bashforth time integration can be easily implemented in the local time stepping framework in order to reduce the computational cost. The main approach to update averaged quantities will be recalled and the general algorithm will be presented.

The local time stepping algorithm

Explicit adaptive schemes are well-known to be time consuming due to a Courant, Friedrichs and Lewy (CFL) stability condition. The CPUtime increases rapidly as the mesh is refined since the CFL stability condition imposes an upper bound on δt h where δt is the time step and h the finest mesh size. Nevertheless, the CPU-time can be significantly reduced using the local time stepping algorithm (see e. g. [START_REF] Altmann | A local time-stepping discontinuous Galerkin algorithm for the MHD system[END_REF][START_REF] Müller | Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping[END_REF][START_REF] Tan | Moving mesh methods with locally varying time steps[END_REF]). In recent years, such numerical schemes have been widely developed for equations which arise in many fluid flows: traffic flows, multiphase flows, multi-layer flows, etc. The local time stepping algorithm has been employed successfully for real two and three dimensional steady and unsteady problems (see for instance [START_REF] Jenny | Adaptive multiscale finite-volume method for multiphase flow and transport in porous media[END_REF][START_REF] Müller | Fully adaptive multiscale schemes for conservation laws employing locally varying time stepping[END_REF][START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF][START_REF] Puppo | Numerical entropy and adaptivity for finite volume schemes[END_REF][START_REF]Dynamic adaptivity and residual control in unsteady compressible flow computation[END_REF][START_REF] Tang | A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space grids[END_REF]).

In this section, we combine the local time stepping algorithm (see for instance [START_REF] Altmann | A local time-stepping discontinuous Galerkin algorithm for the MHD system[END_REF]) with the mesh refinement process presented above. For the sake of clarity, we do not use the binary subscript notation introduced in the previous section.

Let h min be the minimum diameter of the mesh and h max the maximum diameter of the mesh. All cells of the mesh are sorted in groups, called "level", based on their diameter. Thus, the level of refinement L k of the k th cell C k is defined by:

2 N-L k h min h k < 2 N+1-L k h min and N = log 2 h max h min + 1
where N stands for the maximum level. The coarsest cells are therefore of level 1, while the finest cells are of level N.

Remark 7.3.1. For any mesh, the previous inequality holds and in the particular case of a dyadic mesh, we have

h k = 2 N-L k h min .
Let δt n be the minimum time step at time t n according to the CFL condition associated with the smallest cell. The macro time step ∆t n is then defined by ∆t n = 2 N-1 δt n .

We define the level of interface of two adjacent cells L k+1/2 by

L k+1/2 = max (L k , L k+1 ) .
Assuming that the maximum level of refinement is N at the current time t n , following [START_REF] Altmann | A local time-stepping discontinuous Galerkin algorithm for the MHD system[END_REF], the local time stepping algorithm reads as follows:

foreach i ∈ {1, end end Remark 7.3.2. This scheme uses the Osher and Sanders [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] projection. This point is motivated in Section 7.3.2.

Let us illustrate for a global time advancement the local time stepping algorithm for a first order scheme (the principle for a second order is essentially the same since we have to use the value of the flux stored in the two previous local time steps). We consider a mesh composed of 4 cells of levels 1, 2 and 3 respectively, i.e. here N = 3, as displayed in Figure 7 Let us note w n m the state vector at time t n . At the first stage (see Figure 7.3), i = 1, only the small cells are advanced with a step δt n , i.e. w n 1 3 := w 3 (t n + δt n ) and w n 1 4 := w 4 (t n + δt n ) are computed. While i = 2 (see Figure 7.4) concerns the small cells of level L N and the cells immediately larger at level L N-1 . The cells L N travel with a time step δt n , while the cells L N-1 travel with a time step 2δt n , etc.

Computation of the flux

We recall the procedure introduced by Osher-Sanders (OS) [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] and, Tang-Warnecke (TW) [START_REF] Tang | A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space grids[END_REF] for the computation of the flux at the cell interface of two different levels. For the sake of simplicity, we assume here only two time increments which means that we consider cells of only two different kinds, say h(x) = h k 0 = h 0 for x < x k 0 +1/2 and h(x) = h k 0 +1 = h 0 /2 for x > x k 0 +1/2 . Thus, for every k k 0 , all cells C k are of level L k 0 and for every k > k 0 , all cells C k are of level L k 0 + 1 (as displayed in Figure 7.7). Let ∆t n := t n+1t n = 2δt n be the macro time step. To avoid boundary terms, we consider the example with an initial condition with compact support.

A simple projection method

The OS scheme is defined as follows (see Figure 7.7a): ∀k ∈ Z,

w n+1/2 k = w n k if k k 0 , w n k -∆t n 2h k δF n k if k > k 0 w n+1 k = w n k -∆t n 2h k δF n k + δF n+1/2 k . ( 7.6)
It was pointed out by Tang and Warnecke [START_REF] Tang | A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space grids[END_REF] that the present scheme loses the conservative property and it is locally non consistent in the sense of the truncation error. Indeed, at the intermediate time step, one has The fact that the scheme is not consistent in the sense of truncation error can be observed in the simplest case of the one-dimensional transport equation (as done in [START_REF] Tang | A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space grids[END_REF]):

∑ k∈Z w n+1/2 k = ∑ k∈Z w n k + ∆t n 2h k 0 +1 F n k 0 +1/2 (w n k 0 , w n k 0 +1 )
∂ t w + c∂ x w = 0
with a positive constant c. In that case, using the Godunov solver, the scheme (7.6) writes for k = k 0 :

w n+1 k 0 +1 -w n k 0 +1 ∆t n - c 2h k 0 w n k 0 +1 -w n k 0 + w n+1/2 k 0 +1 -w n k 0 = 0 .
Assuming w smooth enough, writing the Taylor expansion of w at point x k 0 +1 at time t n+1 and t n+1/2 , we deduce the following equation

∂ t w + c∂ x w = - c 4 ∆t h ∂ t w + O(h, ∆t).
Thus, if ∆t h converge to a constant d > 0 then the previous equation does not converge to the initial transport equation.

A natural projection method

Contrary to the Osher and Sanders [START_REF] Osher | Numerical approximations to nonlinear conservation laws with locally varying time and space grids[END_REF] approach, Tang and Warnecke [START_REF] Tang | A class of high resolution difference schemes for nonlinear Hamilton-Jacobi equations with varying time and space grids[END_REF] propose the following scheme: ∀k ∈ Z,

w n+1/2 k = w n k -∆t n 2h k δF n k w n+1 k = w n+1/2 k -∆t n 2h k δF n k if k k 0 δF n+1/2 k if k > k 0 .
One can easily check that the global conservativity is also lost and

∑ k∈Z w n+1 k = ∑ k∈Z w n k + ∆t n 2h k 0 +1 F n+1/2 k 0 +1/2 (w n k 0 , w n k 0 +1 ) .
Although the scheme is not conservative, they claim that it is consistent in the sense of the truncation error.

Even if the OS projection leads to non conservative and non consistent scheme at the space-time grid between two levels of refinement, it has the advantage to be less time consuming than the second one. In general, the apparent loss of consistency on the local truncation error seems not to affect the error of the scheme (see [START_REF] Cockburn | A priori error estimates for numerical methods for scalar conservation laws. Part II: Flux-splitting monotone schemes on irregular Cartesian grids[END_REF][START_REF] Simeoni | Remarks on the Consistency of Upwind Source at Interface Schemes on Nonuniform Grids[END_REF] for interesting issues).

Moreover, one can "reinforce" the consistency and the stability of the OS projection by preventing two adjacent cells from having a level difference greater than two. This trick is referred as "smoothing grid" in the following.

Comparing numerically the OS and TW projection methods, in the case of the Sod's shock tube problem for the one-dimensional gas dynamics equations for ideal gas, we get more stable results with the smoothing grid technique and we save the computational time (see [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF] for more details) with respect to the TW projection: CPU-time for TW projection = 288.21 (s), CPU-time for OS projection = 80.04 (s).

data management in high dimension

The multi-dimensional method is almost the same as described in the previous Section 7.1, Section 7.2 and Section 7.3 to the exception of the data management.

Keeping in mind the following requirements: parallel treatment, fast computation and hierarchical grid (mesh level cannot exceed 2 for stability reasons as pointed above in Section 7.3.2 and in [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF]), the main and challenging difficulty is to develop a suitable mesh refinement tool which is related on the design of software data structures and numerical scheme.

As pointed out by several authors, the extension toward threedimensional case is not an easy task. Interesting works have been presented for 2d Cartesian grid or quad-tree [START_REF] Berger | Local adaptive mesh refinement for shock hydrodynamics[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive cartesian grids[END_REF][START_REF] Yiu | Quadtree grid generation: Information handling, boundary fitting and CFD applications[END_REF][START_REF] Zhang | A two dimensional hydrodynamic and sediment transport model for dam break based on finite volume method with quadtree grid[END_REF], octree for 3d simulations [START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF][START_REF] Losasso | Simulating Water and Smoke with an Octree Data Structure[END_REF], and anisotropic AMR [START_REF] Coupez | Solution of high-Reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF][START_REF] Hachem | Immersed stress method for fluid structure interaction using anisotropic mesh adaptation[END_REF].

We aim to generalise the scheme proposed Section 7.1, Section 7.2 and Section 7.3 (see [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF]) with a finite volume solver. This naturally leads to quad-tree meshing in two dimensional problem and octree meshing in three dimensional problem.

In the most of the applications, the presence of a complex moving air-water interface or more generally wave interactions implies to remesh at each time step, which is clearly time consuming. Keeping in mind our willingness to find a relevant compromise between the contradictory aims of physical relevancy and numerical performance, we introduce a Cartesian block-based mesh approach (somehow like in [START_REF] Williamschen | Parallel Anisotropic Block-Based Adaptive Mesh Refinement Algorithm For Three-Dimensional Flows[END_REF][START_REF] Zheng | Block-Based Adaptive Mesh Refinement Finite-Volume Scheme for Hybrid Multi-Block Meshes[END_REF]). The computational domain is split in several blocks. In practice, each block corresponds to the initial unstructured mesh composed of hexahedral cells. According to the mesh refinement procedure (with respect to S n k , S m , the mesh refinement level N), each block is sub-divided in a Cartesian way (2 N-1 n x , 2 N-1 n y , 2 N-1 n z ) where (n x , n y , n z ) stands for the initial block discretisation.

As done in Section 7.1, Section 7.2 and Section 7.3 (see [START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF]) , for each refined cell (or block), averaged values at time t n are projected on each sub-cell and each flux is computed as simply as possible to avoid heavy computation.

The interface defining two blocks is most of the time is a non conforming one as displayed in Figure 7.9b.

Next, in order to spare proportionally the balanced distribution of CPU load, the cells of each block are re-distributed in a fixed number of domains according to the Cuthill-McKee numbering. The number of domain being fixed, each domain are loaded in a given MPI process. The procedure is certainly not optimal but the synchronisation time between each domain are almost the same, hence overall more efficient. Finally, the re-numbering and re-meshing being expensive, the mesh is kept constant on a time interval, called AMR time-step, given by the smallest block (and not by the smallest cell) and the maximum velocity. This procedure provides a good compromise between computational time and accuracy (see for instance [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF][START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -Application to wave breaking[END_REF] and in Section 8.3).

Let us illustrate the algorithm to define unstructured meshes using the BB-AMR scheme in the two-dimensional case (i.e. n z = 1 always). During the refinement process, the number of blocks and domains are kept always constant. As displayed in Figure 7.8a, we consider an initial mesh composed of 9 blocks, assumed to be of level 1 (the coarsest level with n x = n y = 1), and decomposed on 3 domains according to a Cuthill-McKee numbering of the blocks (3 + 3 + 3 internal cells). For example the first domain (blue) is composed of three internal cells and three adjacent "ghost" cells in order to share to send/receive information between domains. After some computations, at the AMR-time T 1 , let us suppose that blocks 1 and 4 have to be refined. As displayed in Figure 7.8b, the blocks distribution into the domains has been modified (5 + 5 + 5 internal cells). One can observed that the new mesh is a non-conforming mesh. For example, the cell belonging to the block 2 is bordered by 7 faces. If, at the AMR-time T 2 , only the block 1 exceeds the mesh refinement criterion, it has to be refined. But, as the level of two adjacent blocks cannot exceeds 2 in order to avoid oscillations, the block 2 must also be refined. As displayed in Figure 7.8c, the blocks distribution has been modified (16 + 7 + 7 internal cells). The distribution is not well balanced because we distribute the blocks and not the cells for sake of simplicity. This drawback is reduced as the number of blocks is much bigger than the number of domains.

The three dimensional case is dealt in the same way. For illustration purpose, Figure 7.9 shows an example of 3d block-based mesh decomposition (64 + 39 + 39 internal cells) with 3 domains and 27 blocks, where the first block in the corner is of mesh refinement level 3.

A U T O M AT I C T H R E S H O L D I N G M E T H O D

In Chapter 7, we have overlooked the mesh refinement threshold which is a problem-dependent parameter whose choice might be the critical weakness of mesh adaptation methods. In this chapter, we address a new construction of an automatic selection of a suitable mesh refinement threshold. It allows to detect almost all relevant regions to refine without any hand-calibration as usually done for mesh adaptation algorithms (see AL 1 and [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF]).

The construction of the automatic selection of the threshold is based on the decreasing rearrangement function of the mesh refinement criterion. If S denotes the mesh refinement criterion, the decreasing rearrangement is the Lebesgue measure of the set {S(x) > α}. This function sorts the local maxima of S from the smallest to the largest. One can therefore set automatically the threshold small enough to detect the most relevant cells to refine. In practice, we are able to reach a good balance between the accuracy and the computational time without any hand-calibration as usually done in mesh adaptation algorithms as shown in Section 8.3 (see also [START_REF] Pons | Adaptive mesh refinement method. Part 2: Application to tsunamis propagation[END_REF]).

In what follows, we define the piecewise mesh refinement criterion as follows :

S(x, t) = ∑ k S n k 1 C k ×[t n ,t n+1 [ (x) (8.1)
where

1 A (x) = 1 if x ∈ A , 0 otherwise.

settings and definition

Let us consider a smooth (at least twice differentiable) mesh refinement criterion S(x, t) ∈ R + , x ∈ [0, L] and t > 0 where L is the length of the domain. The time t being fixed, we write in the sequel S = S(x).

Without loss of generality, we suppose that S(0) = S(L) = S (0) = S (L) = 0 and 0 < S ∞ = max x∈(0,L) S(x) < ∞ (if S ∞ = 0 the numerical solution is identically equal to zero). The function S can be regarded here as the regularisation of the discrete mesh refinement criterion (8.1).

In view of the above assumptions, the set

Z α = {x ∈ (0, L); ϕ α (x) = S(x) -α = 0 and S (x) = 0}
is not empty. Indeed, for each 0 < α < S ∞ , since S ∞ > 0, S has at least one maximum. Then, there exists p α ∈ N * such that the number of elements in the set

Z α is #Z α = 2p α .
Thus, for all α ∈ (0, S ∞ ), one can describe the set Z α as follows

Z α = {x 0 (α) < x 1 (α) < • • • < x 2p α -2 (α) < x 2p α -1 (α)} .
Let us assume that S has p local maxima. Then there exists an increasing sequence α * k 1≤k≤p and a sequence

x * k 1≤k≤p such that ∀k = 1, . . . , p S (x * k ) = 0, S (x * k ) < 0 and S(x * k ) = α * k .
By construction, the sequence α * k 1≤k≤p represents the local maxima sorted from the smallest to the largest with α * p = S ∞ . Let us also remark that x * k / ∈ Z α , ∀k = 1, . . . , p, ∀α ∈ (0, S ∞ ). With these settings, one can define the decreasing rearrangement function of S:

α ∈ [0, S ∞ ] → d(α) :=      L if α = 0 , ∑ p α k=1 x 2k+1 (α) -x 2k (α) if 0 < α < S ∞ , 0 if α = S ∞ . (8.2)
We have

Theorem 8.1.1. Let l ∈ N * , S ∈ C l ([0, L], R + ) be a smooth function such that S(0) = S(L) = S (0) = S (L) = 0, 0 < S ∞ < ∞ with p local maxima. Then 1. S ∞ 0 d(α) dα = L 0 S(x) dx. 2. d ∈ C 0 ([0, S ∞ ], R + ) and a) ∀α ∈ [0, S ∞ ], d (α) < 0 b) ∀k ∈ 0, p , lim α→α * k d (α) = -∞ with the convention α * 0 := 0 3. d ∈ C l (D * , R + ) on the set D * := p-1 k=0 (α * k , α * k+1
) . Let us now illustrate the previous introduced notations using the given mesh refinement criterion displayed in Figure 8.1a. The decreasing rearrangement function d of S is represented in Figure 8.1b. The first and the second order derivatives of d are shown in Figure 8.1c and Figure 8.1d. In these figures, the distribution is normalised by L and we have p = 2, i.e. p local maxima.

In Figure 8.1a, for α = 0.5, the set Z α is composed of 4 = 2p α elements which are approximately

x 0 (α) ≈ 0.8803, x 1 (α) ≈ 1.6196, x 2 (α) ≈ 3.379, x 3 (α) ≈ 4.121 .
The distribution function for α = 0.5 is given by

d(α) = (x 1 (α) -x 0 (α)) + (x 3 (α) -x 2 (α)) L ≈ 0.3 .
In Figure 8.1b (see also Figure 8.1a), each vertical asymptote of d indicates a local maximum of S(x). More precisely, we have

α * 1 = S(x * 1 = 1.25) = 1 and α * 2 = S(x * 2 = 3.75) = 2 = S ∞ where (x * k ) 1≤k≤2
are the local maxima of S. Thus, by construction, (α * k ) 1≤k≤2 is an increasing sequence.

construction of the threshold

As pointed out before in [START_REF] Dannenhoffer | Grid adaptation for complex two-dimensional transonic flows[END_REF][START_REF] Powell | An embedded mesh procedure for leading-edge vortex flows[END_REF], it is difficult to define an automatic threshold from the decreasing rearrangement function (8.2) without using the first order or/and the second order derivatives of d to characterise the local maxima. Moreover, we emphasise in [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] that the use high order derivates requires an accurate interpolation of d to be accurate in the localisation of this local maxima which can be quite costly. To this purpose, in order to avoid these disadvantages, we have proposed the following weight function (8.3). This function is built only from the decreasing rearrangement function (8.2):

f (α) = αd(α) (8.
∀k = 0 . . . p -1, α * * k+1 ∈ (α * k , α * k+1 ) such that d (α * * k+1 ) = 0 .
Moreover, the function f (8.3) has p local maxima α 1 , . . . , α p such that

∀k = 1 . . . p, α k ∈ (α * * k , α * k ) .
As emphasised in [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF], α = S m := 1 |Ω| Ω S(x) dx, i.e. β min = β max = 1, is in general relevant to detect cells to refine for smooth solutions and α < S m , i.e. β < β min = β max = 1, for discontinuous solutions (see AL 1).

Thus, by virtue of the Corollary 8.2.1 and the previous statement, it is natural to look for an automatic threshold as follows:

α PE such that f (α PE ) = max 0<α≤S m f (α) .
In order to illustrate the previous statements, we consider the function S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25 and x B = 3.75. We represent in Figure 8.2 the shape of the function f (8.3) and in Figure 8.3 the mesh refinement criterion for a discontinuous and a smooth solution. As expected, we get • for a discontinuous solution, see Figure 8.2a and 8.3a, α PE < S m .

• for a smooth solution, see Figure 8.2b and 8.3b, α PE = S m .

In both cases, the computed threshold is pertinent to detect the region to refine. Finally, the discretisation of the decreasing rearrangement function (8.2) is performed as follows

d(α) = 1 N M-1 ∑ j=0 d j 1 (α j ,α j+1 ) (α) with d j = #{k ; S n k > α j }
where # is the number of elements in the set {S n k > α j } and where α j is defined by

α j = S m j M , 0 ≤ j ≤ M. (8.4)
In [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF], we show that for M ≥ N, the threshsold α PE is quasiinsensitive and make our method totally parameterless and less expensive that the methods proposed in [START_REF] Dannenhoffer | Grid adaptation for complex two-dimensional transonic flows[END_REF] and [START_REF] Powell | An embedded mesh procedure for leading-edge vortex flows[END_REF].

The algorithm AL 1 is then defined by

α min = β min α PE and α max = β max α PE
where 0 < β min ≤ β max .

numerical results

To illustrate the ability of the method to catch relevant information, we start first with simple waves in the case of a one dimensional Riemann problem in Section 8.3.1. Then, the case of wave interactions is dealt in Section 8.3.2 for the Monai-Walley tsunami run-up onto a complex three dimensional beach. We also show the efficiency of the BB-AMR presented in Section 7.4. For both test case here, we have considered the Saint-Venant equations :

∂ t h + div(hu) = 0 ∂ t (hu) + div hu ⊗ u + g 2 h 2 I = -gh∇Z (8.5)
where the unknowns h(x, t) and u(x, t) = (u 1 (x, t), u 2 (x, t)) are respectively the height of the water and the depth-averaged velocity of the water at a space-time point (x, t), x = (x 1 , x 2 ) ∈ R 2 is the space coordinate, t > 0 is the time, g is the gravitational constant set to 9.81m/s 2 , Z(x) is the bottom surface elevation and I is the identity matrix. I also refer to [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF][START_REF] Ersoy | Adaptive multi scale scheme based on numerical density of entropy production for conservation laws[END_REF][START_REF] Ersoy | Adaptive scheme based on entropy production: robustness through severe test cases for hyperbolic conservation laws[END_REF][START_REF] Golay | Entropy production as mesh refinement criterion -Application to wave breaking[END_REF][START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 2: Application to tsunamis propagation[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -application to wave breaking[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -Application to wave breaking[END_REF] for more test cases.

1d Riemann problem over a flat bottom

We consider the Riemann problem with the initial data ∀x ∈ [0, 80], (h(x, 0), u(x, 0))=

(5.64,

8) if x < 20 , (0.6, 8) if x > 20 (8.6)
where the exact solution is composed of a left-going rarefaction and a right-going shock wave (see for instance Figure 8.4a, black line).

The goal is to show that, whatever the mesh refinement criterion is, the threshold α PE is selected automatically to adapt the mesh without hand-calibration. To this purpose, we propose : criterion 0: Since the exact solution of the Riemann problem (8.5)-(8.6) is known, we consider the following criterion

S n k := |h n k -h ex (x k , t n )|
where h ex stands for the exact solution of (8.5)-(8.6).

criterion 1: We consider the numerical density of entropy production:

S n k = |E n k | where E n k = E(w n+1 k ) -E(w n k ) δt n + G k+1/2 (w n k , w n k+1 ) -G k-1/2 (w n k-1 , w n k ) δx k with E(w) = q 2 2h + gh 2 2 , q = hu and G(w) = E(w) + gh 2 2
u . The numerical flux G k+1/2 is computed with the exact local solution of the 1D Riemann problem at the interface x k+1/2 with the left state w n k and the right state w n k+1 as done in Chapter 7.

criterion 2: Finally, the last mesh refinement criterion is based on the gradient method. It can be constructed as a combination of gradients based on physical quantities to detect pertinent regions to refine or to coarsen [START_REF] Quirk | A parallel adaptive grid algorithm for computational shock hydrodynamics[END_REF]. To make it as simple as possible, we consider the gradient of h which can be an efficient mesh refinement criterion if the threshold is well-chosen. Thus, for instance, we define

S n k = h n k+1 -h n k x k+1 -x k .
Each numerical experiment starts initially with N(0) = 100 cells, the maximum level of refinement is set to l max = 3, M = N(t) for all t ≥ 0 (see Eq. (8.4)) and we prescribe free boundary conditions.

Even if we perform the following strategy for each cell C k :

• if S n k > α PE , the mesh is refined and split,

• if S n k < α PE the mesh is coarsened which is known to be a source of spurious oscillations, the proposed method does not generate spurious oscillations in the computation of the numerical solution (see [START_REF] Pons | Adaptive mesh refinement method. Part 1: Automatic thresholding based on a distribution function[END_REF] for further details) at least for simple waves as shown in Figure 8.4. Moreover, whatever the mesh refinement criterion, we obtain almost the same numerical solution with almost the same order of accuracy as presented in Table 8.1 (see also Figure 8.4). Finally, in comparison with the mean method, if the threshold is not well-calibrated, the region of interests are refined or too refined as displayed in Figure 8.5.

N m (l) → ||h -h ex || 1 N m (l) → ||u -u

Tsunami runup onto a complex three dimensional Monai-Walley beach

We now focus on the overall performance of the BB-AMR (presented in Section 7.4) scheme with automatic thresholding compared to experimental and state-of-the-art numerical results 1 . For each test case, the numerical simulations are carried out on adaptive grids and on uniform grids. The minimum adaptive cell size δx adaptive is equal to the uniform one δx uni f orm . The comparisons of both simulations allow to show the reliability and the efficiency of the BB-AMR method. To simplify the analysis, we have considered only 1 domain, i.e., 1 processor. Interested readers can found a detailed numerical study with several domains in [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF]. For all computations,

• the MUSCL reconstruction and the second order Runge-Kutta scheme (RK2) method are used to achieve a second order scheme in space and time.

• to simplify the analysis and to focus on the automatic threshold, we have chosen β min = β max = 1.

• M = N (see Equation (8.4)).

• the CFL number is set to 0.9,

• For the sake of simplicity, the re-meshing time step is fixed constant to 0.25 s. This choice roughly verifies the macro CFL condition.

This test case concerns the Hokkaido-Nansei-Oki tsunami, in 1993, that struck the Okushiri Island in Japan. The tsunami run-up records was about 30 m height and the currents speed of order 10-18 m/s for which relevant high-quality data were saved. 2 The largest value run-up (32 m) was recorded near the Monai-Walley beach.

Experimental settings

To understand this complex run-up, a 1/400 scale laboratory model of Monai was realized in a large-scale tank of 205 m long, 6 m deep, 1 The numerical soft used for the following test cases is the EOLENS code developed by the Institut de Mathématiques de Toulon (IMATH) and Principia. The code is written in Fortran 90, compiled with gfortran on Linux machine Intel(R) Core(TM) i5-2500 CPU @ 3.30GHz.

2 Several sources and data can be found, see for instance http://nctr. pmel.noaa.gov/benchmark/Laboratory/Laboratory _ MonaiValley/ or http://isec. nacse.org/workshop/2004 _ cornell/bmark2.html and 3.5 m wide at the Central Research Institute for Electric Power Industry (CRIEPI) in Abiko (Japan). The offshore incident wave on a water depth d = 13.5 cm is prescribed. There are reflective side-walls at y = 0 and 3.5 m as for the x boundaries. The bathymetry as well as the coastal topography reproduced in the laboratory experiment are represented 3 in Figure 8.6b and Figure 8.6c .

The input wave at x = 0 m is a leading-depression height of -2.5 mm with a crest of 1.6 cm, as displayed in Figure 8.6a. In the experiment, the waves are measured at thirteen locations and complete time histories are given at three locations, (x, y) = (4.521, 1.196), (4.521, 1.696), and (4.521, 2.196) in meters, see also [START_REF] Liu | Advanced numerical models for simulating tsunami waves and runup[END_REF]. In contrast with the previous experimental test cases, the dispersive phenomena can be fully neglected here.

Numerical vs experimental results.

The initial condition for this test case is a lake at rest and an imposed water height on the left side wall, see Figure 8.6a. Wave gauges data 4 experimental results are given up to time t = 100 s. However, the information concerning the input wave profile (see Figure 8.6a) at x = 0 m is only available up to time t = 22 s. For numerical purpose, we prescribe for t > 22 s, a constant free surface level η = 0 m on the left wall. In practice, due to this assumption, some artificial reflecting waves can occur on the wave-maker beyond t > 22 s and can reach the free surface gauges around t = 30 s. Finally, to avoid any wrong comparisons due to these artificial reflections, the comparisons between the experimental gauges and the numerical results will be only done over the first 30 s of the experiment 5 For computational purpose, we have considered 240 initial blocks composed of 8 000 (varying up to 40 000) cells for the initialisation of the adaptive simulation and 62 000 cells for uniform mesh simulation. The time simulation is 30s. Reflecting boundary conditions are prescribed to walls.

To illustrate this test case we show in Wave gauges data are available at the address http://nctr.pmel.noaa.gov/ benchmark/Laboratory/Laboratory _ MonaiValley/ However, let us note that in [START_REF] Dutykh | The VOLNA code for the numerical modeling of tsunami waves: Generation, propagation and inundation[END_REF][START_REF] Dutykh | Mathematical modelling of tsunami waves[END_REF], the simulation is done over 30 s and can be a point of comparison for future work. At the time of writing [START_REF] Pons | Adaptive mesh refinement method. Part 2: Application to tsunamis propagation[END_REF], we are unaware of Dutykh' results [START_REF] Dutykh | The VOLNA code for the numerical modeling of tsunami waves: Generation, propagation and inundation[END_REF][START_REF] Dutykh | Mathematical modelling of tsunami waves[END_REF]. 8.11c). These figures show that the wave amplitude and the wave shape are accurately predicted at three different location in the computed domain. As emphasised before, dispersive effects being negligible, the numerical results are accurate. However, it can be noticed in Figure 8.11 a small amplitude errors occurring at the beginning of the simulations. This error can be attributed to the lack of accuracy of the prescribed initial condition and boundary condition at x = 0. This problem was already encountered by Popinet [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF].

AMR performance

To highlight the efficiency of the mesh refinement parameter and the automatic selection of the threshold, the numerical density of entropy production (without absolute value) and the associated mesh are displayed at several times in Figure 8.12-Figure 8.15, (time t = 11.25 s in Figure 8.12, for t = 13.25 s in Figure 8.13, for t = 16 s in Figure 8.14 and for t = 17.5 s in Figure 8.15). We observe mainly that the wave front and subsequent waves are very well-captured by the criterion and the mesh refinement threshold is well chosen to refine in those areas (see for instance Figure 8.10 and Figure 8.15). Between 3 and 7 seconds, the tsunami propagates in the domain involving a grow-up of the number of cells. Then the wave passes through the island and reaches the reef around 15 seconds. After this, the wave is reflected and goes back to the wave maker with a second interaction with the island. The maximum number of cells is reached when the wave is passing on the island and begins to run-up the coast (see Figure 8.13), and when the wave goes back to the wave maker. These results are coherent with the physical process involved since it corresponds to the both times where multiple active flow regions are observed. The overall scales of the tsunami propagation with its reflections are very well captured as displayed in Figure 8.12, Compared to the computation on the uniform grid (see Figure 8.11), the adaptive mesh refinement method requires an average of almost 25000 cells against 62000 cells. For the same accuracy, the AMR method allows to save almost 60% of cells with respect to the uniform sim- ulation. Thus, the AMR method allows to speed up the computation almost 3 times.

In this work, a multidimensional finite volume solver is developed to solve non-linear hyperbolic equation system on unstructured meshes.

An AMR is introduced to improve both solution accuracy and CPU performance. It is based on a useful numerical criterion: the numerical density of entropy production. From a computational viewpoint, to make the local time step method more efficient in a parallel context, a new Block-Based AMR technique (BB-AMR) is applied. A series of 1D, 2D and 3D test cases have been performed to test and validate our approach, using several model configurations (uniform or adaptative grids, Adams-Bashforth or Runge-Kutta schemes, etc. ). Moreover, the BB-AMR method is improved using an Automatic Mesh Refinement Threshold based on the decreasing rearrangement of the mesh refinement criterion function yielding to a robust parameterless AMR method.

From an application viewpoint, we have successfully applied and solved several problems [START_REF] Golay | Block-based adaptive mesh refinement scheme using numerical density of entropy production for three-dimensional two-fluid flows[END_REF][START_REF] Pons | Adaptive mesh refinement method. Part 2: Application to tsunamis propagation[END_REF][START_REF] Yushchenko | Entropy production and mesh refinement -Application to wave breaking[END_REF]: dam-break, breaking waves, run-up, reflection and/or flooding, drying phenomena, etc.

We are now interested in the extension of the method to onedimensional and two-dimensional dispersive equations for dispersive waves for tsunami propagation for instance.

Part IV C F D : D G M E T H O D , h p -A M R A N D A P P L I C AT I O N S T O WAV E -D R I V E N G R O U N D WAT E R F L O W S
Les vagues ne naissent pas quand elles déferlent sur la plage.

Elles roulent longtemps en mer avant de se casser sur un repli de sable.

-Geneviève Dormann

In this part, I present the work in collaboration with J.B. Clément, F. Golay and D. Sous for wave-driven groundwater flows in sandy beaches. This work is a part of the PhD of J-B. Clément (see [START_REF] Clément | DISCONTINUOUS GALERKIN METHOD FOR STEADY-STATE RICHARDS EQUA-TION[END_REF][START_REF] Clément | ADAPTIVE DISCON-TINUOUS GALERKIN METHOD FOR RICHARDS EQUATION[END_REF][START_REF] Clément | A Richards's equation based model for wave-resolving variable-saturation groundwater beach dynamics[END_REF][START_REF] Clément | An adaptive strategy for discontinuous Galerkin simulations of Richards' equation[END_REF][START_REF] Clément | Wave-driven Groundwater Flows in Sandy Beaches: A Richards Equation-based Model[END_REF]).

Partially-saturated porous media are described by Richards' equation. Simulations of Richards' equation involve sharp wetting fronts which evolve both in space and in time as well as unsaturated/saturated regions of parabolic/elliptic natures which stands alongside. Getting numerical solutions for Richards' equation can be hard to achieve. In this work, a DG method is chosen since they share advantages both with finite elements methods and finite volumes methods (as developed in Part iii). It allows in principle to reach high-order accuracy schemes, h-adaptivity, p-adaptivity and more generally hp-adaptivity.

I N T R O D U C T I O N

This work enters into a long-term research project aiming to simulate the interactions between waves and groundwater in sandy beaches. Recent studies in laboratory [START_REF] Steenhauer | Subsurface processes generated by bore-driven swash on coarsegrained beaches[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] and field studies [START_REF] Heiss | Coupled surface-subsurface hydrologic measurements reveal infiltration, recharge, and discharge dynamics across the swash zone of a sandy beach[END_REF][START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] provided a good experimental understanding of groundwater dynamics and shown a circulation pattern below the swash zone. Up to now, few models have been suggested and their scope remains limited [START_REF] Li | BeachWin: modelling groundwater effects on swash sediment transport and beach profile changes[END_REF][START_REF] Malott | Dynamic groundwater flows and geochemistry in a sandy nearshore aquifer over a wave event[END_REF].

New numerical developments are necessary to capture a maximum of physical processes, in particular to resolve all along the wave phase the complex dynamics of saturation and pressure fields. The modelling is quite hard to achieve because it requires to solve both nonlinear equations for subsurface water flow and the shallow water equations for surface flow. Additionally, these models has to be coupled in a suitable way in order to respect the physics of the problem. Indeed, on one hand, the dynamics of incoming waves is very fast due to wave run-up and run-down which are covering and uncovering alternatively the beach. But, on the other hand, groundwater circulation is particularly slow [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] because of sediments acting as low-pass filter. Thereby, infiltrations/exfiltrations of the swash zone, sand saturation fluctuations, multiple timescales and various space-scales require a robust and accurate methodology for both numerical schemes and algorithms as well as a careful insight into the modelling.

In this work, we focus on the subsurface flow which is assumed to occur in a partially saturated porous medium and can be described by Richards equation. The latter takes into account the actions of gravity (advection) and capillarity (diffusion) but neglecting the flow of the non-wetting phase, namely the air [START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media[END_REF]. Richards equation has been widely used to simulate water flow in unsaturated porous media [START_REF] Belfort | On equivalent hydraulic conductivity for oscillation-free solutions of Richard's equation[END_REF][START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF][START_REF] Sochala | Mass conservative BDF-discontinuous Galerkin/explicit finite volume schemes for coupling subsurface and overland flows[END_REF]. Nevertheless, its derivation is generally rather roughly introduced, which leads us to first set out the main steps to underline some issues in regards to our model problem. Richards equation is a nonlinear parabolic equation but it degenerates into an elliptic equation under saturated conditions. Nonlinearities are due to constitutive laws and, consequently, the solution includes sharp moving fronts with dynamic smooth regions which challenge numerical methods to obtain robust and efficient solutions of Richards equation [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF]. Extensive efforts have been made to overcome these difficulties, in particular to find efficient nonlinear solvers for Richards equation [START_REF] Cooley | A Finite Difference Method for Unsteady Flow in Variably Saturated Porous Media: Application to a Single Pumping Well[END_REF][START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF][START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF][START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF][START_REF] List | A study on iterative methods for solving Richards' equation[END_REF][START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF][START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF][START_REF] Rubin | Theoretical Analysis of Two-Dimensional, Transient Flow of Water in Unsaturated and Partly Unsaturated Soils1[END_REF][START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF]. Results of existence, regularity and uniqueness were found [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] and solution was assessed through diverse simulations involving different numerical methods.

Current research into Richards equation consist in advance spatial and temporal discretizations [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF] allowing construction of adaptive approximations in size and/or in order for space and time. A primary focus on spatial discretization can be outlined [START_REF] Baron | Adaptive multistep time discretization and linearization based on a posteriori error estimates for the Richards equation[END_REF][START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF][START_REF] Sochala | Mass conservative BDF-discontinuous Galerkin/explicit finite volume schemes for coupling subsurface and overland flows[END_REF] with the use of finite differences, finite volumes, finite elements, mixed finite elements and DG. We choose to use DG method because it offers suitable advantages for adaptivity and the treatment of flow dynamics [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF]. Indeed, DG methods may be seen as a method sharing properties both from finite elements and finite volumes : they are based on a variational formulation but in an element-wise fashion. So they are locally conservative which is crucial in fluid dynamics and explains why the method is more and more used for problems in porous media. Moreover, it enables to change the degree of polynomial approximation and to use non-conforming mesh (hanging node). This is an important benefit since it is possible to handle high-order accuracy and hpadaptation. Various forms of DG methods [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF] can be used for Richards equation which can be discretized in a primal or mixed formulation.

Here, we will focus on primal interior penalty DG methods.

richards equation and constitutive laws

Richards equation is a classic nonlinear parabolic equation to describe flow occurring in the unsaturated zone of an aquifer [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF][START_REF] Szymkiewicz | Modelling Water Flow in Unsaturated Porous Media[END_REF]:

∂ t (θ (ψ)) -∇ • (K (ψ) ∇ (ψ + z)) = 0 (10.1)
where ψ denotes the pressure head (m), K the hydraulic conductivity tensor (m/s), z the elevation (m) and θ the water content. Equation (10.1) is called the mixed form of Richards equation because it involves both variables, θ and ψ. This form was chosen because it is mass conservative and is applicable to a wide range of porous media [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF]. Throughout this work, hydraulic head h (m) and flux q (m/s) are also used: h = ψ + z and q = -K (ψ) ∇h. Equation (10.1) is rewritten as:

∂ t (θ (h -z)) -∇ • (K (h -z) ∇h) = 0 . (10.2)
To close this system of equations, two constitutive laws are needed for water content and hydraulic conductivity. To describe these relations, it exists several models which depends on hydraulic properties of the porous medium. In this work, the model used for the simulations [START_REF] Haverkamp | A Comparison of Numerical Simulation Models For One-Dimensional Infiltration1[END_REF][START_REF] Vauclin | Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils[END_REF] is:

(K(ψ); θ(ψ)) =                  AK s A + |ψ| B ; C (θ s -θ r ) C + |ψ| D + θ r , if ψ ≤ 0 (unsaturated zone) (K s ; θ s ) , if ψ > 0 (saturated zone) (10.3)
where A, B, C and D are empirical parameters, K s is the hydraulic conductivity at saturation, θ s is the saturated water content and θ s is the residual water content.

seepage boundary condition

The seepage boundary condition is specific to subsurface model [START_REF] Scudeler | Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models[END_REF]. It models the interface between porous medium and atmosphere. If the porous medium is saturated and outflow occurs, then water pours out at atmospheric pressure (ψ = 0). Otherwise, the interface acts as an impervious boundary. The difficulty is that the length of seepage face is unknown a priori. It can be seen as a non-linear Robin boundary condition:

a (h) h -b (h) K (h -z) ∇h • n = c (h) where      a (h) = H (h) b (h) = 1 -H (h) c (h) = H (h) z with H (h) = 1 if "seepage", 0 otherwise.
Sometimes, the Heaviside function H is regularized [START_REF] Dolejší | Adaptive higher-order spacetime discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF] to overcome convergence difficulties. The seepage boundary condition is also taken into account as a Signorini's type:

ψ ≤ 0, -K (h -z) ∇h • n ≥ 0, (-K (h -z) ∇ (h) • n) ψ = 0 .
In this work, the seepage boundary condition is not regularized and is treated like in [START_REF] Scudeler | Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models[END_REF] by switching between Dirichlet and Neumann boundary condition for each degree of freedom within the weak formulation of the next section. Therefore, if a seepage boundary condition lies on Γ S , it is written:

h = z if ψ ≥ 0 and -K (h -z) ∇h • n > 0 , -K (h -z) ∇h • n = 0 otherwise.
Paniconi and Putti [START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF] reported that seepage condition can cause oscillations leading to convergence problems if rapid changes of saturation state occur near the seepage face, no matter convergence requirements were relaxed or not.

D G M E T H O D

discretization

First of all, some definitions are set for the DG formulation. More careful developments can be found in Dolejší and Feistauer [START_REF] Dolejší | Discontinuous Galerkin Method[END_REF] or Rivière [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF].

Let Ω be a space of R 2 with ∂Ω = Γ D ∪ Γ N if no seepage boundary condition is considered. From equation (10.2) and by adding initial and boundary conditions, the model problem is: Let be two neighbouring elements E l and E r sharing one face F. There are two traces of a function u along F denoted by u l and u r . The jump across a face is denoted by u = u lu r , the mean by {u} = 1 2 (u l + u r ) and n is the unit normal vector oriented from E l to E r .

           ∂ t θ (h -z) -∇ • (K (h -z) ∇h) = 0, in Ω × [0, T] h = h D , on Γ D × [0, T] -K (h -z) ∇h • n = q N , on Γ N × [0, T] h = h 0 , in Ω × {0}
Several DG methods are possible. Here, the incomplete interior penalty Galerkin (IIPG) method is chosen to avoid the addition of stabilization terms from symmetrization [START_REF] Dolejší | Discontinuous Galerkin Method[END_REF]. Penalization terms are added to enforce Dirichlet boundary condition and continuity constraint through respectively, the boundary penalty σ b and the interface penalty σ i . Problem (11.1) is multiplied by a test function v, integrated on each element E. After using Green's theorem, summing over all ele-dg method ments and performing some calculations, the space semi-discretization is obtained:

s > 3 2 , Find h ∈ H s (E h ) such that ∀v ∈ H s (E h ) , ∑ E∈E h E ∂ t θ (h -z) v dx + a h (h, v) = l h (v) with (11.2) a h (h, v) = ∑ E∈E h E K (h -z) ∇h • ∇v dx -∑ F∈F I h F {K (h -z) ∇h • n F } v ds + ∑ F∈F I h F σ F i |F| h v ds -∑ F∈F D h F K (h -z) (∇h • n F ) v ds + ∑ F∈F D h F σ F b |F| hv ds l h (v) = ∑ F∈F D h F σ F b |F| h D v ds -∑ F∈F N h F q N v ds.
Contrary to the Finite Elements method, there is no continuity constraint between elements for the test functions. Therefore, any polynomial basis functions is suitable what makes DG method flexible. Even if a nodal basis can be chosen, this work is based on a modal basis, the monomial basis, because it is simple and can be extend easily to higher orders.

To get the full discretization, the first term of the weak formulation (11.2) is approximated by a backward Euler scheme which gives:

Find h ∈ H s (E h ) such that ∀v ∈ H s (E h ) , ∑ E∈E h E θ(h n+1 -z)-θ(h n -z) δt v dx + a h (h n+1 , v) -l h v; t n+1 = 0 ⇐⇒ R(h n+1 , v) = 0 (11.3)
where δt is the time step and n is the number of time iteration. The choice is quite natural for backward Euler scheme because its implicit nature provides natural stability properties. Good convergence behaviour for Richards equation was also observed in [START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF]. Together with its simple implementation, this scheme is widespread and use of higher order scheme is questionable.

nonlinear solver with adaptive time stepping

Solving equation (11.3), which represents the residual R, can be challenging because of strong nonlinearities. Newton-Raphson methods converge quadratically but they often fail for parabolic degenerate problem where the Jacobian matrix might become singular and because of poor initializations [START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF][START_REF] List | A study on iterative methods for solving Richards' equation[END_REF][START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF]]. That's why fixed-point iteration methods are robust and popular techniques for Richards equation even if they converge linearly1 .

Let m be the number of nonlinear iteration and J is either the Jacobian matrix computed from (11.3) or is the matrix K(h) coming from the variational formulation (11.3) written under a matrix system K(h)h = f . These two schemes reads:

   J h n+1,m δh = -R h n+1,m h n+1,m+1 = h n+1,m + δh
In this study, the fixed-point iteration method is chosen and constructed like in [START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF] by avoiding the first derivative terms inside the Jacobian matrix. The stopping criterion is ||R(h)|| L 2 (Ω) < ε where ε is a user-defined tolerance. Time adaptation is motivated by the convergence of the nonlinear solver which often failed under too strong boundary and initial conditions. On one hand, transient simulations have difficulties to converge if the time step is too large but, on the other hand, shorter time steps mean more time steps and so, a longer computational time. That is why time adaptation is very attractive and common for Richards equation. However, it exists different strategies to adjust the time step [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF]. First ones are heuristic and mainly based on convergence performance of the nonlinear solver. Second ones are rational and based on error control. The latter are generally much more efficient than the former.

Yet, heuristic methods may outperform them in regards of their simple implementation.

For this study, the time step is adjusted heuristically according to the previous number of iterations N it from the nonlinear solver such as [START_REF] Bergamaschi | Mixed finite elements and Newtontype linearizations for the solution of Richards' equation[END_REF][START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF]. The simulations begins with a time step δt 0 . The next time step δt n+1 is calculated from the previous one δt n according to the time-stepping scheme:

                         δt n+1 =            λ amp δt n if N it ≤ inf lt δt n if inf lt < N it ≤ sup lt λ red δt n if sup lt < N it ≤ max lt δt n = λ red δt n if max lt < N it
or if the solver has failed (time step is started again)

The factors 0 < λ red < 1 and 1 < λ amp as well as the threshold values 0 < inf lt ≤ sup lt < max lt are prescribed by the user. A minimum time step is defined to avoid excessive small time step. With this approach, the nonlinear solver is more robust and computational time is saved compared to a fixed time step. Indeed, whatever δt 0 , the time step is adjusted until the nonlinear solver succeeds to converge. Then, thanks to the coefficients, the time step is adjusted smoothly. Nevertheless, the method depends on many empirical parameters and does not provide an optimal time step because the amplification/reduction coefficients restrict the possible size of the new time step.

adaptive mesh refinement

Solution of Richards equation holds sharp wetting fronts which are moving dynamically through the domain among smooth regions [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF]. They arise typically from boundary conditions and dissipate over time by coming across the saturated zone. Mesh adaptation is a technique which can handle these particularities to reach the desired accuracy [START_REF] Dolejší | Adaptive higher-order spacetime discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF][START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF][START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF]. Recent strategies involve both size and order adaptation, based on a posteriori errors [START_REF] Li | Adaptive local discontinuous Galerkin approximation to Richards' equation[END_REF] or interpolation error [START_REF] Dolejší | Adaptive higher-order spacetime discontinuous Galerkin method for the computer simulation of variably-saturated porous media flows[END_REF].

In this study, an h-adaptation based on a simple a posteriori error is performed like [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF]. Yet, whereas they use changes in effective saturation [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF], that's directly changes in the solution, the hydraulic head h, which serve as error indicator. The latter has the advantage of working both for the saturated and unsaturated regions.

Mesh adaptation is a 2D block-based technique adapted from [START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF] and Part iii. The process of mesh adaptation begins with the evaluation of the criterion defined by:

C E = 1 max E∈E h 1 |E| ||∇h|| L 2 (E) 1 |E| ||∇h|| L 2 (E)
The criterion is normalized so that 0 ≤ C E ≤ 1. The second step consist in deciding to refine, coarsen or do nothing thanks to user-defined threshold values 0 < β min ≤ β max < 1 for the criterion. Elements are merged or split according to the scheme described in Figure 11.1. To avoid excessive refinement and coarsening, minimum and maximum levels are set up. Then, the mesh is smoothed by checking that each neighbouring block of elements cannot have a difference of refinement level superior to one as done in Part iii. In the last step, the solution is projected from the previous mesh to the new one by solving a local DG problem. The element-wise weak formulation corresponding to h E new = h E old reads:

E new h E new v dx = E new h E old v dx
This strategy of adaptation is very simple but it uses an a posteriori indicator and requires empiric parameters which balance computational performance and accuracy. This forces to calibrate each numerical simulation.

N U M E R I C A L R E S U LT S

one-dimensional infiltration numerical test case

This numerical test-case was firstly introduced by [START_REF] Celia | A general massconservative numerical solution for the unsaturated flow equation[END_REF] on the basis of an experiment from [START_REF] Haverkamp | A Comparison of Numerical Simulation Models For One-Dimensional Infiltration1[END_REF]. It was taken afterwards by others [START_REF] Manzini | Mass-conservative finite volume methods on 2-D unstructured grids for the Richards' equation[END_REF] and corresponds to a class of popular test cases, for example, see [START_REF] Miller | A spatially and temporally adaptive solution of Richards' equation[END_REF].

This problem deals with an 1D-infiltration into a soil column. K and θ are taken from (10.3) with A = 1.175 × 10 6 , B = 4.74, C = 1.611 × 10 6 , D = 3.96, K s = 0.00944 cm/s, θ s = 0.287 and θ r = 0.075. The initial condition is h 0 = -61.5 + z cm. On the top and bottom, hydraulic head is kept fixed: h top = 19.3 cm, h bottom = -61.5 cm. Sides are impervious. Initial condition is: h = -61.5 + z cm.

The solutions are in very good agreement to one another as seen for pressure head in the Figure 12.1a. Moreover, h-adaptation is capturing the wetting front efficiently as it is shown in Figure 12.1b.

comparison with water table drainage experiment

Vauclin et al. made a drainage experiment of a slab of soil in which experimental details and positions of water table can be found [START_REF] Vauclin | Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils[END_REF].

A water table depletion occurs after an instantaneous drop of the imposed hydraulic head. K and θ are taken from (10.3) with A = 359720.49, B = 4.5, C = 40000, D = 2.9, K s = 40 cm/h, θ s = 0.3 and θ r = 0. This experiment is used as a test-case to show a transient simulation with seepage boundary condition. The domain is initially in hydrostatic equilibrium: h 0 = 145 cm. Boundary conditions on the right, at the bottom, on the top and on the left for z > 145 cm are ∇h • n = 0. On the left, the level is kept constant in the ditch for z ≤ 75 cm, h = 75 cm, and, the seepage boundary condition is prescribed for 45 cm < z ≤ 145 cm.

Numerical simulation is performed with linear approximation p = 1, σ i = σ b = 100 on a mesh with ∆x = 10 cm and ∆z = 5 cm. Water table location matches ψ = 0. Numerical positions of water table are displayed in Figure 12.2 together with the experimental data from [START_REF] Vauclin | Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils[END_REF]. Results presented in the Figure 12.2 show that the model is able to reproduce satisfactory the transient water table dynamics. Disagreements near the seepage boundary are found between numerical and experimental profiles. Three reasons can explain it. One is the lack of spatial resolution. Refining the grid would improve the computed solution and catch better the nonlinear seepage boundary condition. Another reason is in relation with the experimental set-up because water should not pour out completely freely through the experimental device in such a way that seepage boundary condition does not model it properly. Finally, the determination of hydraulic properties was done by correlation from measurements with some fairly important dispersion, due to natural heterogeneity, neglecting hysteresis [START_REF] Vauclin | Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils[END_REF]. This can explain some discrepancies. Other numerical studies based on this experiment are reported [START_REF] Clement | A physically based, twodimensional, finite-difference algorithm for modeling variably saturated flow[END_REF][START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF]. They also struggle to match water table profiles near the seepage.

idealized beach test case

In this test case, the Simulation Waves till SHores (SWASH) code is used to simulate waves in coastal waters and to produce data which enforce Richards equation through suitable boundary conditions. This weak coupling can take into account seepage along the beach face. This idealized case is constructed to investigate the beach groundwater response to the action of swash: infiltration/exfiltration, water table position, hydraulic head and saturation distribution as well as pore velocities are outlined.

12.3.1 Description of the idealized beach case Initial conditions are hydrostatic with water table level at same height of surface water: h 0 = 1.5 m. The bottom is an impervious boundary. The water table and head are kept fixed on the inland boundary: h = 1.5 m. The beach face is monitored by the SWASH code. In SWASH, a sinusoidal time-series of free surface elevation is imposed at the offshore boundary to simulate a long infragravity-like wave of period 50 s and amplitude 0.3 m. These types of waves have been observed to be a major driver of beach groundwater dynamics of the Rousty microtidal dissipative beach [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF]. This results in two swash events on the beachface. Time simulated lasts 120 s.

Two waves are generated by SWASH and their characteristics are: wavelength, frequency. This results in two swash events on the beach face. Time simulated lasts 7.2 s. The mesh used for the simulation is displayed on the Figure 12.3.

Results

An overview of groundwater features during a swash cycle is presented. Then, results focus on uprush and backwash events. Some selected snapshots from the simulation show interesting observations for the idealized beach case. Surface water and water table are depicted by blue and white solid lines respectively. Hydraulic head contours are displayed together with velocities vectors.

overall pattern for swash cycle:

A first observation concerns wave propagation. When a wave is moving, groundwater flows in the saturated zone respond immediately. Hydraulic head spreads vertically quite sharply under the wave front. Infiltration/exfiltration occurs at the back/front of the bore respectively. It is particularly visible in the Figure 12.4 even if this phase is completely ideal because it's about the first wave moving through hydrostatic conditions.

Throughout the simulation, a nearly hydrostatic state is observed in the upper beach area. Moreover, flow divides into seaward and landward components near the uprush limit. Velocities magnitude is rather small everywhere, about 1 cm/h, except for uprush infiltration at the top where it reaches up to 160 cm/h. Around exit point of water table, a local overpressure is observed. It is probably caused by numerical artefacts due to the lack of mesh resolution for the seepage boundary condition. An overall pattern can be depicted for swash cycles. Two recirculation zones are attached respectively to the back and the front of the swash tongue. These two groundwater cells are interacting together and are evolving during uprush and backwash, growing or diminishing. Sometimes one vanishes completely, sometimes another one gets trapped by recharge induced by water table elevation. pattern during uprush event:

During uprush, there is downward infiltration while the swash tongue is covering the beach. It is mainly concentrated in the upper zone, through the boundary in contact with the fully saturated part of the beach. Infiltration induces a rise of the water table near the exit point, with, as seen in Figure 12.5.

pattern during backwash event:

During backwash, a low hydraulic head zone emerges at the base of the swash zone where exfiltration occurs. Velocities are mainly seaward under the swash tongue moving back. Meanwhile, a small region of high hydraulic head persists in the corner of the saturated zone where water table have been raised by the precedent uprush. It results in recharge and, as water is percolating inside the porous medium, flow divides between seaward and landward direction. Moreover, above the unsaturated zone which has been covered by uprush, water is infiltrating downward in a sharp wetting front. While the swash toe retreating to the sea, a seepage face is appearing which drops, locally but quickly, the water table what gives form to a small hump-like structure. Figure 12.6 presents these observations. pattern between backwash and uprush events:

The next incoming wave arrives just after the backwash. It pushes the region of low hydraulic head up. As this zone swept away, it reaches the remaining small region of high hydraulic head induced by the water table elevation and makes it disappear temporally. The incoming wave is bringing a new high hydraulic head front. This causes infiltration which arises water table again, ending seepage and erasing the hump-like structure of water table. Then, a new cycle begins. This intermediate state is shown in Figure 12.7. pattern between backwash and uprush events:

The next incoming wave arrives just after the backwash. It pushes the low hydraulic head region up. As it swept away, the zone concentrates more and more and reaches the remaining small high hydraulic head region around hump-like structure to make it disappear. The wave which is coming goes with a new high hydraulic head front. This causes infiltration which arises again water table, ending seepage and erasing the hum-like structure. The situation is shown in Figure 12.7.

Discussion

In the simulation, the offshore water level is similar to the inland water table, therefore no steady pressure gradient is observed across the beach-face by contrast to field or laboratory studies where crossbeach gradients can be observed depending on boundary conditions [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. Within the beach, overall velocity magnitude is about 1 cm/h which is fairly comparable to measurements made in [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF]. The only exception is for infiltration/exfiltration where maximum velocity order is about 100 cm/h. This seems to be less compared with the experiments of [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. Further tests are now engaged in the exact configurations of the laboratory [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] and field [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] experiments to allow direct comparisons. As described by [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF], the moving bore is carrying a sharp hydraulic front. Exfiltration occurs at its front and infiltration at its back when it meets the groundwater conditions of the precedent swash. A local superelevation of the water table has also been revealed in these observational studies due to the water table charge by short or IG wave forcing under the swash zone. [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF] showed a domed water table in their study. The present benchmark case will be extended to a longer duration to test if the model is able to represent such feature. Water table elevation induces charge in this zone as long as the low-pressure zone does not affect it with the ongoing uprush. As a consequence, infiltration with flow division occurs near the uprush limit. Laboratory experiments and numerical simulations [START_REF] Li | Wave-induced beach groundwater flow[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] both predict this feature but [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] do not observe it. During the simulation, seepage face has small length because water is able to pour out from sand quite fast following surface water. As a consequence, a small hump-like structure appears during backwash because of seepage along the boundary. This last feature is not described in any study found, but further tests are ongoing with more realistic configurations and longer wave forcings to better understand involved processes. The two circulation zones with infiltration/exfiltration in the upper/lower part of the beach are also identified in the simulation from [START_REF] Li | Wave-induced beach groundwater flow[END_REF] and experiments from [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. Infiltration and exfiltration involve greater pore velocities than in other parts of the beach. Like it is depicted by [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF], gradients brought by the wave are much larger than the average values in the beach due to larger differences in water depth across the wave.

In this work, we have developed a full numerical code on unstructured meshes for Richards equation through a DG.

This numerical code, associated with a backward Euler scheme, is able to simulate variably saturated flows in porous media for a wide range of problems steady and unsteady problems, including seepage boundary condition. Richards equation gives rise to some issues concerning non-linearities and convergence robustness. In order to provide a computationally simple but efficient algorithm in terms of solution accuracy and robustness, the nonlinear solver is based on a fixed-point iteration method with adaptive time stepping. Besides, an adaptive mesh refinement technique is employed to achieve accuracy and deal with steep gradients. Through several test cases ([35, 36, 39]) numerical results present the functioning of the method and its potential.

Further research is underway to extend adaptive mesh refinement to an efficient hp adaptation, assess more carefully the convergence behaviour of seepage boundary condition and make time step adaptation less heuristic. -

Mehmet Ersoy

In this part, I focus on the model reduction in the shallow water limit for compressible and incompressible flows with hydrostatic and non-hydrostatic pressure distribution. Opposite to the CFD viewpoint (see Part iii and Part iv), the idea in this part is to reduce the dimension of the initial model preserving some of the main physical and geometrical features. This is one of the most challenging issues that we address with the obvious consequence to decrease the computational time.

In this part, in collaboration with M.A. Debyaoui, we present the first non-linear and weakly dispersive extension of the FS-model [START_REF] Ersoy | A Construction of Biorthogonal Wavelets With a Compact Operator[END_REF] (or the Serre Green Naghdi (SGN) equations) for natural rivers and channel hydraulics ( [START_REF] Debyaoui | Generalised Serre-Green-Naghdi equations for open channel and for natural river hydraulics[END_REF]). This work is a part of the PhD of M.A. Debyaoui (see also [START_REF] Debyaoui | A Generalised Serre-Green-Naghdi equations for variable rectangular open channel hydraulics and its finite volume approximation[END_REF][START_REF] Debyaoui | Generalised Serre-Green-Naghdi equations for open channel and for natural river hydraulics[END_REF]).

In Chapter 14, we first motivate our work. Then, in Chapter 15, we present the geometrical settings and the Euler equations. Chapter 16 is devoted to the asymptotic limit of the Euler equations. In Chapter 17, we present the new model for arbitrary non rectangular channel/river. Finally, the rectangular case is developed in Chapter 18 and a finite numerical scheme is proposed.

I also refer to [START_REF] Ersoy | Dimension reduction for incompressible pipe and open channel flow including friction[END_REF][START_REF] Ersoy | Dimension reduction for compressible pipe flows including friction[END_REF] for a rigorous model reduction including friction and to [START_REF] Ersoy | A Saint-Venant shallow water model for overland flows with precipitation and recharge[END_REF] for overland flows including rainfall and infiltration (or exfiltration). These works are not presented here.

I N T R O D U C T I O N
The modelling of hydrology of catchment basins and rivers holds a central place in environmental sciences, particularly in connection with water availability, urban sewer systems, flood risks and in particular for tsunamis. Indeed, rivers are known to be the tsunami highways.

Waves penetrate through rivers much faster inland than the coastal inundation reaches over the ground, and may lead flooding in lowlying areas located several km away from the coastline [START_REF] Tsuji | Tsunami ascending in rivers as an undular bore[END_REF]. This is important today in understanding and forecasting the impact of climate variability on the human and natural environment. Modelling these processes and predicting the motion of water is a difficult task for which substantial effort has been devoted [START_REF] Esteves | Overland flow and infiltration modelling for small plots during unsteady rain: numerical results versus observed values[END_REF][START_REF] Grace | The modeling of overland flow[END_REF][START_REF] Graf | Hydraulique fluviale: écoulement et phénomènes de transport dans les canaux à géométrie simple[END_REF][START_REF] Rousseau | Study of overland flow with uncertain infiltration using stochastic tools[END_REF][START_REF] Singh | Numerical Modeling of Rainfall-Generated Overland Flow Using Nonlinear Shallow-Water Equations[END_REF][START_REF] Weill | A generalized Richards equation for surface/subsurface flow modelling[END_REF][START_REF] Woolhiser | Unsteady, one-dimensional flow over a plane -The rising hydrograph[END_REF][START_REF] Zhang | Modeling of two-dimensional overland flow[END_REF].

One of the most widely used models to describe the channel and river motion of watercourses is the section-averaged free surface model [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Decoene | Asymptotic derivation of the section-averaged shallow water equations for natural river hydraulics[END_REF][START_REF] Ersoy | Dimension reduction for incompressible pipe and open channel flow including friction[END_REF] which is a generalisation of the well-known Saint-Venant system (introduced by Adhémar Jean Claude Barré de Saint-Venant in the 19th Century [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF][START_REF] De Saint-Venant | Théorie et equations générales du mouvement non permanent des eaux courantes[END_REF]):

     ∂ t A + ∂ x Q = 0 , ∂ t Q + ∂ x Q 2 A + I 1 (x, A) = I 2 (x, A) . ( 14.1)
In these equations, A is the wet area of fluid cross-section, Q is the water discharge, I 1 is the hydrostatic pressure and I 2 is the hydrostatic pressure source term which takes in account of the variation of the section. The model (14.1) reduces to the well-known one-dimensional Saint-Venant equations for uniform rectangular section. The free surface model is the first order shallow water approximation of the sectionaveraged Navier-Stokes or Euler equations under suitable assumptions on the horizontal and the vertical scales (see, e. g. , [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Decoene | Asymptotic derivation of the section-averaged shallow water equations for natural river hydraulics[END_REF][START_REF] Ersoy | Dimension reduction for incompressible pipe and open channel flow including friction[END_REF][START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF][START_REF] Gouta | A finite volume solver for 1D shallowwater equations applied to an actual river[END_REF] and the reference therein).

Thanks to the hyperbolic structure of these equations, sharp transitions between two different flow states result in a discontinuous solution, both in the water surface and in the velocity. These discontinuous solutions (called shocks and also referred as bores) are well-suited to approximate breaking waves with turbulent rollers for large transitions of the Froude's number. However, for small or moderate transitions, the advancing wave front can be followed by a train of free-surface undulations, sometimes called "whelps". This phenomenon, called undular bore (also often called dispersive shock waves), is induced by a non hydrostatic pressure distribution [START_REF] Peregrine | Calculations of the development of an undular bore[END_REF]. As a consequence, wave solutions spread out in space as they evolve in time, i. e. waves of different wavelengths travel with different speeds. This is the so-called dispersive effect. Consequently, undular bores are not reproductible with the non-dispersive free surface system and non hydrostatic pressure is required.

Dispersive equations were first introduced by Boussinesq [START_REF] Boussinesq | Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire[END_REF] in 1872 to mathematically justify the existence of solitary waves observed by Russell's experiments in 1834. These equations enters in the framework of shallow water equations. They can be obtained as the second order asymptotic approximation in µ 2 , with µ = H L 2 , of the depth-averaged Euler equations where H represents a characteristic water depth and L is a characteristic horizontal scale.

The Boussinesq type equations, which are weakly non-linear and weakly dispersive, are defined by an additional small non-linearity parameter ε = a H = O(µ) 1 where a is the order of the free surface amplitude. These equations for a flat bottom, for instance in 1d, are given by

∂ t ξ + ∂ x (hu) = O(µ 2 ) ∂ t u + εu∂ x u + ∇ξ + µD = O(µ 2 ) (14.2)
where h is the water depth, ξ is the free surface elevation and u is the depth-averaged velocity. The term D represents the dispersive term. In 1877, the KdV equation was discovered by Boussinesq [START_REF] Boussinesq | Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire[END_REF] and was later derived by Korteweg and Gustav de Vries (KdV) [START_REF] Korteweg | XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF]. This equation approximately describes the evolution of long, onedimensional waves in many physical settings, including shallow-water waves with weakly non-linear behaviour. This equation can also be obtained in the case of unidirectional waves for which the theoretical framework allows to compute an analytical solution such as 1d solitary wave propagation on flat bottom (see e.g. Boussinesq [START_REF] Bayraktar | A new numerical scheme for improved Boussinesq equations with surface pressure[END_REF][START_REF] Boussinesq | Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire[END_REF][START_REF] Korteweg | XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves[END_REF]). We have also the Benjamin, Bona, and Mahony (BBM) equation which is an improvement of the KdV equation [START_REF] Benjamin | The stability of solitary waves[END_REF]. However, in practice and especially in coastal engineering applications, the nearshore wave dynamics being often varying-bottom dependent, dispersive and nonlinear, the Boussinesq type equations are not appropriate.

In 1967, Peregrine [START_REF] Peregrine | Long waves on a beach[END_REF] introduced the first weakly non-linear two dimensional Boussinesq type equations for non flat bottom. Witting [START_REF] Witting | A unified model for the evolution nonlinear water waves[END_REF] proposed a method based on Padé expansion to improve the frequency dispersion of the Boussinesq-type equations. From this method, several equations of order O(µ 2 ) with improved dispersion characteristic have been proposed, see for instance [START_REF] Madsen | A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF][START_REF] Schäffer | Further enhancements of Boussinesqtype equations[END_REF].

In 1953, a one-dimensional, fully non-linear (ε = O(1)) and weakly dispersive equations for flat bottom were derived by Serre [START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF], independently of Su and Gardner [START_REF] Su | Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries equation and Burgers equation[END_REF] with

D(v) = 1 h ∂ x h 3 3 D(v) (14.3)
where

D(v) = (∂ x v) 2 -∂ x ∂ t v -v∂ xx v . ( 14.4)
This model was then extended to non flat bottom by Seabra-Santos et al [START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF], see also [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: boundary conditions and validation[END_REF]. Finally, Green and Naghdi [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF] derived the two dimensional fully nonlinear dispersive equations for uneven bottom which are the extension of Serre equations. In the literature, this system is often called SGN equations. Further interesting extensions are also presented in [START_REF] Clamond | Conservative modified Serre-Green-Naghdi equations with improved dispersion characteristics[END_REF][START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF][START_REF] Pons | Modélisation des tsunamis : Propagation et impact[END_REF][START_REF] Richard | Modelling turbulence generation in solitary waves on shear shallow water flows[END_REF][START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF] based on either Euler equations, water waves equations or variational principles. All of the previous dispersive models are obtained either from 3d-2d or 2d-1d reduction, but, up to our knowledge, the 3d-1d reduction has never been done before. Thus, our main goal is to derive from the three-dimensional incompressible and irrotational Euler equations with suitable boundary conditions, a model akin to (14.2) -(14.4) via section averaging under the shallow water assumption. The section averaged model that we obtain extends the section-averaged free surface model (14.1) and the SGN equations (14.2) -(14.4):

         ∂ t A + ∂ x Q = 0 ∂ t Q + ∂ x Q 2 A + I 1 (x, A) + µ 2 ∂ x (G(x, A)D(u)) = I 2 (x, A) + G(x, A, Q) + O(µ 2
2 ) where A is the wet area, Q = Au is the water discharge, u is the sectionaveraged velocity, I 1 (resp. I 2 ) is the hydrostatic pressure (resp. source) term, G(x, A) generalises h 3 3 in (14.3), the term G(x, A, Q) extends the uneven bottom source term, D is the term given by (14.4) and therefore D = ∂ x (G(x, A)D(u)) represents the dispersive term. where y * describes the transversal variation of the channel with respect to the main channel/river direction.

We define the local height of the water by h(t, x, y) := η(t, x, y)d(x, y).

The boundary of the section can be also described as a function of y by

∀x ∈ [0, L c ], ϕ(x, z) = α(x, z) if ϕ(x, z) < y * (x) β(x, z) otherwise.
(15.1)

the three-dimensional incompressible euler equations

The width of the section at the elevation z is given by an increasing function z → σ(x, z) = β(x, z)α(x, z) .

The wet region is defined as the region in which the fluid resides at each time t ∈ [0, T]

Ω(t) = 0≤x≤L Ω(t, x)
with its global counterpart

Ω := 0≤t≤T Ω(t). (15.2)
Ω(t, x) is the cross-section of fluid at the position x of the channel C:

Ω(t, x) = (y, z) ∈ R 2 ; α(x, z) ≤ y ≤ β(x, z) , d(x, y) ≤ z ≤ η(t, x, y)} (15.

boundary conditions

We assume that the flow is governed, on the space-time domain Ω, by the three-dimensional irrotational and incompressible Euler equations

divρ 0 u = 0, ∂ ∂t (ρ 0 u) + divρ 0 u ⊗ u + ∇p -ρ 0 F = 0 (15.5)
where u = (u, v, w) is the velocity field, ρ 0 is the density of the fluid (taken to be constant), F = (0, 0, -g) is the external force of gravity with constant g and p is the pressure. These equations are completed by the irrotational equation:

curl [u] = 0 . (15.6)

boundary conditions

For a given time t ∈ [0, T] and x ∈ [0, L c ], the boundary of the fluid cross-section Ω(t, x) is composed of a free surface Γ fs (t, x) = {(y, z) ∈ R 2 ; z = η(t, x, y)} and a wet boundary (part of the boundary in contact with water) Γ wb (x) = {(y, z) ∈ R 2 ; z = d(x, y)} such that

∂Ω(t, x) = Γ fs (t, x) ∪ Γ wb (x) .
We prescribe a kinematic boundary condition at the free surface boundary and a no-penetration condition at the wet boundary as described below.

Free surface boundary conditions

Assuming a kinematic boundary condition, we set for all points m(t, x, y) = (x, y, η(t, x, y)) ∈ Γ fs (t, x),

u • n fs = ∂ t m • n fs
where

n fs (t, x, y) = 1 1 + (∂ x η(t, x, y)) 2 + ∂ y η(t, x, y) 2    -∂ x η(t, x, y) -∂ y η(t, x, y) 1   
is the unit outward normal vector to the free surface. This leads to the following explicit form of the kinematic boundary condition

∂ t η + u∂ x η + v∂ y η = w on z = η . ( 15.7)
We also assume that the pressure at the free surface level is equal to the atmospheric pressure p 0 p = p 0 on z = η . (15.8) In the sequel, without loss of generality, we set p 0 = 0.

Wet boundary conditions

On the wet boundary, i. e. , the part of the boundary in contact with water, prescribing a no-penetration condition, we set for all points m(x, y) = (x, y, d(x, y)) ∈ Γ wb (x), u • n wb = 0 which leads to the following explicit form:

u∂ x d + v∂ y d = w on z = d (15.9)
where

n wb (x, y) = 1 1 + (∂ x d(x, y)) 2 + ∂ y d(x, y) 2    ∂ x d(x, y) ∂ y d(x, y) -1   
is the unit outward normal vector to the wet boundary. In view of the definition of the function ϕ, see Equation (15.1), the no-penetration condition can be also expressed as a function of (x, z) by

u∂ x ϕ + w∂ z ϕ = v on y = ϕ (15.10)
where the unit outward normal vector to the wet boundary is

n wb (x, z) = 1 1 + (∂ x ϕ(x, z)) 2 + (∂ z ϕ(x, z)) 2    ∂ x ϕ(x, z) -1 ∂ z ϕ(x, z)    .
In this work, we neglect some physical process arising in river flows: sedimentation, exchange between groundwater flows and subsurface flows, porosities, etc. However, these phenomena can be easily integrated to this work, see for instance [START_REF] Ersoy | Modeling, mathematical and numerical analysis of various compressible or incompressible flows in thin layer[END_REF][START_REF] Ersoy | A Saint-Venant shallow water model for overland flows with precipitation and recharge[END_REF] In this section, we present the strategy to derive the section-averaged non-linear dispersive model (see Chapter 17). To this end, we first introduce in Section 16.1 the dimensionless problem by introducing the classical dispersive parameter µ 2 . We also introduce a dispersive parameter µ 1 but in the transversal direction. In this work, we assume that µ 1 < µ 2 to obtain the section-averaged dispersive model with respect to µ 2 as usually done. Due to the structure of the equations, we cannot obtain the model by a direct section-averaging. We need to develop first suitable asymptotic expansions in two steps, briefly summarized below:

• in Section 16.2, the Euler equations (15.5)- (15.6) are width-averaged to get the following asymptotic expansion of the horizontal fluid velocity u(t, x, y, z) = u (t, x, z) + O(µ 1 )

where u is the width-averaged velocity.

• in Section 16.3, the width-averaged Euler equations allows to obtain the asymptotic expansion of the horizontal width-averaged fluid velocity

u (t, x, z) = u(t, x) + µ 2 f (u(t, x), Ω(t, x)) + O(µ 2 2 )
where ū is the section-averaged velocity for some function f given later on.

Thus, the asymptotic expansion of u up to order O(µ 2 2 ), the so-called "motion by slices" (see [START_REF] Gerbeau | Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation[END_REF]), can be written as follows:

u(t, x, y, z) = u(t, x) + µ 2 f (u(t, x), Ω(t, x)) + O(µ 2
2 ) .

Finally, using these asymptotic expansions, we are able in Chapter 17 to section-average the Euler equations (15.5)- (15.6) to obtain the new one-dimensional non-linear dispersive equations.

dimensionless euler equations

Let us consider the following scales involved in the wave motion: L a characteristic wave-length in the longitudinal direction, H 2 a characteristic water depth, H 1 a characteristic scale of the channel width and h 1 a characteristic wave-length in the transversal direction. We then define the classical dispersive parameter µ 2 (see e. g. [START_REF] Lannes | The water waves problem: mathematical analysis and asymptotics[END_REF])

µ 2 = H 2 2 L 2
and

µ 1 = h 2 1
L 2 where µ 1 is also a dispersive parameter but in the transversal direction.

In the following, we consider the asymptotic regime:

h 1 < H 1 = H 2 L
such that the following inequality holds

µ 1 < µ 2 2 .
Under these assumptions, we get the following ordering:

µ 2 1 < µ 2 1 µ 2 < min µ 2 1 µ 2 2 , µ 1 µ 2 < max µ 2 1 µ 2 2 , µ 1 µ 2 < µ 1 < min µ 1 µ 2 , µ 2 2 < max µ 1 µ 2 , µ 2 2 < µ 2 .
We also introduce U = (U,

V = √ µ 1 U, W = √ µ 2 U) the scale of fluid velocity so that V < W < U. The time scale is T = L U . We set P = p ρ 0 and we define P = U 2 .
This allows us to introduce the dimensionless quantities of time t, space ( x, y, z), pressure P, depth d, water elevation η and velocity field ( U, V, W), via the following scaling relation

x = x L , P = P P , ϕ = ϕ h 1 , y = y h 1 , u = u U , d = d H 2 , z = z H 2 , v = v V , η = η H 2 , t = t T , w = w W . (16.1)
Finally, we define the non-dimensional Froude's number

F r = U gH 2 .
For the sake of clarity and simplicity dropping •, dividing the Euler equations (15.5) by ρ 0 , using the dimensionless variables (16.1), and reordering the terms with respect to the powers of µ 1 and µ 2 , the dimensionless incompressible Euler system (15.5) reads as follows:

∂ x u + ∂ y v + ∂ z w = 0 ∂ t u + u∂ x u + v∂ y u + w∂ z u + ∂ x P = 0 µ 1 ∂ t v + u∂ x v + v∂ y v + w∂ z v + ∂ y P = 0 µ 2 ∂ t w + u∂ x w + v∂ y w + w∂ z w + ∂ z P = - 1 F r 2
The boundary conditions (15.7), (15.8), (15.9) and (15.10) read ∂ t η + u∂ x η + v∂ y η = w and P = 0 on z = η u∂ x d + v∂ y d = w on z = d or u∂ x ϕ + w∂ z ϕ = v on y = ϕ and the irrotational equation (15.6) becomes

∂ y u = µ 1 ∂ x v, µ 1 ∂ z v = µ 2 ∂ y w, ∂ z u = µ 2 ∂ x w .
(16.9)

3d-2d model reduction and asymptotic expansions

In this first step, we focus on the width-averaging of System (16.2)-(16.5). In particular, we compute the asymptotic expansions of the velocity (u, v) and the pressure P of System (16.2)-(16.5) as a function of their width-averages. We obtain what we call "flat free surface approximation" property which means that the variations of the free surface following y can be neglected in the three-dimensional Euler system (16.2)-(16.5). Finally, the asymptotic approximation of the width-averaged Euler system is obtained. This averaged model is the starting point for the second step (see Section 16.3).

Asymptotic expansions of the fluid velocity

Given a function (t, x, y, z) → X(t, x, y, z) and (x, z) → q(x, z), we define X q as follows X q (t, x, z) := X (t, x, q(x, z), z)

and we use the notations

w = u w , div x,z [w] = ∂ x u + ∂ z w and ∇ x,z (u) = ∂ x u ∂ z u
to represent the vector w, the two dimensional divergence operator and the gradient operator with respect to the variable x and z.

Integrating the two first equations of the irrotational equations (16.9) for s ∈ [α(x, z), y], we get

u(t, x, y, z) = u α (t, x, z) + µ 1 y α(x,z) ∂ x v
= u α (t, x, z) - µ 1 2 ∂ x div x,z w α (t, x, z)(y -α(x, z)) 2 +O µ 2 1 µ 2 (16.14) and w(t, x, y, z) = w α (t, x, z) - µ 1 2µ 2 ∂ z div x,z w α (t, x, z)(y -α(x, z)) 2 +O µ 2 1 µ 2 2 .
(16.15)

Width-averaged Euler equations

Given a function (t, x, y, z) → X(t, x, y, z), we define the width-average of X by the quantity

X (t, x, z) := 1 σ(x, z) β(x,z) α(x,z)
X(t, x, y, z) dy where σ(x, z) = β(x, z)α(x, z) is the width of the section at the elevation z.

We average the equations (16.2)-(16.5) for y ∈ [α(x, z), β(x, z)] using Leibniz integral rule, to get the width-averaged Euler system:

                           ∂ ∂x (σ u ) + ∂ ∂z (σ w ) = 0 ∂ ∂t (σ u ) + ∂ ∂x (σ u 2 ) + ∂ ∂z (σ uw ) + ∂ ∂x (σ P ) = P β ∂β ∂x -P α ∂α ∂x µ 2 ∂ ∂t (σ w ) + ∂ ∂x (σ uw ) + ∂ ∂z (σ w 2 ) + ∂ ∂z (σ P ) = - σ F r
(x, z) u (t, x, z) = σ(x, z)u α (t, x, z) - µ 1 6 ∂ x div x,z w α (t, x, z)σ(x, z) 3 +O µ 2 1 µ 2 , σ(x, z) w (t, x, z) = σ(x, z)w α (t, x, z) - µ 1 6µ 2 ∂ z div x,z w α (t, x, z)σ(x, z) 3 +O µ 2 1 µ 2 2 , σ(x, z) u 2 (t, x, z) = σ(x, z)u 2 α (t, x, z) + O µ 2 1 µ 2 , σ(x, z) w 2 (t, x, z) = σ(x, z)w 2 α (t, x, z) + O µ 2 1 µ 2 2 , σ(x, z) uw (t, x, z) = σ(x, z)u α (t, x, z)w α (t, x, z) -u α (t, x, z) µ 1 6µ 2 ∂ z div x,z w α (t, x, z)σ(x, z) 3 +O µ 2 1 µ 2 2 .

irrotationality:

We lose the irrotational condition (16.9) by width-averaging, since we get

∂ ∂z (σ u ) = µ 2 ∂ ∂x (σ w ) + (u β ∂ z β -w β ∂ x β -u α ∂ z α + w α ∂ x α) .
However, using the last equation in (16.9) together with Equations (16.14) and (16.15), we obtain

∂u α ∂z = µ 2 ∂w α ∂x + O(µ 1 ) . ( 16.22)
The approximated irrotational condition can be also written as a function of ( u , w ) by means of Equations (16.17) and (16.18).

asymptotic expansion of the pressure:

The last equation in System (16.16) allows to write

P β = P α + O(µ 1 ) . ( 16.23)
Thanks to (16.23), on one hand, the terms in the right hand side of the second and third equations of System (16.16) can be simplified in

P β ∇ x,z (β) -P α ∇ x,z (α) = P α ∇ x,z (
y -(t,x) η(t, x, y) -d(x, y) dy = y + (t,x) y -(t,x) η eq (t, x) -d(x, y) dy + O(µ 1 ) = η eq (t,x) d * (x) σ(x, z) dz + O(µ 1 ) = Ω eq (t,x) dy dz + O(µ 1 ) = A eq (t, x) + O(µ 1 ) (16.28)
Therefore, in the following, we consider the fluid cross-section Ω eq instead of Ω for which computations are easier. In the following, we set η eq as η eq (t, x) = η(t, x, y * (x)) .

width-averaged equations:

Using all the previous equations (16.17

                                 div x,z [σw α ] + O µ 2 1 µ 2 2 = µ 1 6µ 2 ∂ ∂z σ∂ z div x,z w α σ 3 ∂ ∂t (σu α ) + div x,z [σu α w α ] + ∂ ∂x (σP α ) + O µ 2 1 µ 2 2 = P α ∂σ ∂x + µ 1 6µ 2 ∂ x u α ∂ z div x,z w α σ 3 µ 2 ∂ ∂t (σw α ) + div x,z [σw α w α ] + ∂ ∂z (σP α ) = - σ F r 2 +P α ∂σ ∂z + O(µ 1 ) (
Ω (t, x) = z ∈ R; d * (x) ≤ z ≤ η eq (t, x) .
For a given time t ∈ [0, T] and x ∈ [0, L c ], the boundary of the fluid cross-section ∂ Ω (t, x) is composed of a free surface Γ fs (t, x) = {z ∈ R; z = η eq (t, x)} and a wet boundary

Γ wb (x) = {z ∈ R; z = d * (x)} such that ∂ Ω (t, x) = Γ fs (t, x) ∪ Γ wb (x)
The unit outward normal vector to the free surface is

n fs (t, x) = 1 1 + ∂ x η eq (t, x) 2 -∂ x η eq (t, x) 1
and the unit outward normal vector to the wet boundary is σ(x, s) ds the wet area of water between d * (x) and z.

n wb (x) = 1 1 + (∂ x d * (x)) 2 ∂ x d * (x) -1 .
In the sequel, we consider the fluid cross-section (16.26) Ω eq (t, x) instead of Ω(t, x) thanks to (16.28) (see also (15.3)). Thus, we define the local height of the water in Ω eq (t, x) by h eq (t, x, y) = η eq (t, x)d(x, y) .

In this new setting, the section-average feq of a function f can be defined by feq = 1 A eq (t, x)

η eq (t,x) d * (x) β(x,z) α(x,z) f (t, x, y, z) dydz .

Asymptotic expansion of the fluid velocity

Let us first integrate Equation (16.22) .

for s ∈ [d * (x), z], to get u α (t, x, z) = u b (t, x) + µ 2 z d * (x) ∂w α ∂x ds + O(µ 1 ) . ( 16
By means of Equation (16.32), the previous equation becomes:

w α (t, x, z) = - 1 σ(x, z) ∂ ∂x (u b (t, x)S(x, z)) + O(µ 2 ) .
By injecting Equation (16.33) in (16.32), we get

u α (t, x, z) = u b (t, x) -µ 2 z d * (x) ∂ x S(u b , x, s) ds + O(µ 2 2 ) (16.34) where S(u, x, z) = 1 σ(x, z) ∂ ∂x (uS(x, z)) . (16.35)
As a consequence, we have ūeq = 1 A eq (t, x)

η eq (t,x) d * (x) β(x,z) α(x,z) u(t, x, y, z) dy dz = 1 A eq (t, x) η eq (t,x) d * (x) σ(x, z) u (t, x, z)dz = 1 A eq (t, x) η eq (t,x) d * (x) σ(x, z)u α (t, x, z) dz + O(µ 1 )
and through Equation (16.34), we obtain

u b (t, x) = ūeq (t, x) + µ 2 A eq (t, x) η eq (t,x) d * (x) σ(x, z) z d * (x) ∂ x S(u b , x, s) ds dz + O(µ 2
2 ) .

(16.36)

Using the fact that

u b (t, x) = ūeq (t, x) + O(µ 2 )
in Equation (16.36) allows to write u b (t, x) as a function of ūeq :

u b (t, x) = ūeq (t, x) + µ 2 A eq (t, x) η eq (t,x) d * (x) σ(x, z) z d * (x) ∂ x S(
= u α (t, x, z) + O(µ 1 ) = u b (t, x) -µ 2 z d * (x) ∂ x S(u b , x, s) ds + O(µ 2 2 ) = ūeq (t, x) + µ 2 B 0 ( ūeq , x, z) + O(µ 2 2 ) (16.38)
where where

B 0 ( ūeq , x, z) = 1 A eq (t, x) η eq (t,x) d * (x) σ(x, z) z d * (x) ∂ x S( ūeq , x, s) ds dz - z d * (x) ∂ x S(
D(u) = (∂ x u) 2 -∂ t ∂ x u -u∂ xx u , (16.43)
and

Q(u, S, σ) = u 2 σ(x,s) ∂ x S(x,s)∂ x σ(x,s) σ(x,s) -∂ xx S(x, s) +∂ x u 2 2 S(x,s)∂ x σ(x,s) σ(x,s) 2 -(∂ t u + u∂ x u) ∂ x S(x,s) σ(x,s) . (16.44)
The non-hydrostatic pressure (16.42) is more suitable to have a formulation akin to the SGN equations (14.2)-(14.4).

the incompressibility condition, Φ satisfies the following indicator transport equation

∂ t Φ + ∂ x (Φu) + div y,z [Φv] = 0 on Ω eq (t) (17.2)
where div y,z [v] = ∂ y v + ∂ z w with v = (v, w).

equation of the conservation of the mass

Integrating Equation (17.2) over the section Ω eq (t, x) and using Leibniz integral rule, we get the following conservation of the mass equation

Ω eq (t,x) ∂ t Φ + ∂ x (Φu) + div y,z [Φv] dy dz = ∂ t A eq + ∂ x Q eq = 0 (17.3)
where Q eq is the water discharge

Q eq = A eq ūeq (17.4)
and ūeq is Ω eq (t, x) section-averaged velocity given by ūeq = 1 A eq (t, x) Ω eq (t,x) u(t, x, y, z) dydz .

equation of the conservation of the momentum

In order to get the momentum equation of the section-averaged free surface model, we integrate each terms of (16.3) along the section Ω eq (t, x) as follows:

Ω eq (t,x) ∂ t (u) a 1 + ∂ x (u 2 ) a 2 + div y,z [uv] a 3 + ∂ x P a 4
dy dz = 0 .

2 ), the non-linear term u 2 can be written as:

u 2 (t, x, y, z) = ūeq (t, x) 2 + 2µ 2 ūeq B 0 ( ūeq , x, z) + O µ 2 2 .
Therefore, we have

Ω eq (t,x) ∂ x (u 2 ) dy dz = ∂ x Ω eq (t,x) u 2 dy dz - ∂Ω eq (t,x) u 2 ∂ x m • n ds i. e. Ω eq (t,x) ∂ x (u 2 ) dy dz = ∂ x (Q eq 2 /A eq ) + µ 2 ∂ x Q eq B(x, Q eq /A eq ) A eq - ∂Ω eq (t,x) u 2 ∂ x m • n ds + O µ 2 2 (17.6)
where B stands for the quadratic part of the fluid velocity, defined via B 0 (16.39), and B(x, Q eq /A eq ) = Ω eq (t,x) 2B 0 (Q eq /A eq , x, z) dy dz

= 2 η eq (t,x) d * (x)
σ(x, z)B 0 (Q eq /A eq , x, z) dz = 0 .

Ω eq (t,x) a 1 + a 2 + a 3 dy dz = ∂ t Q eq + ∂ x (Q eq 2 /A eq ) +µ 2 ∂ x Q eq B(x, Q eq /A eq ) A eq - ∂Ω eq (t,x) u (∂ t m + u∂ x m -v) • n ds +O µ 2 2 = ∂ t Q eq + ∂ x (Q eq 2 /A eq ) + O µ 2 2 . ( 17
∈ [d * (x), η eq (t, x)] → m(x, z) = ϕ(x, z) z
where ϕ is given by (15.1). The unit speed curve parameterisation is

∂ y m(x, z) = 1 + ∂ y ϕ(x, z) 2 .
The outward unit normal is given by

n =            1 ||n|| -1 ∂ z α if ϕ(x, z) < y * (x) 1 ||n|| 1 ∂ z β otherwise where ||n(x, z)|| = ∂ y m(x, z) .
Thus, we obtain

Γ wb (x) P∂ x m • n ds = η eq (t,x) d * (x) P(t, x, z)∂ x σ(x, z) dz .
hydrostatic part: (see equation (16.41)) Using P = P h in Equation (17.9), we get:

Ω eq (t,x) Each term DI j for j = 1, 2 is the non-hydrostatic counterpart of the term I j for j = 1, 2.

Remark 17.2.1. The terms I j (t, x) and DI j (t, x), respectively, can be written I j (x, A eq (t, x)) and DI j (x, A eq (t, x), Q eq (t, x) for j = 1 and j = 2. In what follows, we make use of these notations.

Finally, gathering Equation (17.10) and Equation (17.13), noting G(x, A eq , Q eq ) = -∂ x η eq d * (x) σ(x, z)

η eq z Q Q eq
A eq , S, σ ds dz +DI 2 , (17.16) we obtain Ω eq (t,x) ∂ x P dy dz = ∂ x (I 1 (x, A eq )) -I 2 (x, A eq ) +µ 2 ∂ x (G(x, A eq )D( ūeq )) -µ 2 G(x, A eq , Q eq ) (17.17)

Gathering results (17.8) and (17.17), we get the equation of the conservation of the momentum.

the dispersive model for arbitrary non rectangular channel/river

From now on, we omit the notations X eq and X for the sake of readability.

Gathering Equations (17.3), (17.8) and (17.17), we present the new one-dimensional dispersive model for open channel/river flows

         ∂ t A + ∂ x Q = 0 ∂ t Q + ∂ x Q 2 A + I 1 (x, A) + µ 2 ∂ x (G(x, A)D(u)) = I 2 (x, A) +µ 2 G(x, A, Q) + O(µ 2 2 ) (17.18)
where A is the wet area (16.27), Q is the water discharge (17.4), I 1 (resp. I 2 ) is the hydrostatic pressure (resp. source) terms (17.11) (resp. (17.12)), G(x, A) (17.14) generalises h 3 3 in the classical SGN equations (see (14.3)), G(x, A, Q) (17.16) extends the uneven bottom source term in the classical SGN equations and ∂ x (G(x, A)D(u)) is the dispersive term where D is given by (16.43).

The new section-averaged model extends the section-averaged free surface model for open channel flows [START_REF] Bourdarias | A mathematical model for unsteady mixed flows in closed water pipes[END_REF][START_REF] Ersoy | Modeling, mathematical and numerical analysis of various compressible or incompressible flows in thin layer[END_REF] by taking µ 2 = 0 and extends the SGN equations for uneven bottom [START_REF] Barthélemy | Nonlinear shallow water theories for coastal waves[END_REF][START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF][START_REF] Seabra-Santos | Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle[END_REF][START_REF] Serre | Contribution à l'étude des écoulements permanents et variables dans les canaux[END_REF] to arbitrary channel/river section. We consider the motion of an incompressible and irrotational fluid with constant density ρ 0 > 0 in a three dimensional domain (see Figure 18 where h(t, x) is width-averaged of the local height of the water h(t, x, y) := η(t, x, y)d(x), i. e. , h(t, x) = η(t, x)d(x).

The boundary of the domain Ω(t) is defined by ∂Ω(t) and is decomposed into four parts: the free surface Γ fs (t), the wet boundary Γ wb (t), the inflow boundary Γ i (t) and the ouflow boundary Γ o (t). The wet boundary can be decomposed itself in three parts: the bottom Γ b (t), the left lateral boundary Γ lb (t) and the right one Γ rb (t).

Thus, in contrast with the non-rectangular case, the hydrostatic (17.12) and the non-hydrostatic (17.15) Thus, the model (17.18) is the same as before but with the formula (18.1) and (18.2) instead of (17.12) and (17.15).

reformulated model for rectangular case

The main drawback of Equations (17.18) is that it has third order terms in space which may lead to instabilities at the numerical level. Therefore, we first propose a more stable formulation of Equations (17.18) before to present its numerical approximation. In particular, one can explicitly compute those operators:

• if σ ∈ R + * and d ∈ R are constant then we recover the standard one dimensional SGN equations (see for instance [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF][START_REF] Lannes | Nonlinear wave-current interactions in shallow water[END_REF] • if σ ∈ R + * is constant and d = d(x) then we recover the standard one dimensional SGN equations (see for instance [START_REF] Lannes | Derivation of asymptotic two-dimensional time-dependent equations for surface water wave propagation[END_REF][START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF][START_REF] Lannes | Nonlinear wave-current interactions in shallow water[END_REF]) over uneven bottom with -u 2 σ (x)

L[A, d, σ](u) = L 1 [A, d, σ](u) = L 0 [A, σ](u) -1 σh ∂ x
σ 2 h 3 3 + u 2 d (x) h 2 2 +∂ x u 2 2 σ (x) σ h 3 3 -u 2 d (x) σ (x) σ h 2 2 -u 2 σ (x) (d (x)) 2 h + u 2 σ (x)d (x) h 2 2 -∂ x u 2 2 σ (x)d (x) h 2 2 .
It is known that third order derivatives involved in the initial model (17.18) may create high frequencies instabilities, but the presence of the I dµ 2 L[A, d, σ] -1 in the second equation of (18.3) stabilises the equations with respect to these perturbations. Therefore, in the following, we construct a numerical scheme for Equations (18.3) instead of Equations (17.18).

a well-balanced finite volume approximation

This section is devoted to the numerical method to solve the reformulated dispersive model (18.3). It is rather natural to split the hyperbolic part to the dispersive one as done by several authors (see for instance [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF][START_REF] Chazel | Numerical simulation of strongly nonlinear and dispersive waves using a Green-Naghdi model[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: boundary conditions and validation[END_REF]).

Let N ∈ N * . Let us consider the following uniform mesh on [0, L c ]. Cells are denoted for every i ∈ [0, N + 1], by m i = (x i-1/2 , x i+1/2 ) with x i =

x i-1/2 +x i+1/2 2 the cell center and δx = x i+1/2x i-1/2 the space mesh. The interfaces x 1/2 = 0 and x = x N+1/2 denote the upstream and the downstream ends. We also consider a time discretisation t n defined by t n+1 = t n + δt n where the time step δt n is computed through a CFL condition related to the hyperbolic part.

Let us first highlight that the still water steady state for Equations (18.3) is independent of µ 2 . Indeed, one has ∀µ 2 > 0, the still water steady state equation reads u = 0, A σ + d = h 0 for some positive h 0 . As a consequence, the construction of a wellbalanced scheme can be easily achieved considering only the hyperbolic part of Equations ( 18.3), for instance, by the use of the hydrostatic reconstruction (see for instance [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]). where U i = (A i , A i u i ) T ≈ 1 δx m i (A, Au) T dx. Let us introduce the flux Then, one can write System (18.3) as follows:

F 1 (U) = Q, F 2 (U) = Q 2 /A ,
∂ t A + ∂ x F 1 (U) = 0 I d -µ 2 L[A, d, σ] (∂ t Q + ∂ x F 2 (U)) + ∂ x F 3 (x, U) -S(x, U) = 0
With these settings, we define the following numerical scheme: such that whenever µ 2 = 0, we recover the classical numerical scheme1 for the hyperbolic part

A n+1 i = A n i -δt n δx F 1 U -,n i+1/2 , U +,n i+1/2 -F 1 U -,n i-1/2 , U +,n i-1/2 Q * i = Q n i -δt n δx F 2 U -,n i+1/2 , U +,n i+1/2 -F 2 U -,n i-1/2 , U +,n i-1/2 Q n+1 i = Q * i -δt n δx (Y n ) i where A n Y n = F 3 x i+1/2 , U -,n i+1/2 , U +,n
U n+1 i = U n i -δt n δx
× F (x i+1/2 , U -,n i+1/2 , U +,n i+1/2 ) -F (x i-1/2 , U -,n i-1/2 , U +,n i-1/2 ) with F (x, U, V) = (F 1 (U, V), F 2 (U, V) + F 3 (x, U, V)). In these expressions,

s i+1/2 = max j=1,2 λ j (x i+1/2 , U -,n i+1/2 ) , λ j (x i+1/2 , U +,n i+1/2 )
where λ j (x, U) = Q/A + (-1) j A σ(x)F 2 r , j = 1, 2 are the eigenvalues of the Jacobian matrix of (F 1 , F 2 + F 3 ) T .

The numerical scheme is consistent and stable under the CFL condition max 1≤i≤N

(|λ 1 (x i , U n i )| , |λ 2 (x i , U n i )|)
δt n δx ≤ 1 .

numerical test case

In this section, we test the accuracy of the method and we show numerically the influence of the section variation in the case of the propagation of a solitary wave. For this purpose, we consider the exact solitary wave solutions of the Green-Naghdi equations in the one-dimensional setting over a flat bottom (see [START_REF] Lannes | A new class of fully nonlinear and weakly dispersive Green-Naghdi models for efficient 2D simulations[END_REF]), given in variables with dimensions, by η(t, x) = asech 2 (k(xct)), u(t, x) = c η(t,x) η(t,x)+z 0 with k = √ 3a 2z 0 √ z 0 +a and c = g(z 0 + a) (18.4) where z 0 is the depth of the fluid and a is the relative amplitude.

accuracy: The propagation of the solitary wave (18.4) is initially centered at x 0 = 10 m with a relative amplitude a = 0.2 m over a constant water depth z 0 = 1 m. The computational domain is L c = 100 m and it is discretized with N cells. The single solitary wave propagates from left to right. In this test, since the solitary wave is initially far from boundaries, the boundary conditions do not affect the computation, thus we choose to impose free boundary conditions at the downstream and upstream ends. The exact solution is computed in a channel of width σ = 1.

In what follows, we quantify the numerical accuracy of our numerical scheme by computing the numerical solution for this particular test case for an increasing number of cells N over a duration T = 20 s. Starting with N = 100 number of cells, we successively multiply the number of cells by two. For all n, we compare, in Figure 18.2, M n := max 0≤i≤N+2 (h n i ) of our numerical solution provided by Equations (18.3) with the exact one M(t n ) := max h(t n , x) x∈[0,L c ] = 2.2 given by (18.4). One can easily remark that the first order discretisation is not accurate for long time simulation due to the numerical dissipation. However, to limit the numerical dissipation of the first order numerical scheme, one can either limit the simulation time or consider a very large number of cells. However, it is better to increase the order of the numerical scheme but this is left to a future work. Therefore, in what follows, we consider a shorter simulation time and a large number of cells, just to illustrate the influence of the variation of the channel. influence of the section variation: We consider again the propagation of a solitary wave initially centered at x 0 = 10 m of relative amplitude a = 0.2 m, over a constant water depth z 0 = 1 m onto a computational domain of L c = 50 m and discretized with N = 5000 cells. Initially starting with (η(0, x), u(0, x)) (see Equations (18.4)), we compute the numerical simulation for the channels defined by σ(x; ε) = β(x; ε)α(x; ε)

with β = 1 2 - ε 2 exp -ε 2 x -L/2) 2
and α = -β with = 0, = 0.1, = 0.2, = 0.3 and = 0.4. The obtained results are presented in Figure 18.3. In Figure 18.3a, for each geometry, we show the evolution of the maximum of the water level M n := max 0≤i≤N+2 (h n i ). As expected, since the first part for x ≤ 25 is linearly converging, the water level increases while for x > 25, the channel is linearly diverging and therefore, the amplitude of the water level decreases. Moreover, in all numerical simulations, the mass is conserved. Indeed, for each value of ε, we have displayed in Figure 18.3b, the ratio of m n m 0 where m n := 1 N+2 ∑ N+1 i=0 A n i is the mass of water at time t n . The ratio m n m 0 is almost equal to 1, up to the order of accuracy of the numerical scheme. In what follows, we quantify the numerical accuracy of our numerical scheme. Starting with N = 100 number of cells, we successively multiply the number of cells by two. The errors on the water surface deformation are presented in Table 18.1 and in Table 18.2. These errors are computed at t = 8 s using the L 2 :

η num -η ref 2 = δx ∑ i η num i (t = 8) -η ref (t = 8, x i )
and the L ∞ norms where η ref is the exact solution in the case ε = 0 and is a reference one computed with 10 000 cells for ε = 0.4. Since the results are almost the same whatever ε is, we have decided to present only the results for ε = 0 and ε = 0.4 in Table 18.1 and in Table 18.2.

As expected, the obtain numerical order is slow because of the numerical dissipation of the solitary wave (as already pointed out in several works, see for instance [START_REF] Bourdarias | A numerical scheme for an improved Green-Naghdi model in the Camassa-Holm regime for the propagation of internal waves[END_REF] for which we obtain almost the same order of convergence in the case of uniform section). Moreover, as expected, cross-sectional variations have no influence on the convergence rate. Let us just emphasise that the convergence rates are slightly better in the case of non-uniform section because we are comparing our results to a reference solution and not to the exact one.

N

||η In this part, we have presented the first non-linear and weakly dispersive model for natural rivers and channel hydraulics. The obtained model is a generalisation of the well-known one dimensional Serre-Green-Naghdi equation.

In the specific case of rectangular channels, we have proposed a stable reformulation of the model for numerical purpose. To this end, we have constructed a first order well-balanced finite volume numerical scheme for which we have proposed some simple computations to study the influence of the section variation. Our numerical scheme is stable under the classical CFL condition of the hyperbolic part of the equations. However, it is clear that, unless to use fine grids, the numerical scheme is too dissipative to capture correctly the dispersive waves.

On one hand, the natural extension is to increase the order of this numerical scheme, not only for the rectangular sections but for arbitrary ones. The accuracy can be improved by using Finite Volume adaptive methods as done in Part iii or Continuous Galerkin methods1 , or even better Discontinuous Galerkin methods (see Part iv), especially in the context of Dissipative SGN models, see for instance [START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Favrie | A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves[END_REF].

On the other hand, at least for the performed numerical test case, it seems, that the operator is invertible. Theoretical mathematical and numerical framework are addressed in a forthcoming works.

Figure 7 . 7 Figure 7 . 8

 7778 Figure 2.1état à surface libre p(x, A, 0) (gauche), état en charge avec surpression p(x, A, 1) > 0 (milieu), état en charge avec dépression p(x, A, 1) < 0 (droite) où S section d'eau physique, S section de la conduite, R rayon de la conduite. . . . . . Figure2.2 Cas tests numériques. . . . . . . . . . . . . . . . . Figure 3.1 Illustration numérique de l'analyse de sensibilité de (v ε , δϕ ε ) par rapport à g pour le problème avec f (u) = u 2 2 . . . . . . . . . . . . . . . . . . . . . Figure 3.2 Exemple 1D et 2D. . . . . . . . . . . . . . . . . . . Figure 3.3 Exemple 3D. . . . . . . . . . . . . . . . . . . . . . Figure 3.4 Exemples de résultats de simulation (1). . . . . . Figure 3.5 Exemples de résultats de simulation (2). . . . . . Figure 3.6 Cas test 1d d'infiltration. . . . . . . . . . . . . . . Figure 5.1 Numerical illustration of the sensitivity of (v ε , δϕ ε ). 46 Figure 6.1 Multiple steady shocks (dashed line: source term, solid line: numerical solution). . . . . . . . . . . . Figure 6.2 A steady solution emerging in finite time with g c (dashed line: source term, solid line: numerical solution). . . . . . . . . . . . . . . . . . . . . . Figure 6.3 A local finite time emerging steady solution with g dc (dashed line: source term, solid line: numerical solution). . . . . . . . . . . . . . . . . . Figure 7.1 Example of hierarchical dyadic tree. . . . . . . . Figure 7.2 Mesh refinement and mesh coarsening: construction of the fluxes. . . . . . . . . . . . . . . . Figure 7.3 Illustration of a global time advancement with a first order scheme at time t n + δt n . . . . . . . . Figure 7.4 Illustration of a global time advancement with a first order scheme at time t n + 2δt n . . . . . . . . Figure 7.5 Illustration of a global time advancement with a first order scheme at time t n + 3δt n . . . . . . . . Figure 7.6 Illustration of a global time advancement with a first order scheme at time t n + 4δt n . . . . . . . . Figure 7.7 Projection's methods. . . . . . . . . . . . . . . . .

  . . . . . . . . . . . . . . . . . . . A C R O N Y M S Part I I N T R O D U C T I O N L'alchimie est à mi-chemin entre la poésie et les mathématiques, entre le monde du symbole et celui du nombre.

3 4 5

 345 Fluide compressible, Fluide incompressible, Couplage compressible-incompressible, Modèle hydrostatique, Modèle non-hydrostatique, Modèle dispersif, Modèle à faible Mach, Existence et stabilité de solutions faibles, Analyse asymptotique, Contrôle, Problème de Riemann, Volumes Finis, Schéma cinétique, Schéma équilibre, Raffinement de maillage adaptatif h et hp. Le mémoire de thèse "Modélisation, analyse mathématique et numérique de divers écoulements compressibles ou incompressibles en couche mince" est composé de 253 pages réparties en une préface, de trois parties totalisant cinq chapitres, d'une annexe et d'un index. Il s'articule autour du plan suivant : partie i Modélisation des équations en conduite fermée chapitre 1 Les équations Pressurised and Free Surface (PFS) chapitre 2 Approximation des équations PFS par Volumes Finis partie ii Autour des équations Primitives Compressibles chapitre Existence et stabilités de solutions faibles partie iii Travaux en cours et perspectives chapitre Vers la cavitation chapitre Modélisation de la sédimentation

  d'eau mixtes en conduites fermées.

Figure 2 .

 2 1).

Figure 2 . 1 :

 21 Figure2.1: état à surface libre p(x, A, 0) (gauche), état en charge avec surpression p(x, A, 1) > 0 (milieu), état en charge avec dépression p(x, A, 1) < 0 (droite) où S section d'eau physique, S section de la conduite, R rayon de la conduite.

  Comparaison de notre modèle (CSV, schéma cinétique) avec le modèle uniforme à section équivalente (CE). conduite (b) Cas test de Wiggert.

  Cas test de Baines (schéma cinétique).

Figure 2 . 2 :

 22 Figure 2.2: Cas tests numériques.

Figure 3 . 1 :

 31 Figure 3.1: Illustration numérique de l'analyse de sensibilité de (v ε , δϕ ε ) par rapport à g pour le problème avec f (u) = u 2 2 .

  alors les deux cellules sont fusionnées pour former une cellule C k b .

Figure 3 . 2 :

 32 Figure 3.2: Exemple 1D et 2D.

Figure 3 . 3 :

 33 Figure 3.3: Exemple 3D.

Figure 3 . 4 :

 34 Figure 3.4: Exemples de résultats de simulation (1).

  (a) Rupture de barrage 2D (Eqs Euler bi-fluide) : (a) Maillage; (b) Densité (air-bleu, eau-rouge); (c) Densité de production numérique d'entropie (valeurs vert-zero, bleu-negatif); (d) Niveau de raffinement par bloc (1 à 5); (e) Expérience Koshizuka (95); (f) Critère de raffinement par bloc. (b) Tsunami (vallée de la Monai) : temps t = 16s (Eqs Saint-Venant). Les couleurs correspondent aux valeurs de la densité de production numérique.

Figure 3 . 5 :

 35 Figure 3.5: Exemples de résultats de simulation (2).

3. 3

 3 dynamique des écoulements dans les plages sableuses 3.3 dynamique des écoulements dans les plages sableuses motivations.

  où ψ indique la pression (m), K le tenseur de conductivité hydraulique (m/s), z l'élévation (m) et θ la teneur en eau. Ce système est fermé par deux lois constitutives sur la teneur en eau et la conductivité hydraulique. La méthode Galerkin discontinue (ou Discontinuous Galerkin (DG) method) est utilisée car elle offre des avantages appropriés pour l'adaptation hp. Les méthodes discontinues de Galerkin peuvent être considérées comme une méthode partageant des propriétés à la fois des éléments finis et des volumes finis : elles sont basées sur une formulation variationnelle mais de facon élémentaire. Elles sont donc localement conservatrices, ce qui est crucial en dynamique des fluides ([35-39] et Figure 3.6).

  (a) Profile de charge. (b) Raffinement de maillage adaptatif autour du front.

Figure 3 . 6 :

 36 Figure 3.6: Cas test 1d d'infiltration.

Corollary 4 . 1 . 1 .

 411 Let the flux f be a locally Lipschitz function on R. Let v ∈ L ∞ (R) be an entropy solution to the stationary problem (4.1) and u an entropy solution to the evolutionary problem (4.2), with u(0

Figure 5 . 1 :

 51 Figure 5.1: Numerical illustration of the sensitivity of (v ε , δϕ ε ).

Figure 6 . 1 :

 61 Figure 6.1: Multiple steady shocks (dashed line: source term, solid line: numerical solution).

Figure 6 . 2 :

 62 Figure 6.2: A steady solution emerging in finite time with g c (dashed line: source term, solid line: numerical solution).

Figure 6 . 3 :

 63 Figure 6.3: A local finite time emerging steady solution with g dc (dashed line: source term, solid line: numerical solution).

Figure 7 . 1 :

 71 Figure 7.1: Example of hierarchical dyadic tree.

Figure 7 . 2 :

 72 Figure 7.2: Mesh refinement and mesh coarsening: construction of the fluxes.

. 3 ,

 3 Figure 7.4, Figure 7.5 and Figure 7.6.

Figure

  Figure 7.3: Illustration of a global time advancement with a first order scheme at time t n + δt n .

  (a) Osher and Sanders projection. (b) Tang and Warnecke projection.

Figure 7 . 7 :

 77 Figure 7.7: Projection's methods.

2 k

 2 but the global conservativity is lost.

Figure 7 . 8 :

 78 Figure 7.8: Example of two dimensional BB-AMR with 3 domains (blue-yellowgreen) and 9 blocks.

Figure 7 . 9 :

 79 Figure 7.9: Example of three dimensional block-based mesh with 3 domains and 27 blocks.

Figure 8 . 1 :

 81 Figure 8.1: The distribution function d, the first order d and the second order d derivatives for the mesh refinement criterion S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25, x B = 3.75, a = 2 and b = 10.

  f (α) with a = 200 and b = 1000 f (α) with with a = 2 and b = 10

Figure 8 . 2 :

 82 Figure 8.2: The function f for the mesh refinement criterion S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25 and x B = 3.75.

Figure 8 . 3 :

 83 Figure 8.3: The mesh refinement criterion S(x) = a exp(-b(xx B ) 2 ) + exp(-5(xx A ) 2 ) where x A = 1.25 and x B = 3.75.

2 Figure 8 . 4 : 7 Figure 8 . 5 :

 284785 Figure 8.4: Numerical results for the water height at time t = 2 s using the automatic threshold α PE . N m is the mean number of cells used during the simulation.

Figure 8 . 7 -

 87 Figure 8.10, the propagation of the wave at time t = 11.25 s (Figure 8.7), t = 13.25 s sources and pictures are available at http://isec.nacse.org/workshop/2004 _ cornell/bmark2.html

Figure 8 . 6 :

 86 Figure 8.6: Experimental settings.

Figure 8 . 7 :

 87 Figure 8.7: Numerical result for the water height at time t = 11.25s. The colours correspond to the magnitude of the kinetic energy.

Figure 8 . 9

 89 and Figure 8.10 show respectively the large flooding and run-up of the wave and its reflection on the cliff. The numerical solutions are now quantitatively compared to the experimental measures. The comparisons are shown in Figure 8.11 at the three waves gauges WG 1 (see Figure 8.11a), 2 (see Figure 8.11b) and 3 (see Figure

Figure 8 . 8 :

 88 Figure 8.8: Numerical result for the water height at time t = 13.25s. The colours correspond to the magnitude of the kinetic energy.

Figure 8 . 9 :

 89 Figure 8.9: Numerical result for the water height at time t = 16s. The colours correspond to the magnitude of the kinetic energy.

Figure 8 . 10 :

 810 Figure 8.10: Numerical result for the water height at time t = 17.5s. The colours correspond to the magnitude of the kinetic energy.

Figure 8 .

 8 14 and Figure 8.15.

3 Figure 8 . 11 :

 3811 Figure 8.11: Free surface results at different positions : experimental data versus numerical simulation with and without mesh adaptivity.

Figure 8 . 12 :

 812 Figure 8.12: Adaptive mesh : time t = 11.25s. The colours correspond to the numerical density of entropy production S.

Figure 8 . 13 :

 813 Figure 8.13: Adaptive mesh : time t = 13.25s. The colours correspond to the numerical density of entropy production S.

Figure 8 . 14 :

 814 Figure 8.14: Adaptive mesh : time t = 16s. The colours correspond to the numerical density of entropy production S.

Figure 8 . 15 :

 815 Figure 8.15: Adaptive mesh : time t = 17.5s. The colours correspond to the numerical density of entropy production S.

Chapter 10 introduces

 10 Richards' equation to get the model problem. In Chapter 11, Richards' equation discretization is presented through a DG framework. Then, several numerical experiments are discussed in Chapter 12.

(11. 1 )

 1 with T the total duration and h 0 the initial condition. Ω is subdivided into elements E which form the partition denoted by E h . Elements are made of faces F. F I h stands for the set of interior faces, F D h for the set of Dirichlet boundary faces and F N h for the set of Neumann boundary faces.

Figure 11 . 1 :

 111 Figure 11.1: Refinement and coarsening process for one triangular and quadrangular element.

Figure 12 . 1 :

 121 Figure 12.1: Results for the one-dimensional infiltration numerical test case.

Figure 12 . 2 :

 122 Figure 12.2: Numerical simulation of two-dimensional drainage with watertable position compared to experimental data from[START_REF] Vauclin | Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils[END_REF].

Figure 12 . 3

 123 Figure 12.3 shows a schematic profile of the beach considered. Y-axis was multiplied by 20 to make it more visible. It's a straight slope of ratio 1:35 which spans over 105 meters. Hydraulic properties for constitutive laws of Equation (11.1) are written in the Table 12.1.

Figure 12 . 3 :

 123 Figure 12.3: Topography of the idealized beach and mesh used for the simulation.

Figure 12 . 4 :

 124 Figure 12.4: Snapshot of hydraulic head contours at t = 12.5s during the first wave propagation on calm waters.

Figure 12 . 5 :

 125 Figure 12.5: Snapshot of hydraulic head contours at t = 37.5 s for the end of the first uprush. Velocities range from about 0.3 to 0.6 cm/h in the beach. Maximum velocity is 35 cm/h for infiltration.

Figure 12 . 6 :

 126 Figure 12.6: Snapshot of hydraulic head contours at t = 57.5 s for the end of the first backwash. Most velocities are less than 1 cm/h.

Figure 12 . 7 :

 127 Figure 12.7: Snapshot of hydraulic head contours at t = 67.5 s for the next incoming uprush. Most velocities are less than 0.8 cm/h.

  E N S I O N R E D U C T I O N : F R O M H Y D R O S TAT I C T O N O N H Y D R O S TAT I C M O D E L SL'art de réduire n'est pas une dépréciation mais une volonté de compréhension.

T

  H E T H R E E -D I M E N S I O N A L I N C O M P R E S S I B L E E U L E R E Q U AT I O N S We start in Section 15.1 by reviewing the irrotational and incompressible Euler equations in the special geometric setting, describing the physics with a wet boundary on the bottom of the water course and a free surface on the top. Boundary conditions are presented in Section 15.2. 15.1 geometric set-up and the euler equations Let T > 0 be an arbitrary time. We consider an incompressible and irrotational fluid moving in the time-space box [0, T] × C with typical point (t, (x, y, z)) where C is the geometrical definition of a convex (non-rectangular) channel/river C = (x, y, z) ∈ R 3 ; x ∈ [0, L c ], α(x, z) ≤ y ≤ β(x, z) and d(x, y) ≤ z} . L c > 0 is the horizontal length of the domain and for z ≥ d(x, y), α(x, z) (respectively β(x, z)) is the left (respectively right) boundary point at the elevation z as displayed in Figure 15.1. The height of the surface of the water level and the boundary of the section are modelled, respectively, by the functions η(t, x, y) and d(x, y) with respect to a reference horizontal height z = 0. For all x ∈ [0, L c ], we assume that the function y → d(x, y), has a global minimum at y * (x) i. e. d * (x) = d(x, y * (x)) = min y d(x, y)

  [START_REF] Arnold | Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems[END_REF] and A(t, x) = |Ω(t, x)| is the wet area of the cross-section of fluid:A(t, x) = Ω(t,x) dω = y + (t,x) y -(t,x) η(t, x, y)d(x, y) dy = max y η(t,x,y) d * (x) β(x,z) α(x,z) 1 {d(x,y)≤z≤η(t,x,y)} dydz (15.4) where y -(t, x) := min{y ∈ R; η(t, x, y) = d(x, y)} and y + (t, x) := max{y ∈ R; η(t, x, y) = d(x, y)} stands for the left and the right transversal limit at z = η as shown in Figure 15.1.

Figure 15 . 1 :

 151 Figure 15.1: Geometric set-up.

  and Part iv, by considering the following boundary condition u • n wb = (∂ t m + I) • n wb where ∂ t m models the evolution in time of the bed and I is the infiltration function. I models the amount of water that leaves (I > 0) or enters (I < 0) the flow per elementary boundary element. W I D T H -AV E R A G E D A N D D E P T H -AV E R A G E D A S Y M P T O T I C E X PA N S I O N S

Figure 16 . 1 :

 161 Figure 16.1: Equivalent geometric set-up.

  )-(16.21) and (16.24)-(16.25) in System(16.16), we obtain the width-averaged Euler equations expressed as the couple of unknown (w α , P α ):

Figure 16 . 2 :

 162 Figure 16.2: Width-averaged geometric set-up.

T

  H E D I S P E R S I V E M O D E L F O R R E C TA N G U L A R S E C T I O NFor the specific case of the rectangular section, almost all the previous computations are the same. The changes are mainly in the geometrical definition and the boundaries of the channel.

  .[START_REF] Altazin | Numerical investigation of BB-AMR scheme using entropy production as refinement criterion[END_REF])Ω(t) = (x, y, z) ∈ R 3 ; x ∈ [0, L c ], α(x) ≤ yy * (x) ≤ β(x), d(x) ≤ z ≤ η(t, x)}where y * (x) describes the transversal variation of the channel with respect to the main channel direction, defined byy * (x) = α(x) + β(x) 2where α and β are the transversal limit of the channel. Here the bottom d is now a function of x only and d * (x) = d(x) by definition. Thus, in view of the definition of the fluid domain, since σ = σ(x) = β(x)α(x) and y -(x) = α(x) and y + (x) = β(x), the wet-area can be simply defined by A(t, x) = σ(x) h(t, x)

∂x u 2 2 S

 2 Skipping the technical details, defining a linear operator L (where L is defined below)L[A, d, σ](u) = AL[A, d, σ] u A ,one can show that System (17.18) can be written:t A + ∂ x (Au) = 0 I dµ 2 L[A, d, σ] ∂ t (Au) + ∂ x Au 2 +∂ x I 1 (x, A) + µ 2 AQ[A, d, σ](u) = I 2 (x, A) + O(µ 2 2 )(18.3)where I d stands for the identity operator, L is a linear operatorL[A, d, σ](u) = 1 A ∂ x T [A, d, σ] (u, σ) -T [A, d, σ] (u, ∂ x σ) + 1 A σ(x)d (x) T [A, d, σ, z = d(x)] (u)and Q is a quadratic operatorQ[A, d, σ](u) = 1 A ∂ x G[A, d, σ] (u, σ) -G[A, d, σ] (u, ∂ x σ) + 1 A σ(x)d (x) G[A, d, σ, z = d(x)] (u)with T , G are given byT [A, d, σ, z](u) = ∂ x (u) , d, σ, z](u) = η z 2 (∂ x u) 2 S(x,s) σ(x) + u 2 σ(x) ∂ x S(x,s)∂ x σ(x) σ(x) -∂ x ∂ x S(x, s) +∂ (x,s)∂ x σ(x) σ(x) 2ds with X [A, d, σ](u, ψ) = η d(x) ψX [A, d, σ, z](u) dz .

1 σh ∂ x 2 3 σh 3 (

 13 ) over flat bottom withL[A, d, σ](u) = L 0 [A, σ](u) , d, σ](u) = Q 0 [A, σ](u) = ∂ x u) 2 .

σh 2 2 1 σh ∂ x σ h 2 2 u 2 d 3 + u h 2 3 σu 2 σ h 3 3 + 1 3 3+∂ x ∂ x u 2 2 2 h 3 3

 21233313223 ud (x) + h 2 ∂ x ud (x)u (d (x)) 2 and Q[A, d, σ](u) = Q 1 [A, d](u) = Q 0 [A, σ](u) + (x) + h (∂ x u) 2 d (x) +u 2 d (x) d (x) .• if σ = σ(x) and d = d(x) then we get the generalised one dimensional SGN equations for non uniform rectangular channel over uneven bottom withL[A, d, σ](u) = L 1 [A, d, σ](u) + 1 σh ∂ x σ (x) h 3 3 u σ (x) σ ∂ x u h 2 (x) σu h 2 d (x) and Q[A, d, σ](u) = Q 1 [A, d, σ](u) + 1 σh ∂ x (σ (x)) 2 σh ∂ x d (x)σ (x)u 2 h 2 2 -1 σh ∂ x σ (x)u 2 h σ (x) h 3 3 -1 σh σ (x)R[A, d, σ](u) with R[A, d, σ](u) = (∂ x u) 2 h 3 3 + u 2 σ (x) σ + u 2 σ (x) σ d (x) h 2 2

F 3

 3 (x, U) = I 1 (x, A) + µ 2 G[A, d, σ] (u, σ) and S(x, U) = I 2 + µ 2 G[A, d, σ] (u, ∂ x σ) -µ 2 σ(x)d (x) G[A, d, σ, z = d(x)] (u) .

i+1/ 2 -

 2 F 3 x i-1/2 , U -,n i-1/2 , U +,n i-1/2 + µ 2 N n i 1≤i≤N.The matrix A n is the cell-centered approximation of the linear operatorI dµ 2 L[A, d, σ] and N n i is the cell-centered approximation of -G[A, d, σ] (u, ∂ x σ) + σ(x)d (x) G[A, d, σ, z = d(x)] (u).

Figure 18 . 2 :

 182 Figure 18.2: M n := max 0≤i≤N+2 (h n i ).

4 (Figure 18 . 3 :

 4183 Figure 18.3: Influence of σ.

  introduction 10.1 Richards equation and constitutive laws . . . . . . . . . . 10.2 Seepage boundary condition . . . . . . . . . . . . . . . . 11 dg method 11.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Nonlinear solver with adaptive time stepping . . . . . . 11.3 Adaptive mesh refinement . . . . . . . . . . . . . . . . . . 12 numerical results 12.1 One-dimensional infiltration numerical test case . . . . . 12.2 Comparison with water table drainage experiment . . . 12.3 Idealized beach test case . . . . . . . . . . . . . . . . . . . Geometric Equation of the conservation of the momentum . . . . . 17.3 The dispersive model for arbitrary non rectangular channel/river . . . . . . . . . . . . . . . . . . . . . . . . . . . . Reformulated model for rectangular case . . . . . . . . . 18.2 A well-balanced finite volume approximation . . . . . . . 18.3 Numerical test case . . . . . . . . . . . . . . . . . . . . . . Cas tests numériques. . . . . . . . . . . . . . . . . Figure 3.1 Illustration numérique de l'analyse de sensibilité de (v ε , δϕ ε ) par rapport à g pour le problème avec f (u) = u 2

	L I S T O F F I G U R E S
	Figure 2.1	état à surface libre p(x, A, 0) (gauche), état en
		charge avec surpression p(x, A, 1) > 0 (milieu),
		état en charge avec dépression p(x, A, 1) < 0
		(droite) où S section d'eau physique, S section
		de la conduite, R rayon de la conduite. . . . . .
	Figure 2.2
		13 concluding remarks and extensions
		v dimension reduction: from hydrostatic to non
		hydrostatic models
		14 introduction
		15 the three-dimensional incompressible euler equa-
		tions
		15.1 18 the dispersive model for rectangular section
		18.1 19 concluding remarks and extensions
		bibliography

set-up and the Euler equations . . . . . . . . . 15.2 Boundary conditions . . . . . . . . . . . . . . . . . . . . . 16 width-averaged and depth-averaged asymptotic expansions 16.1 Dimensionless Euler equations . . . . . . . . . . . . . . . 16.2 3d-2d model reduction and asymptotic expansions . . . 16.3 2d-1d like model reduction and asymptotic expansions . 17 a new non-linear dispersive model 17.1 Equation of the conservation of the mass . . . . . . . . . 17.2 2 . . . . . . . . . . . . . . . . . . . . . Figure 3.2 Exemple 1D et 2D. . . . . . . . . . . . . . . . . . . Figure 3.3 Exemple 3D. . . . . . . . . . . . . . . . . . . . . . Figure 3.4

Table 8 . 1

 81 Numerical result for the water height at time t = 13.25s. The colours correspond to the magnitude of the kinetic energy. . . . . . . . . . . . . . . . . Figure8.9 Numerical result for the water height at time t = 16s. The colours correspond to the magnitude of the kinetic energy. . . . . . . . . . . . . . . . . Figure 8.10 Numerical result for the water height at time t = 17.5s. The colours correspond to the magnitude of the kinetic energy. . . . . . . . . . . . . . . . . Figure 8.11 Free surface results at different positions : experimental data versus numerical simulation with and without mesh adaptivity. . . . . . . . . . . . Adaptive mesh : time t = 11.25s. The colours correspond to the numerical density of entropy production S. . . . . . . . . . . . . . . . . . . . . Figure 8.13 Adaptive mesh : time t = 13.25s. The colours correspond to the numerical density of entropy production S. . . . . . . . . . . . . . . . . . . . . Figure 8.14 Adaptive mesh : time t = 16s. The colours correspond to the numerical density of entropy production S. . . . . . . . . . . . . . . . . . . . . Figure 8.15 Adaptive mesh : time t = 17.5s. The colours correspond to the numerical density of entropy production S. . . . . . . . . . . . . . . . . . . . . Figure 11.1 Refinement and coarsening process for one triangular and quadrangular element. . . . . . . . . Snapshot of hydraulic head contours at t = 57.5 s for the end of the first backwash. Most velocities are less than 1 cm/h. . . . . . . . . . . . . . . Figure 12.7 Snapshot of hydraulic head contours at t = 67.5 s for the next incoming uprush. Most velocities are less than 0.8 cm/h. . . . . . . . . . . . . . . . Figure 15.1 Geometric set-up. . . . . . . . . . . . . . . . . . . Convergence rate of the L 1 discretisation error obtained with the refinement level l = 1, 2, 3 and 4. . Table 12.1 Parameters of constitutive laws for the idealized beach case [156]. . . . . . . . . . . . . . . . . . . . . . Table 18.1 Convergence rate of the L 2 error for ε = 0. The order is computed through the first order interpolation polynomial . . . . . . . . . . . . . . . . . . . . . . . .

	Figure 8.12

.6 Experimental settings. . . . . . . . . . . . . . . . . Figure 8.7 Numerical result for the water height at time t = 11.25s. The colours correspond to the magnitude of the kinetic energy. . . . . . . . . . . . . . . . . Figure 8.8 Figure 12.1 Results for the one-dimensional infiltration numerical test case. . . . . . . . . . . . . . . . . . . . Figure 12.2 Numerical simulation of two-dimensional drainage with water table position compared to experimental data from [157]. . . . . . . . . . . . . . . . Figure 12.3 Topography of the idealized beach and mesh used for the simulation. . . . . . . . . . . . . . . . Figure 12.4 Snapshot of hydraulic head contours at t = 12.5s during the first wave propagation on calm waters. . . . . . . . . . . . . . . . . . . . . . . . . Figure 12.5 Snapshot of hydraulic head contours at t = 37.5 s for the end of the first uprush. Velocities range from about 0.3 to 0.6 cm/h in the beach. Maximum velocity is 35 cm/h for infiltration. . . . . . Figure 12.6 n := max 0≤i≤N+2 (h n i ). . . . . . . . . . . . . . . Figure 18.3 Influence of σ. . . . . . . . . . . . . . . . . . . . . xvii L I S T O F TA B L E S

Table 18 . 2

 182

  où x, y dénote la variable spatiale horizontale et verticale, ρ est la densité, (u, v) est le vecteur vitesse avec u (resp. v) la composante horizontale (resp. verticale), p est une loi de pression donnée par la loi d'état p(ρ) = c 2 ρ. La constante c est définie par c 2 = RT où R

constante universelle des gaz parfaits de l'air et T la température, supposée constante. Le profil de viscosité est choisi sous la forme stratifiée suivante ν 1 (t, x, y) = ν 0 e -g/c 2 y où ν 0 ∈ R et ν 2 est une fonction donnée.

  la hauteur d'eau, q = hu est le débit d'eau, q b est le débit de sédiments (ou le flux de transport solide) et ξ = 1/(1 -

	ψ) est un coefficient de porosité associé à la couche de sédiments.
	Il existe principalement trois modes de transport : par suspension,
	par charriage et/ou par saltation. Chaque mode est régi par des lois
	empiriques valable pour une plage de donnée. Par conséquent, il
	est difficile d'établir une loi générale décrivant les trois modes de
	transport. La réduction d'un modèle avec une loi générale (issue d'une
	description microscopique, par exemple) est un des enjeux majeurs
	pour la modélisation de la sédimentation.

  , x) est le champ de vitesses et ρ w (t, x) la densité. Le terme p(t, x) est une pression barotrope non homogène tenant en compte l'effet des sédiments dans le fluide. Le terme σ(ρ w , u) est le tenseur de contrainte -pI 3 + 2Σ(ρ w ) : D(u) + λ(ρ w )div(u) I 3 où I 3 correspond à la matrice identité et Σ(ρ w ) est un tenseur de viscosité anisotrope.L'interaction fluide/sédiments est décrite au niveau du terme de pression et du terme source. Le terme de pression s'écrit p(t, x)

  As a well-known result, the uniqueness of the (weak) solution of (7.1) is lost even if the initial data are smooth. It can be recovered by completing Equation (7.1) with Equation (7.3). Numerical approximation of Equations (7.1) with Equation (7.

Remark 7.1.1.

  2 N-1 } do Let j be the biggest integer such that 2 j divides i foreach interface x k+1/2 such that L k+1/2 Nj do 1. compute the integral of F k+1/2 (t) on the time interval 2 N-L k+1/2 δt n , 2. distribute F k+1/2 (t n ) to the two adjacent cells, 3. update only the cells of level greater than Nj.

  7.3: Illustration of a global time advancement with a first order scheme at time t n + δt n . Figure 7.4: Illustration of a global time advancement with a first order scheme at time t n + 2δt n . Figure 7.5: Illustration of a global time advancement with a first order scheme at time t n + 3δt n . Figure 7.6: Illustration of a global time advancement with a first order scheme at time t n + 4δt n .

  which shows that this scheme is locally non conservative. At the next intermediate step, the conservativity is preserved, i.e.

Table 8 . 1 :

 81 ex || 1 Convergence rate of the L 1 discretisation error obtained with the refinement level l = 1, 2, 3 and 4.

	Criterion 0	2.0136	2.1342
	Criterion 1	2.0491	2.0901
	Criterion 2	2.1109	2.1789

Table 12

 12 .1.

	Parameters	Hydraulic conductivity Water content
	A/C	2.99 × 10 6	40 000
	B/D	5.0	2.90
	Value at saturation	35.0	0.30

Table 12 .

 12

1: Parameters of constitutive laws for the idealized beach case

[START_REF] Vauclin | Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem[END_REF]

.

  ∂ t v + u∂ x v + v∂ y v + w∂ z v .

			σ
				2
	+P β P β = P α -µ 1 ∂β ∂z ∂α -P α ∂z	Dv Dt
				(16.16)
	where	Dv Dt	stands for the material derivative
	Dv Dt := asymptotic expansions of the width-averaged terms in
	system (16.16):
	Thanks to the expressions (16.14) and (16.15), the average of the
	terms in System (16.16) can be written

  σ) + O(µ 1 ) . , since ∂ z P |z=η = 0, we get ∂ y η = O(µ 1 ), i. e. , we obtain the flat surface approximation η(t, x, y) = η eq (t, x) + O(µ 1 )

	with
		A eq = Ω eq (t, x)			(16.27)
	On the other hand, coming back to System (16.2)-(16.5) together with thanks to
	Equation (16.23), integrating Equation (16.4) for s ∈ [α(x, z), y], we show that P(t, x, y, z) = P α (t, x, z) -µ 1 y α(x,z) Dt (16.24) ds Dv A(t, x) = Ω(t,x) dy dz = y + (t,x)
		= P α (t, x, z) + O(µ 1 ) .
	As a consequence, the width-averaged pressure P can be approxi-
	mated by P ∂ ∂z w Next, integrating this equation for s ∈ [z, η(t, x, y)] using the boundary P α (t, x, z) = -1 F r 2 -µ 2 D Dt
	condition (16.6) with Equation (16.24) lead to
	P α (t, x, z) =	η(t, x, y) -z F r 2	+ µ 2	η(t,x,y) z	D Dt	w α (t, x, s) ds + O(µ 1 ) .
	Thus, taking the y-derivative of this expression yields to
	0 = ∂ y η w Consequently(16.25) 1 F 2 r + µ 2 D Dt
	for some function η eq defined hereafter. It means that one can neglect
	the y-variations of the free surface of the three-dimensional model
	(16.2)-(16.5) (see Figure 16.1). In other words, the fluid cross-section
	Ω(t, x) (15.3) and the wet area A(t, x) (15.4) can be simplified in
	Ω (16.26)

α at order O(µ 1 ), i. e. P (t, x, z) = P α (t, x, z) + O(µ 1 ) .

asymptotic expansion of the free surface:

Using Equations (16.15) and

(16.24)

in Equation (

16

.5), we can write the z-gradient of the pressure as α (t, x, z)

+ O(µ 1 ) . α |z=η + O(µ 1 ) = -∂ y η ∂ z P |z=η + O(µ 1 ) . eq (t, x) = (y, z) ∈ R 2 ; α(x, z) ≤ y ≤ β(x, z)

and d * (x) ≤ z ≤ η eq (t, x)

  In what follows, we note f b (t, x) = f α (t, x, d

	These equations are completed with the irrotational condition (16.22)
	and the following boundary conditions
	∂η eq ∂t	+ u α	∂η eq ∂x	= w α + O	µ 1 µ 2	and P α = O(µ 1 ) on z = η eq
					u α ∂ x d = w α + O	µ 1 µ 2	on z = d
	which are obtained through Equations (16.6) and (16.7) using Equa-
	tions (16.13), (16.14), (16.15).
	These end the first step.
	16.3 2d-1d like model reduction and asymptotic expan-
	sions

* (x)) for a given function f and S S(x, z) = z d * (x)

  ūeq , x, s) ds dz + O(µ 2 2 ) .

	(16.37)
	Finally, thanks to (16.10), (16.34), (16.37), by linearity of S (16.35),
	we obtain
	u(t, x, y, z)

  Computation of the term Ω eq (t,x) a 1 dy dz Computation of the term Ω eq (t,x) a 2 dy dz

	17.2.2 Thanks to the asymptotic approximation of the term u (16.38) up to
	order O(µ
	We have
	Ω eq (t,x)	∂ t (u)dy dz = ∂ t	Ω eq (t,x)	u dy dz -	∂Ω eq (t,x)	u∂ t m • n ds . (17.5)

17.2.1

  17.2.2.1Computation of the term Ω eq (t,x) a 3 dy dz Computation of the term Ω eq (t,x) a 1 + a 2 + a 3 dy dz

	We have
	Ω eq (t,x)	div y,z [uv] dy dz =	∂Ω eq (t,x)	uv • n ds .	(17.7)
	17.2.3 Gathering the results of the computations (17.5), (17.6) and (17.7),
	using the boundary conditions (17.1), we get

  Computation of the term Ω eq (t,x) a 4 dy dzWe apply the Leibniz rule to the gradient of the pressure and we obtain

							.8)
	17.2.4 Ω eq (t,x)	∂ x P dy dz = ∂ x	Ω eq (t,x)	P dy dz -	∂Ω eq (t,x)	P∂ x m • n ds (17.9)
	where
		∂Ω eq (t,x)	P∂ x m • n ds =	Γ wb (x)	P∂ x m • n ds
	thanks to the free surface boundary condition (15.8).
	To compute the above boundary integral, we consider the parametri-
	sation
		z

  pressure source terms readI 2 (t, x) = σ (x) h eq (t, x) 2 nh (t, x, z)σ (x) dzσ(x)P nh (t, x, d(x))d (x) .

		2F r	2	-σ(x)	h eq (t, x) F r 2	d (x)	(18.1)
	and
	DI 2 (t, x) =	η eq (t,x) d(x)
							(18.2)

P

  Let us define d i+1/2 = max(d i , d i+1 ) where d i = 1 = max(σ i , σ i+1 ) where σ i = 1

				δx m i	d(x)dx,
	σ i+1/2 δx m i the reconstructed states	σ(x)dx and let us define
	A -i+1/2 = σ i+1/2	A i σ i	+ d i -d i+1/2 ,
	A + i+1/2 = σ i+1/2	A i+1 σ i+1	+ d i+1 -d i+1/2
	with
	U -i+1/2 = (A -i+1/2 , A -i+1/2 u i ), U + i+1/2 = (A + i+1/2 , A + i+1/2 u i+1 )

  numη exact || 2 ||η numη exact || ∞

	100	0.0789	0.0449
	200	0.0497	0.0288
	400	0.0304	0.0180
	800	0.0198	0.0116
	1600	0.0153	0.0081
	3200	0.0138	0.0062
	Order	0.53	0.58

Table 18 . 1 :

 181 Convergence rate of the L 2 error for ε = 0. The order is computed through the first order interpolation polynomialN η numη ref 2 η numη ref ∞

	100	0.05212	0.02533
	200	0.02096	0.01082
	400	0.01079	0.00554
	800	0.00748	0.00503
	1600	0.00635	0.00412
	3200	0.00505	0.00300
	Order	0.64	0.56

Table 18 . 2 :

 182 Convergence rate of the L 2 error for ε = 0.4. The reference solution is computed with 10 000 cells. The order is computed through the first order interpolation polynomial

travaux de recherche : période postdoctorale

We have used the standard Godunov scheme to compute each numerical solution on a fine grid with N = 10000 and free boundary conditions.

11.2 nonlinear solver with adaptive time stepping

Let us note that recent progress in[START_REF] Zeng | A Newton's Iteration Converges Quadratically to Nonisolated Solutions Too[END_REF] shows that the Jacobian rank deficiency is not an obstacle anymore for the convergence of Newton-Raphson methods. Even stronger, the quadratic convergence may be preserved for rank-deficient systems.

For the sake of simplicity and clarity, we have presented the finite volume method using the Rusanov solver but the method is not limited to this one.

[START_REF] Dougalis | Theory and numerical analysis of Boussinesq systems: A review[END_REF] demonstrated all the advantages and accuracy of the continuous Galerkin methods for the nonlinear dispersive wave propagation.

mois. Elle m'a aidé à regarder dans la bonne direction, à voir le bon côté de la vie après

A N E W N O N -L I N E A R D I S P E R S I V E M O D E L In this section,we assume that the flow is governed, on the spacetime domain Ω eq = 0≤t≤1 Ω eq (t), instead of Ω (see Def. (15.2)), by the three-dimensional incompressible and irrotational Euler equations (16.2)-(16.5), (16.9) where the wet region is defined as the region in which the fluid resides at each time t ∈ [0, 1], Ω eq (t) = 0≤x≤1 Ω eq (t, x) with Ω eq (t, x) the fluid cross-section defined by (16.26).

These equations (16.2)-(16.5) are completed with the kinematic boundary condition at the free surface (16.6) Γ fs (t, x) = {z ∈ R; z = η eq (t, x)}, and a no-penetration condition on the wet boundary (16.7)

These boundary conditions can be written under the following form:

where m(t, x, y, z) = (y,

stands for the outward unit normal vector in the Ω eq (t, x)-plane (see Figure 16.1) and v = v w .

To work with the wet region, we introduce its indicator function

x) for all t, x, y, z ∈ R where

The function Φ is advected by the flow so its material derivative, with respect to the flow u, must therefore be zero. Moreover, thanks to The kinematic free surface condition and the wet boundary condition are the same, except the definition of the outward unit normal vector which is now

Thus, we have now

Almost all the computations are the same to the exception of the computation of the Ω eq (t,x) a 4 dy dz which yields to different expression of the hydrostatic pressure source term (17.12) and the non hydrostatic pressure source term (17.15). Indeed, since we have now d * (x) = d(x) and σ = σ(x), thus A eq (t, x) = σ(x)h eq (t, x). Applying successively the one-dimensional Leibniz integral rule in the y-direction and in the z-direction, we obtain Ω eq (t,x)