Les calottes polaires, ou inlandsis, sont parmi les principaux contributeurs à la montée des océans. Ces systèmes dynamiques gagnent de la masse par accumulation de neige, et en perdent par fonte au contact de l'océan et à la surface, ainsi que par le vêlage d'icebergs. Depuis plus de trois décennies, les observations ont montré que les calottes polaires de l'Antarctique et du Groenland perdent plus de masse qu'ils n'en gagnent. L'évolution des glaciers suite à ce déséquilibre de masse est devenue aujourd'hui l'une des problématiques les plus importantes des implications du changement climatique. Le Groupe d'experts intergouvernemental sur l'évolution du climat (GIEC) a identifié la contribution des glaciers comme l'un des facteurs clés d'incertitude de prédiction de l'élévation du niveau des mers. La modélisation numérique est le seul outil efficace pour répondre à cette question. Cependant, modéliser l'écoulement de glace à l'échelle du Groenland ou de l'Antarctique représente un défi à la fois scientifique et technique. Deux aspects clés de l'amélioration de la modélisation des glaciers sont abordés dans cette thèse. Le premier consiste à déterminer certaines propriétés non mesurables de la glace par méthode inverse. La friction ou la rigidité des barrières de glace, sont des paramètres qui ne peuvent être mesurés directement et doivent donc être déduits à partir d'observations par télédétection. Nous appliquons ici ces inversions pour trois modèles d'écoulement de glace de complexité croissante: le modèle bidimensionnel de MacAyeal/Morland, le modèle dit d'ordre supérieur de Blatter/Pattyn et le modèle full-Stokes. Les propriétés ainsi calculées sont ensuite utilisées pour initialiser des modèles grande-échelle et pour déterminer le degré de complexité minimum nécessaire pour reproduire correctement la dynamique des glaciers. Le second aspect abordé dans ce travail est l'amélioration de la consistance des données pour la modélisation numérique. Les données disponibles sont souvent issues de campagnes de mesures s'étalant sur plusieurs années et dont résolutions spatiales varient, ce qui rend leur utilisation pour des simulations numériques difficiles. Nous présentons ici un algorithme basé sur la conservation de la masse et les méthodes inverses pour construire des épaisseurs de glace qui sont consistantes avec les mesures de vitesse. Cette approche empêche la redistribution artificielle de masse qu'engendrent généralement les autres méthodes de cartographie de l'épaisseur de glace, ce qui améliore considérablement l'initialisation des modèles d'écoulement de glace. Les avancées présentées ici sont des étapes importantes afin de mieux caractériser de manière précise les glaciers et de modéliser leur évolution de manière réaliste.

Les calottes polaires, ou inlandsis, sont parmi les principaux contributeurs à la montée des océans. Ces systèmes dynamiques gagnent de la masse par accumulation de neige, et en perdent par fonte au contact de l’océan et à la surface, ainsi que par le vêlage d’icebergs. Depuis plus de trois décennies, les observations ont montré que les calottes polaires de l’Antarctique et du Groenland perdent plus de masse qu’ils n’en gagnent. L’évolution des glaciers suite à ce déséquilibre de masse est devenue aujourd’hui l’une des problématiques les plus importantes des implications du changement climatique. Le Groupe d’experts intergouvernemental sur l’évolution du climat (GIEC) a identifié la contribution des glaciers comme l’un des facteurs clés d’incertitude de prédiction de l’élévation du niveau des mers. La modélisation numérique est le seul outil efficace pour répondre à cette question. Cependant, modéliser l’écoulement de glace à l’échelle du Groenland ou de l’Antarctique représente un défi à la fois scientifique et technique. Deux aspects clés de l’amélioration de la modélisation des glaciers sont abordés dans cette thèse. Le premier consiste à déterminer certaines propriétés non mesurables de la glace par méthode inverse. La friction ou la rigidité des barrières de glace, sont des paramètres qui ne peuvent être mesurés directement et doivent donc être déduits à partir d’observations par télédétection. Nous appliquons ici ces inversions pour trois modèles d’écoulement de glace de complexité croissante: le modèle bidimensionnel de MacAyeal/Morland, le modèle dit d’ordre supérieur de Blatter/Pattyn et le modèle full-Stokes. Les propriétés ainsi calculées sont ensuite utilisées pour initialiser des modèles grande-échelle et pour déterminer le degré de complexité minimum nécessaire pour reproduire correctement la dynamique des glaciers. Le second aspect abordé dans ce travail est l’amélioration de la consistance des données pour la modélisation numérique. Les données disponibles sont souvent issues de campagnes de mesures s’étalant sur plusieurs années et dont résolutions spatiales varient, ce qui rend leur utilisation pour des simulations numériques difficiles. Nous présentons ici un algorithme basé sur la conservation de la masse et les méthodes inverses pour construire des épaisseurs de glace qui sont consistantes avec les mesures de vitesse. Cette approche empêche la redistribution artificielle de masse qu’engendrent généralement les autres méthodes de cartographie de l’épaisseur de glace, ce qui améliore considérablement l’initialisation des modèles d’écoulement de glace. Les avancées présentées ici sont des étapes importantes afin de mieux caractériser de manière précise les glaciers et de modéliser leur évolution de manière réaliste. / Ice sheets are amongst the main contributors to sea level rise. They are dynamic systems; they gain mass by snow accumulation, and lose it by melting at the ice-ocean interface, surface melting and iceberg calving at the margins. Observations over the last three decades have shown that the Greenland and Antarctic ice sheets have been losing more mass than they gain. How the ice sheets respond to this negative mass imbalance has become today one of the most urgent questions in understanding the implications of global climate change. The Intergovernmental Panel on Climate Change (IPCC) has indeed identified the contribution of the ice sheets as a key uncertainty in sea level rise projections. Numerical modeling is the only effective way of addressing this problem. Yet, modeling ice sheet flow at the scale of Greenland and Antarctica remains scientifically and technically very challenging. This thesis focuses on two major aspects of improving ice sheet numerical models. The first consists of determining non-observable ice properties using inverse methods. Some parameters, such as basal friction or ice shelf hardness, are difficult to measure and must be inferred from remote sensing observations. Inversions are developed here for three ice flow models of increasing complexity: MacAyeal/Morland’s shelfy-stream model, Blatter/Pattyn’s higher order model and the full-Stokes model. The inferred parameters are then used to initialize large-scale ice sheet models and to determine the minimum level of complexity required to capture ice dynamics correctly. The second aspect addressed in this work is the improvement of dataset consistency for ice sheet modeling. Available datasets are often collected at different epochs and at varying spatial resolutions, making them not readily usable for numerical simulations. We devise here an algorithm based on the conservation of mass principle and inverse methods to construct ice thicknesses that are consistent with velocity measurements. This approach therefore avoids the artificial mass redistributions that occur in existing algorithms for mapping ice thickness, hence considerably improving ice sheet model initialization. The advances made here are important steps towards the ultimate objective of accurate characterization of ice sheets and the realistic modeling of their evolution.

Des schémas numériques précis et robustes sont proposés pour modéliser les écoulements souterrains et leur couplage avec le ruissellement surfacique. Les écoulements souterrains sont d´écrits par l'équation de Richards (instationnaire) qui est discrétisée par une méthode BDF en temps et une méthode de Galerkine discontinue à pénalisation intérieure symétrique en espace. Des cas tests sur des colonnes d'infiltration confirment la robustesse des schémas choisis. Dans un premier temps, nous considérons des conditions de Signorini pour l'équation de Richards afin de modéliser la présence de drains en fond d'aquifère ou l'affleurement de la nappe en négligeant le ruissellement, c'est-à-dire en supposant que l'eau exfiltrée est immédiatement évacuée du système. Dans un second temps, nous prenons en compte le ruissellement par le biais de conditions de couplage qui imposent l'égalité des flux d'eau échangés et la continuité de la pression à l'interface. Les écoulements superficiels sont d´écrits par l'équation de l'onde cinématique qui constitue une approximation des équations de Saint-Venant. L'équation de l'onde cinématique est discrétisée par une méthode de Godunov. Les deux schémas, pour l'écoulement souterrain et pour l'écoulement superficiel, sont conservatifs et peuvent être utilisés dans des algorithmes de couplage faisant intervenir un ou plusieurs pas de temps. Pour assurer la conservation de la masse d'eau totale du système couplé, les flux à l'interface doivent être convenablement choisis. Nous donnons en particulier la construction de ces flux pour les schémas BDF1 et BDF2. La précision et la robustesse de nos schémas sont évaluées sur plusieurs cas tests dont le drainage d'une lame d'eau, deux cas d'exfiltration de nappe (l'un provoqué par la pluie et l'autre par une injection en fond d'aquifère) et un ruissellement hortonien. Enfin, nous présentons une application concrète portant sur le fonctionnement hydrologique d'un petit bassin versant drainé.

La méthode de projection et l'approche variationnelle de Sasaki sont deux techniques permettant d'obtenir un champ vectoriel à divergence nulle à partir d'un champ initial quelconque. Pour une vitesse d'un vent en haute altitude, un champ de vitesse sur une grille décalée est généré au-dessus d'une topographie donnée par une fonction analytique. L'approche cartésienne nommée Embedded Boundary Method est utilisée pour résoudre une équation de Poisson découlant de la projection sur un domaine irrégulier avec des conditions aux limites mixtes. La solution obtenue permet de corriger le champ initial afin d'obtenir un champ respectant la loi de conservation de la masse et prenant également en compte les effets dûs à la géométrie du terrain. Le champ de vitesse ainsi généré permettra de propager un feu de forêt sur la topographie à l'aide de la méthode iso-niveaux. L'algorithme est décrit pour le cas en deux et trois dimensions et des tests de convergence sont effectués. / The Projection method and Sasaki's variational technique are two methods allowing one to extract a divergence-free vector field from any vector field. From a high altitude wind speed, a velocity field is generated on a staggered grid over a topography given by an analytical function. The Cartesian grid Embedded Boundary method is used for solving a Poisson equation, obtained from the projection, on an irregular domain with mixed boundary conditions. The solution of this equation gives the correction for the initial velocity field to make it such that it satisfies the conservation of mass and takes into account the effects of the terrain. The incompressible velocity field will be used to spread a wildfire over the topography with the Level Set Method. The algorithm is described for the two and three dimensional cases and convergence tests are done.

La méthode de projection et l'approche variationnelle de Sasaki sont deux techniques permettant d'obtenir un champ vectoriel à divergence nulle à partir d'un champ initial quelconque. Pour une vitesse d'un vent en haute altitude, un champ de vitesse sur une grille décalée est généré au-dessus d'une topographie donnée par une fonction analytique. L'approche cartésienne nommée Embedded Boundary Method est utilisée pour résoudre une équation de Poisson découlant de la projection sur un domaine irrégulier avec des conditions aux limites mixtes. La solution obtenue permet de corriger le champ initial afin d'obtenir un champ respectant la loi de conservation de la masse et prenant également en compte les effets dûs à la géométrie du terrain. Le champ de vitesse ainsi généré permettra de propager un feu de forêt sur la topographie à l'aide de la méthode iso-niveaux. L'algorithme est décrit pour le cas en deux et trois dimensions et des tests de convergence sont effectués. / The Projection method and Sasaki's variational technique are two methods allowing one to extract a divergence-free vector field from any vector field. From a high altitude wind speed, a velocity field is generated on a staggered grid over a topography given by an analytical function. The Cartesian grid Embedded Boundary method is used for solving a Poisson equation, obtained from the projection, on an irregular domain with mixed boundary conditions. The solution of this equation gives the correction for the initial velocity field to make it such that it satisfies the conservation of mass and takes into account the effects of the terrain. The incompressible velocity field will be used to spread a wildfire over the topography with the Level Set Method. The algorithm is described for the two and three dimensional cases and convergence tests are done.

L’objectif de notre étude est de modéliser un réservoir géothermique. Si nous supposons que le réservoir géothermique n’est composé que d’eau pure le transfert de matière et d’énergie est classiquement décrit par deux équations de conservation : la conservation de la matière et la conservation de l’énergie. À ces deux équations vient s’ajouter la vitesse du fluide classiquement donnée par la loi de Darcy tandis que les propriétés thermodynamiques, obtenues grâce à des équations théoriques ou empiriques (les équations d’état), ferment le modèle mathématique. Dès lors, ce modèle fermé, il existe différents schémas de résolution. Le premier est de résoudre en pression et température puis de procéder à un changement de variables lors du passage de monophasique à diphasique ou de diphasique à monophasique. TOUGH2 utilise le couple pression-saturation dans la zone diphasique.La seconde approche est de résoudre en pression et enthalpie afin d’accroître la stabilité lors de la transition entre l’état monophasique et l’état diphasique (voir Hydrotherm). Nous avons adopté la seconde option, résoudre en pression et enthalpie. De plus la résolution spatiale est faite avec les volumes finis.La modélisation d’un réservoir géothermique fait intervenir des équations fortement dépendantes l’une de l’autre. Cependant nous avons fait le choix de découpler la résolution afin de se libérer de la complexité de la résolution du système couplé. En effet, cette méthode possède l’avantage d’être moins consommatrice de mémoire puisque nous travaillons toujours avec le même nombre de données, mais dans une matrice deux fois moins importante. Nous montrerons que cette méthode demeure suffisamment précise pour une utilisation aussi bien dans le domaine industriel que dans celui de la recherche.Nous offrons à l’utilisateur une grande liberté grâce à l’implémentation de plusieurs méthodes : Euler implicite, explicite, Runge-Kutta ou BDF2 pour les solveurs temporels ou GMRES et BICGSTAB pour les solveurs linéaires. Nous pouvons gérer des conditions aux limites très variées telles que des flux nuls (décrivant une frontière qui n’échange pas de matière avec l’extérieur) ou une condition mixte (un Dirichlet sur la pression et un Dirichlet ou condition « sortie libre » sur la température… Cette dernière situation décrit une zone de recharge ou de décharge. Nous avons développé un outil multilangage : Python,Fortran et C++ (une implémentation de l’IAPWS provenant du projet freesteam incluant la zone supercritique). Tous ces langages sont orientés objet. L’IAPWS est l’outil permettant de calculer les propriétés physiques inconnues et par conséquent il ferme le système.Enfin nous avons appliqué le modèle sur le bassin parisien, France, sur plusieurs systèmes 1D et un autre système 2D réalisés par Coumou avec la plateforme CSMP++. Le bassin parisien est un réservoir exploité pour produire de la chaleur par le biais du pompage d’une eau à 70 _C et réinjecté à 40 _C. Les simulations 1D permettent de visualiser le déplacement d’un front de chaleur en haute enthalpie. La simulation 2D montre la convection naturelle de l’eau dans une faille. Chaque simulation a été comparée aux résultats obtenus avec un autre code (CSMP++, HYDROTHERM ou TOUGH2) et les résultats sont en accord. / The purpose of our study is to model a geothermal reservoir. When geothermal reservoirs are assumed to be composed of pure water, the transfer of mass and energy is classically described by two balance equations: the mass balance equation and the energy balance equation. In addition to those equations, fluid velocity is classically given by the Darcy law while thermodynamic properties, inferred from theoretical or empirical equations of state, are used to close the mathematical system. Once this system is closed, there exist different solutions. The first one is to solve for pressure and temperature with a variable switch to saturation in the two-phase region (e.g. TOUGH2). The second one is to solve for pressure and enthalpy to increase the stability of phase transition between single and two-phase states (e.g. Hydrotherm). We adopted the second option. We solve the system by using a splitting method — to get rid of the complexity of coupling equations — and a finite volume method for the spatial discretization. We offer some freedom to users thanks to the implementation of several methods like explicit or implicit Euler, Runge-Kutta or BDF2 for time solvers or GMRES and BICGSTAB for the linear solver. We can handle several boundary conditions like no-flow — describing a boundary which cannot exchange matter withthe exterior — or like a mixed-therm condition — a Dirichlet condition to the pressure and a Dirichlet or an outflow condition to the temperature in order to describe a recharge or a discharge zone — …Selecting object-oriented languages, we developed a multi-language framework, combining Python, Fortran and a C++ implementation of IAPWS (from the freesteam project) including the supercritical equations. To close the system physical propertiesare determined by the IAPWS- IF97 thermodynamic formulation. We’ve applied this simulation model to the dogger in Paris, France, to several onedimensional systems and a two-dimensional one made by Coumou with the CSMP++platform. The dogger is a reservoir exploited to produce heat by pumping water at 70 _C and reinjecting it in the reservoir at 40 _C. In the one-dimensional systems we wanted to observe the process of heat transfer from a higher temperature boundary to a smaller one in a high-energy domain. The last simulation shows the natural convection of water in a fault. For every simulation we compared the solutions we found with another code (TOUGH2 or CSMP++) and they agreed.