Keywords: glaciology, glaciers, ice sheet flow, numerical modeling, hybrid model, Tiling method, multi-model

Ice flow numerical models are essential for predicting the evolution of ice sheets in a warming climate. Recent research emphasizes the need for higher-order and even full-Stokes flow models instead of the traditional Shallow-Ice Approximation whose assumptions are not valid in certain critical but spatially limited areas. These higher-order models are however computationally intensive and difficult to use at the continental scale. The purpose of this work, therefore, is to develop a new technique that reduces the computational cost of ice flow models while maximizing their accuracy. To this end, several ice flow models of varying order of complexity have been implemented in the Ice Sheet System Model, a massively parallelized finite element software developed at the Jet Propulsion Laboratory. Analysis and comparison of model results on both synthetic and real geometries shows that sophisticated models are only needed in the grounding line area, transition between grounded and floating ice, whereas simpler models yield accurate results in most of the model domain.

There is therefore a strong need for coupling such models in order to balance computational cost and physical accuracy. Several techniques and frameworks dedicated to model coupling already exist and are investigated. A new technique adapted to the specificities of ice flow models is developed: the Tiling method, a multi-model computation strategy based on the superposition and linking of different numerical models. A mathematical analysis of a mixed Tiling formulation is first performed to define the conditions of application. The treatment of the junction between full-Stokes and simpler models that decouple horizontal and vertical equation is then elaborated in order to rigorously combine all velocity components. This method is finally implemented in the Ice Sheet System Model to design hybrid models that combine several ice flow approximations of varying order of complexity. Following a validation on synthetic geometries, this method is applied to real cases, such as Pine Island Glacier, in West Antarctica, to illustrate its relevance. Hybrid models have the potential to significantly improve physical accuracy by combining models in their domain of validity, while preserving the computational cost and being compatible with the actual computational resources.

First and foremost, I would like to express my sincere gratitude to Eric Rignot for his continuous support during my PhD, for his motivation and enthusiasm. He always advised me, shared his interesting ideas and had the patience to review the countless versions of my first paper.

I am indebted to Eric Larour for offering hundred hours of his time to teach me computer science and the art of developing a beautiful code, for the endless discussions and for the genuine interest in sharing his knowledge.

I am very grateful to Hachmi Ben Dhia who accepted to supervise this PhD; this research would not have been possible whithout his support. Despite the very long distance and the difficulties, he has always shared his ideas and advice. Such a work would never have been possible without our meetings in Paris, Los Angeles, Mineapolis and Saint-Bon (where glaciers are also retreating).

I also thank Denis Aubry who accepted to co-supervise my thesis and has always been very supportive.

I am greateful to my colleagues at JPL for the great moments we shared together. Special thanks to Ala and Nicole for the time they spent correcting my English.

I am thankful to the Ecole Centrale Paris, laboratoire MSSMat and all their staff who always found solutions to my rather special situation and made this project possible.

I thank Caltech's Jet Propulsion Laboratory and NASA for providing the all facilities and funds for this project.

I am grateful to my parents for the understanding and endless support throughout my education.

Above all, I would like thank Mathieu for sharing this experience with me and holding my hand through the hard times.

Résumé

La modélisation numérique des écoulements de glace est indispensable pour prédire l'évolution des calottes polaires suite au réchauffement climatique. De récentes études ont souligné l'importance des modèles d'écoulement dits d'ordre supérieur voir même de Stokes au lieu de la traditionnelle approximation de couche mince dont les hypothèses ne sont pas valables dans certaines zones critiques mais à l'étendue limitée. Cependant, ces modèles d'ordre supérieur sont difficiles à utiliser à l'échelle d'un continent en raison de leurs temps de calculs prohibitifs. Ce travail de thèse propose une nouvelle technique qui permet de réduire les temps de calculs tout en maximisant la précision des modèles. Plusieurs modèles d'écoulement de glace de complexité variables ont été mis en place dans ISSM (Ice Sheet System Model), un code élément fini massivement parallèle développé par le Jet Propulsion Laboratory. L'analyse et la comparaison des différents modèles, à la fois sur des cas théoriques et réels, montrent que l'utilisation des modéles les plus complets est principalement nécessaire au voisinage de la zone d'échouage, transition entre les parties flottantes et posées de la glace, mais aussi que des modèles plus simples peuvent être utilisés sur la majeure partie des glaciers. Coupler différents modèles présente donc un avantage significatif en terme de temps de calcul mais aussi d'amélioration de la physique utilisées dans les modèles. Plusieurs méthodes de couplage de modèles existent et sont présentées dans ce manuscrit. Une nouvelle technique, dite de tuilage, particulièrement adaptée au couplage de modèles d'écoulement de glace est décrite ici : son principe repose sur la superposition et le raccordement de plusieurs modèles mécaniques. Une analyse mathématique est effectuée afin de définir les conditions d'utilisation de cette méthode de tuilage. Le traitement du couplage entre un modèle de Stokes et des modèles simplifiés, pour lesquels le calcul des vitesses horizontales et verticales est découplé, est ensuite présenté. Cette technique a été mise en place dans ISSM afin de pouvoir créer des modèles hybrides combinant plusieurs modèles d'écoulement de complexité variable. Après avoir été validée sur des cas synthétiques, cette technique est utilisée sur des glaciers réels comme Pine Island Glacier, dans l'Antarctique de l'Ouest, afin d'illustrer sa pertinence. Les modèles hybrides ont le potentiel d'améliorer la précision des résultats en combinant différents modèles mécaniques, utilisés chacun dans les zones où leurs approximations sont valides, tout en réduisant les temps de calcul et en étant compatibles avec les ressources informatiques actuelles.

Mots-clés: glaciologie, glaciers, écoulement des calottes polaires, modélisation numérique, modèle hybride, méthode de tuilage, multi-modèle. Ice covers around 10% of the earth land surface [START_REF] Paterson | The Physics of Glaciers[END_REF] and is present on every continent. During the last 65 million years, the extent of ice coverage has considerably changed over time. Ice sheets are relatively recent features that appeared in Antarctica as climate changed from temperate to polar about 35 million years ago. The Antarctic ice sheet formed about 34 million years ago [START_REF] Mayewski | State of the Antarctic and Southern ocean climate system[END_REF]. This ice sheet was dynamic and fluctuating in response to variations in the Earth orbit. Only about 14 millions years ago did the Antarctic ice sheet become persistent, thicker, cooler and somewhat similar to the actual ice sheet [START_REF] Flower | The middle Miocene climatic transition -East Antarctic icesheet development, deep-ocean circulation and global carbon cycling[END_REF]. It is thought to have persisted through the early Pliocene warming. The Greenland ice sheet developed about 7 million years ago [START_REF] Mayewski | State of the Antarctic and Southern ocean climate system[END_REF], and the first ice sheets in northern Europe and North America occured around 2.5 million years ago, following a global cooling around 3 million years ago [START_REF] Shackleton | Oxygen isotope calibration of the onset of ice-rafting and history of glaciation in the North-Atlantic region[END_REF].

Over the last million years, the Earth regularly alternated between glacial, where ice covered most of the Northern Hemisphere, and interglacial periods [START_REF] Mayewski | State of the Antarctic and Southern ocean climate system[END_REF][START_REF] Cuffey | The Physics of Glaciers[END_REF]. Studies of ice cores [START_REF] Petit | Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica[END_REF]EPICA Community members, 2004, 2006] reveal the response of ice sheets to modifications in insolation patterns caused by orbital forcing and the association between atmospheric greenhouse gases and temperatures. Glacial and interglacial periods alternated over periods of about 100,000 years, consisting of long glacial cycles lasting about 90,000 years with temperatures 5 • C colder and sea level 120 meters lower than present, and short warm interglacial periods with a sea level close to the present one, lasting about 10,000 years [START_REF] Mayewski | State of the Antarctic and Southern ocean climate system[END_REF]. The warm events were more pronounced and shorter during the last 500,000 years than during the previous 500,000 ones [START_REF] Petit | Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica[END_REF]. Climatic events happening in the north and south polar regions seems to be correlated even when warm events in Antarctica preceded those in Greenland [EPICA Community members, 2006].

The last 50 years

Mass balance of a glacier is defined by the difference between the mass inputs and losses. The mass inputs are caused by surface accumulation due to snow or ice precipitations, snow drifting or any solid deposition (like subsurface accumulation caused by water freezing or basal accretion in the case of ice shelves). Mass losses are caused by melting, surface sublimation and ice calving [START_REF] Paterson | The Physics of Glaciers[END_REF].

To evaluate the mass balance of large glaciers and ice caps, several methods exist. The first possibility is to use the mass flux methods that estimate the difference between accumulation and depth-averaged ice flux at gates (usually the grounding line) during a certain amount of time [START_REF] Whillans | Mass balance of Ice Stream B, West Antarctica[END_REF]]. The second one consists in performing repeated measurements of the surface elevations with airborne or spaceborne altimeters [START_REF] Zwally | Mass changes of the Greenland and Antarctic ice sheets and shelves and contributions to sea-level rise: 1992-2002[END_REF] and then converting this volume change into a mass change by considering changes in firn compaction and crustal isostatic adjustment. Another possibility is to use spaceborne gravimetry [START_REF] Velicogna | Greenland mass balance from GRACE[END_REF] to measure spatial and temporal variations in the Earth's gravity field with tandem satellites.

Recent observations point out an acceleration of ice loss in Greenland with a widespread acceleration of glaciers below 66 • north between 1996 and 2000, which expandeded northward in 2005, resulted in an estimated increase in mass loss from 83 Gt/yr in 1996 to 183 Gt/yr in 2005 using the mass flux method [START_REF] Rignot | Changes in the velocity structure of the Greenland ice sheet[END_REF]. Gravity measurements between 2002 and 2009 conclude to a similar acceleration of mass loss from 104 to 246 Gt/yr [START_REF] Velicogna | Increasing rates of ice mass loss from the Greenland and Antarctic ice sheets revealed by GRACE[END_REF]. Combining these observations, [START_REF] Rignot | Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise[END_REF] estimate the ice discharge rate of the Greenland Ice Sheet between 1992 and 2009 to 250 ± 40 Gt/yr and its acceleration to 18 Gt/yr 2 .

In West Antarctica, ice discharge from the main glaciers (Pine Island, Thwaites and Smith glaciers) is estimated to be around 227 Gt/yr and the mass loss has increased 170% following the acceleration of these glaciers, or from 39 ± 15 Gt/yr to 105 ± 27 Gt/yr between 1996 and 2007 [START_REF] Thomas | Accelerated sea-level rise from West Antarctica[END_REF][START_REF] Rignot | Changes in West Antarctic ice stream dynamics observed with ALOS PALSAR data[END_REF]. Inflow of warm ocean waters is thought to have caused these modifications in ice dynamics [START_REF] Shepherd | Warm ocean is eroding West Antarctic Ice Sheet[END_REF][START_REF] Payne | Recent dramatic thinning of largest West Antarctic ice stream triggered by oceans[END_REF].

Major changes are also happening in the Antarctic Peninsula where thousand years old ice shelves are not only retreating but also experience rapid collapse in response to regional warming. Seven out of twelve ice shelves have significantly retreated or been almost entirely lost and 87% of the 244 marine glacier fronts have retreated during the past 61 years [START_REF] Cook | Retreating glacier fronts on the Antarctic Peninsula over the past half-century[END_REF][START_REF] Cook | Overview of areal changes of the ice shelves on the Antarctic Peninsula over the past 50 years[END_REF]. The collapse of these ice shelves have resulted in significant acceleration of their tributary glaciers. After the desintegration of Larsen A in 1995, its tributary glaciers experienced acceleration, thinning and retreat [START_REF] Rott | Northern Larsen Ice Shelf, Antarctica: further retreat after collapse[END_REF][START_REF] De Angelis | Glacier surge after ice shelf collapse[END_REF]. Similarly for the collapse of a large section of Larsen B ice shelf in 2002, the lower parts of four glaciers flowing into the collapsed portion of the shelf accelerated and began to thin within months after the collapse; two tributary glaciers that were flowing into the surviving part of the ice sheld experienced no changes [START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF][START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF].

East Antarctica seems to be less affected, which could be due to a modest increase in snowfall that could lead to a small thickening [START_REF] Alley | Ice-sheet and sea-level changes[END_REF]. However, two regions, 6

February 10, 2012

1.1. HISTORY OF GLACIERS AND ICE CAPS Totten and Cook glaciers, have significantly negative mass balance conditions, proving that East Antarctica is not immune to changes [START_REF] Rignot | Changes in ice dynamics and mass balance of the Antarctic ice sheet[END_REF]. Overall, [START_REF] Rignot | Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise[END_REF] estimate the average ice discharge rate in Antarctica between 1992 and 2009 to 200 ± 150 Gt/yr and its acceleration to 14 Gt/yr 2 .

An uncertain future

Without human intervention, a climate similar to the present one would extend well into the future, but the predicted increase in greenhouse gases makes this scenario very unlikely [EPICA Community members, 2004]. Anthropogenic effects are affecting the climate, and model suggests that over the 21st century, the Antarctic interior will warm by 3.4 ± 1 • C and sea ice extent will decrease by about 20% [START_REF] Mayewski | State of the Antarctic and Southern ocean climate system[END_REF].

Global warming is happening faster than anticipated by models projections. One of the main effects of this warming is the rise of sea level. IPCC-AR4 [2007] forecast the sea level rise to be less that 1 meter by 2100 in response to thermal extension of the ocean and the melting of glaciers at mid-latitude and in polar regions. But ice dynamics is not included in this estimate. If mountain glaciers melting and ocean dilation were the main causes of sea-level rise during the twentieth century, influx of fresh water from the Antarctic and Greenland ice sheets will soon overtake steric effects as the most important contributions to sea level rise.

A recent study by [START_REF] Rignot | Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise[END_REF] based on observations over the last decade concludes that the increase in sea level caused by ice discharge from the Antarctic and Greenland Ice Sheets will reach 56 cm if the current rate of acceleration in ice sheet loss (36.5 Gt/yr 2 ) remains constant during the next century. However a constant acceleration is unlikely. The potential increase in sea level due to the collapse of the West Antarctic Ice Sheet is estimated to 3.3 meters [Bamber et al., 2009b]. If all the ice of the Greenland Ice Sheet were to melt, the sea level rise following this event would reach 7 meters; for Antarctic Ice Sheet, the sea level rise would increase by 65 meters (5 meters for West Antarctica and 60 meters for East Antarctica).

Importance of ice dynamics

Several processes drive changes of the ice sheet evolution. They can be separated into surface mass balance and dynamic changes. Surface mass balance changes include modifications in the amount of precipitations or changes in the amount of precipitation that is effectively transformed into ice. Snowfalls tend to increase when climate warms in cold polar regions but changes can also be caused by shifts in atmospheric circulation [ [START_REF] Van Den Broeke | The semiannual oscillation and Antarctic climate, part 5: impact on the annual temperature cycle as derived from NCEP/NCAR re-analysis[END_REF]. Surface mass balance can be reduced by increased summer air temperature that enhance surface melt and decrease the amount of precipitation reaching the surface as snow rather than rain [START_REF] Cuffey | The Physics of Glaciers[END_REF].

Dynamic changes are caused by modifications in ice properties, glaciers' geometry or forces applied to the ice. Thinning and retreat of marine terminated ice streams reduce the buttressing forces and cause ice streams to accelerate [START_REF] Schmeltz | Sensitivity of Pine Island Glacier, West Antarctica, to changes in ice-shelf and basal conditions: a model study[END_REF]. Warming of oceans increases ice shelves bottom melting and lead to ice shelves thinning and grounding line retreat. [START_REF] Gagliardini | Coupling of ice-shelf melting and buttressing is a key process in ice-sheets dynamics[END_REF] showed that not only the melt rate but also its repartition influence the ice sheet response. Any mechanism that reduces the ice viscosity or increases the bed slipperiness tends to thin the ice sheet [START_REF] Cuffey | The Physics of Glaciers[END_REF]. Indirect mechanisms also affect the ice sheet evolution as glaciers and ice sheets interact with the rest of the Earth system.

Detailed and realistic modeling of the Antarctic and Greenland Ice Sheets is therefore needed to improve our understanding of the ice sheets evolution and to make projections of sea level rise in a warming climate. As pointed out by the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report [IPCC-AR4 , 2007], uncertainties associated to ice sheets evolutions dominate the projections of sea level rise. IPCC-AR4 [2007] forecasts are therefore conservative as they do not take into account the ice dynamics. Modeling the dynamics of an ice sheet consists not only of modeling the evolution of its velocity fields, but also the evolution of its temperature and geometry. Therefore, accurate modeling of ice sheets response requires a good representation of the physics of ice motion, well-constrained boundary conditions and computationally scalable software packages.

Several initiatives are being developed and funded to improve projections and produce reasonable estimates of ice sheet contribution to sea level. The two largest efforts are the European initiative ice2sea and the American effort SeaRISE. However, these programs rely on ice flow models that need to be improved to better include the physical processes that drive ice sheets evolution.

Ice flow models

Ice flow models are based on the momentum balance, mass continuity and heat budget equations. Ice flow models can be divided into two categories: flow line/flow band models that only include one horizontal dimension and planview/3d models that include the two horizontal dimensions. Flow line or flow band models requires lateral shear to be parameterize. For this reason, they are difficult to use to model dynamic changes and will not be discussed in this thesis. We will focus on planview and 3d models for this work.

Beginnings of ice sheet modeling

Ice sheet modeling started in the late 1970's with [START_REF] Mahaffy | 3-dimensional numerical model of ice sheets -tests on the barnes ice cap, northwest territories[END_REF] and [START_REF] Jenssen | A three-dimensional polar ice sheet model[END_REF] but was limited by the computational power available at that time. Following the increase in computational power and the rigorous establishment of the thermo-dynamic equations for ice flow [START_REF] Fowler | On the flow of polythermal glaciers. 2. Surface-wave analysis[END_REF][START_REF] Morland | A model for compaction and shear hysteresis in saturated granular-materials[END_REF], ice sheet models started to perform continental-scale simulations [START_REF] Huybrechts | Evolution of the East Antarctic ice sheet: a numerical study of thermo-mechanical response patterns with changing climate[END_REF][START_REF] Huybrechts | A 3-D model for the Antarctic ice sheet: a sensitivity study on the glacialinterglacial contrast[END_REF]. These models were mainly based on the Shallow Ice Approximation, also called shallow ice sheet approximation as it is a good approximation to model grounded ice with almost no sliding [Hutter , 1982b[START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF][START_REF] Morland | A model for compaction and shear hysteresis in saturated granular-materials[END_REF].

These models have been used to reconstruct glacial/interglacial cycles and understand the 8 February 10, 2012 1.2. ICE FLOW MODELS sensitivity of both Greenland Ice Sheet [Greve, 1997b;[START_REF] Ritz | Sensitivity of a greenland ice sheet model to ice flow and ablation parameters: Consequences for the evolution through the last climatic cycle[END_REF] and Antarctic Ice Sheet [START_REF] Pollard | Modelling West Antarctica ice sheet growth and collapse through the past five million years[END_REF] to changes. The majority of these models use the finite difference method to solve the thermo-dynamic equations on regular grids, even though finite element models are becoming more popular. The typical horizontal grid size vary between 20 and 50 km. However, if these models are able to perform paleoclimate reconstructions, the results after thousand of years of spin-up do not coincide with the present-day conditions accurately enough to use this type of initialization to model the future evolution of ice sheets for the next hundreds of years.

Recent efforts and initiatives

The Shallow Ice Approximation allows to correctly reproduce the motion of grounded ice with little sliding and is valid for the majority of areas on ice sheets, but breaks down in critical areas like the ice divides, the shear margins or the grounding line, and is not appropriate to model floating ice. For this reason, new approximations have been developed to model ice flow. The most common ones include a bidimensional vertically integrated model, called shallow ice stream approximation or shallow ice shelf approximation developed by [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment -Theory and application to Ice Stream-B, Antarctica[END_REF] and [START_REF] Morland | Plane and radial ice-shelf flow with prescribed temperature profile[END_REF] and the higher order three dimensional model proposed by [START_REF] Blatter | Velocity and stress-fields in grounded glaciers: A simple algorithm for including deviatoric stress gradients[END_REF] and [START_REF] Pattyn | A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes[END_REF]. Both models decouple the full-Stokes equations into horizontal and vertical equations, in order to limit the computational ressources required. Several recent models also solve the full-Stokes equations for three dimensional models [START_REF] Martin | Three-dimensional modelling of the dynamics of Johnsons Glacier, Livingston Island, Antarctica[END_REF][START_REF] Zwinger | A full Stokes-flow thermomechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka[END_REF]Pattyn, 2008;[START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF], but these models are difficult to use for continental scale simulations.

These new approximations are much more computationally intensive than the Shallow Ice Approximation and cannot be used for paleoclimate reconstruction. Unlike models whose initialization is based on long spin-up, these new models must rely on data assimilations to reproduce the present day conditions. This method was introduced in ice sheet modeling by [START_REF] Macayeal | The basal stress distribution of Ice Stream E, Antarctica, Inferred by Control Methods[END_REF][START_REF] Macayeal | A tutorial on the use of control methods in ice-sheet modeling[END_REF] to infer unknown parameters, like basal drag friction or ice viscosity, and correctly reproduce the measured velocity. This method has been extensively used for the shelfy-stream approximation [START_REF] Rommelaere | Large-scale rheology of the Ross Ice Shelf, Antarctica, computed by a control method[END_REF][START_REF] Vieli | Application of control methods for modelling the flow of Pine Island Glacier, West Antarctica[END_REF][START_REF] Joughin | Basal shear stress of the Ross ice streams from control method inversions[END_REF][START_REF] Larour | Modélisation numérique du comportement des banquises flottantes, validée par imagerie satellitaire[END_REF][START_REF] Khazendar | Larsen B Ice Shelf rheology preceding its disintegration inferred by a control method[END_REF] and recently applied to higher-order and full-Stokes models by [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF]; Jay-Allemand et al. [2011].

Several large initiatives are underway to develope models able to improve projections of sea level rise. Some of the most active members of the modeling glaciology community include the Parallel Ice Sheet Model (PISM) developed by the University of Alaska [Bueler et al., 2005;Bueler and Brown, 2009], the Community Ice Sheet Model (CISM), based on the Glimmer community model and that is part of the Community Earth System Model [START_REF] Rutt | The Glimmer community ice sheet model[END_REF], SICOPOLIS [Greve, 1997a], one of the only model that treat ice as a polythermal fluid, and Elmer [START_REF] Zwinger | A full Stokes-flow thermomechanical model for firn and ice applied to the Gorshkov crater glacier, Kamchatka[END_REF], a software based on finite elements.

A fully operational, three dimensional, thermomechanically coupled, evolutive ice sheet model with higher-order or full-Stokes dynamics is not yet available. However, new computational ressources and efficiency of numerical methods have progressed to the point where February 10, 2012 such a model is likely to be possible in the coming years, taking advantage of adaptive grids to obtain high spatial resolution where required [START_REF] Blatter | A short history of the thermomechanical theory and modelling of glaciers and ice sheets[END_REF].

Further improvements

Although many improvements were made in ice sheet modeling during the last decade, many more still need to be done to improve the accuracy and reliability of ice sheet projections.

A major challenge facing the models is the quality of available datasets. The new higherorder models involve greater computational costs and long spin-ups become prohibitive, so they need to rely on data assimilation for their initialization. The consistency between the different present-day datasets becomes crucial as they drive the response of ice sheet models. As pointed out by [START_REF] Rasmussen | Bed topography and mass-balance distribution of Columbia Glacier, Alaska, USA, determined from dequential serial-photography[END_REF], the inconsistency of the datasets limits the reliability of the projection models. Models that are initialized with data assimilation to reproduce the present-day conditions start by artificially redistributing the glacier mass, not as a realistic projection but to reconcile the inconsistencies when the datasets are not consistent. We showed [START_REF] Seroussi | Ice flux divergence anomalies on 79north Glacier, Greenland[END_REF] that at a high-resolution, ice sheet models are fundamentally limited by the inconsistencies between ice thicknesses that are measured at low spatial resolutions (several kilometers) along flight tracks and surface velocities that are derived from InSAR (Interferometric Synthetic Aperture Radar) at very high resolutions (50 to 300 meters). Although classical ice sheet models were not affected by this problem, because of (1) their low spatial resolution and ( 2) long spin-ups, the quality of current ice thicknesses makes new ice sheet models extremely difficult to use. New ice sheet models are required to improve estimates of the ice sheet contribution to sea level rise but can be affected by sparse ice thickness observations. Full-Stokes models are thought to be essential to model the grounding line area, transition between the grounded and floating ice [START_REF] Nowicki | Conditions for a steady ice sheet-ice shelf junction[END_REF][START_REF] Durand | Full Stokes modeling of marine ice sheets: influence of the grid size[END_REF][START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] as other models fail to reproduce the ice behavior in this critical area at which most changes are occurring. However using full-Stokes models at a continental scale is so computationally intensive that this is a very difficult goal to achieve. Furthermore, a very high spatial resolution is needed in these critical areas [START_REF] Nowicki | Conditions for a steady ice sheet-ice shelf junction[END_REF][START_REF] Durand | Full Stokes modeling of marine ice sheets: influence of the grid size[END_REF] making the cost of full-Stokes models even more prohibitive. First applications using a three dimensional full-Stokes model with grounding line evolution have been made recently but are limited to idealized geometries [START_REF] Favier | A three-dimensional full Stokes model of the grounding line dynamics: effect of a pinning point beneath the ice shelf[END_REF].

Overview of this thesis

Ice sheet modeling is at the frontier of many disciplines: glaciology, computational mechanics, applied mathematics and computer science but also to a smaller extent remote sensing, data processing and all the Earth system components that interact with the ice sheets. This manuscript reflects this diversity and includes elements of both glaciology and numerical modeling. The main goal of this thesis is to propose some improvements to the numerical modeling of ice sheet thermo-dynamics, using enhanced multi-model formulations.

10

February 10, 2012

OVERVIEW OF THIS THESIS

The second chapter is a brief state of the art of the physics and numerics of ice sheet systems. We describe the different parts of an ice sheet system and the notations used in this thesis. We then detail the conservation of mass, momentum and energy principles and apply them to the case of ice, modeled as a viscous and incompressible fluid. We present the boundary conditions of the system for each model and derive the general thermo-dynamic problem describing our system. We finally list several methods that have been developed to combine different physical or mechanical models describing the same system.

The third chapter introduces several common ice flow approximations, all derived from the full-Stokes problem. We rigorously establish the equations and boundary conditions of these mechanical models from the weak formulation of the full-Stokes problem. We then discuss the domain of validity for each of these approximations.

The fourth chapter details the numerical implementation strategies as well as the algorithms used in this work. We explain the choices made for the space and time discretizations, present the weak formulation and Galerkin approximation of both mechanical and thermal models, present the algorithms used to solved the non linearity of these two problems and introduce the iterative schemes used to perform steady-state and transient simulations.

We then describe the strategies chosen for the software architecture, the parallelization and partitioning of the domain and present some simple applications to square geometries commonly used to compare model results.

The following chapter focuses on the coupling of different mechanical models. We first detail the tiling method, a new method adapted from the Arlequin framework that is adapted to multi-model formulations that can be advantageously used to combine ice flow models, and present its main characteristics. We then apply this method to couple three of the most common ice flow models: the shelfy-stream approximation, a higher-order model and the full-Stokes equations. We emphasize the specificities of the coupling between these three models.

The last chapter proposes some applications and results on both idealized and real geometries, to validate the technique and its implementation. We first validate the method on square ice sheets and ice shelves, then on an ice sheet with a bumpy bed that is not transmitted to the surface (an interesting test case as full-Stokes and simpler models behave differently) and finally apply the method to the case of Pine Island Glacier, West Antarctica, as previous studies showed that a full-Stokes model is necessary in the critical region of the grounding line. This chapter introduces the thermo-mechanical equations governing the physics of ice sheet systems. The mechanical quantities (momentum, energy, mass, ...) are subject to fundamental conservation laws. In the context of continuum mechanics, conservation laws lead to a set of field equations that must be satisfied for every point of Ω and at all time. Equations governing ice flow are the general equations used for any continuum solid or fluid and are derived from conservation principles. Ice sheets interact with other components of the Earth system: ocean, atmosphere and land. A detailed description of the thermal and mechanical properties at these interfaces is therefore needed to describe ice flow. We conclude this chapter by presenting some coupling methods that allow combining different physical or mechanical models and present the requirements for coupling different ice flow models of varying order of complexity.

Ice sheet system

We detail here the different parts of an ice sheet system and provide some definitions. We then describe the same system from a modeling perspective and define the notations. very large masses of ice that cover areas over 50 000 km 2 . In some regions, the ice slides over the underlying bedrock and the ice velocity can reach several thousands of meters per year. These rivers of ice are called ice streams. If they terminate into the ocean and start to float, the floating part is called ice shelf. Marine terminated glaciers are ice sheets whose bedrock rests below sea level. No floating part develops for this type of glacier, so ice directly calves into the sea to form icebergs when they reach the ocean. The grounding line marks the transition between the grounded ice and the floating ice. The almost vertical ice cliff at the seaward end of an ice shelf or a marine terminated glacier is the ice front. Ice sheet systems interact with the ocean, the atmosphere and the bedrock upon which they lay.

Components of ice sheet systems

Three-dimensional geometry and notations

To model this complex system, we adopt the following representation (Figure 2.2) that shows a cross-section of the system. Let Ω be the domain we consider, which represents indistinctly ice sheets, ice streams or ice shelves.

Bedrock

Ice Ω The border of the domain ∂Ω is divided into several sections to represent the interfaces between ice and other media (ocean, atmosphere, ...). We use the following notations:

Air

Ocean Γ u Γ i Γ b Γ w Γ s
• Γ s is the ice/air interface, the upper surface of the domain Ω

• Γ b is the ice/bedrock interface, which is part of the lower surface of Ω

• Γ i is the ice front, the almost vertical cliff where a glacier thins and ends. At the ice front, ice can be in contact with ocean or air.

• Γ w is the ice/water interface below ice shelves, it does not include the ice front

• Γ u is the lateral part of the border, ∂Ω, that is not an ice front. It includes points that are neither on the ice/air interface, nor on the ice/ocean interface, nor on the ice/bedrock interface, nor on the ice front:

Γ u = ∂Ω \ (Γ s ∪ Γ w ∪ Γ b ∪ Γ i ) • s (x, y) is the upper surface elevation • b (x, y) is the lower surface elevation • H (x, y) is the ice thickness (H = s -b)

Plane view geometry and notations

Some approximations and equations are bidimensional, in the horizontal plane. The notations for these bidimensional problems are slightly different (see figure 2.3). Ice occupies the domain ω of the plan, which represents indistinctly ice sheets, ice streams and ice shelves.

γ i Ice ω γ u Bedrock Ocean Figure 2.3: Bidimensional modeled ice sheet system (top view)
The border of the domain ∂ω is divided into:

• γ i the ice front • γ u the part of the border ∂ω that is not on the ice front: γ u = ∂ω \ γ i Notations for the ice thickness, H, the upper surface elevation, s, and the lower surface elevation, b, remain the same. For a plane view problem, the ice/bedrock, ice/air and ice/ocean interfaces do not exist contrary to the three dimensional problem; these conditions are included in the description of the plane view models.

Mass balance 2.2.1 Mass balance equation

Law of conservation of mass

Conservation of mass leads to the general mass balance equation or continuity equation:

∂ρ ∂t + ∇ • ρv = ṡ (2.1)
where ρ is the mass density of the ice, v its velocity, ∇ the divergence operator and ṡ a function describing the production and destruction of ice (source and sink). This equation is independent of the coordinate system used.

Ice mass balance

We consider that the accumulation of ice only occurs at the surface of the glacier and that melting only occurs at the glacier's base, so there is no local production or destruction of ice. The accumulation and melting that take place on the upper and lower surface are not considered in the local mass balance equation:

∂ρ ∂t + ∇ • ρv = 0 (2.2)

Ice incompressibility

For an incompressible body, the density remains unchanged during the motion:

dρ dt = 0 (2.
3)

It follows that the mass balance equation (Eq. 2.2) is reduced to:

∇ • v = 0 (2.4)
The incompressibility of the ice is a reasonable assumption made in all ice sheet models [START_REF] Hooke | Principles of Glacier Mechanics[END_REF], p10). This is not true near the surface of a glacier where snow and firn are undergoing compaction, but it is a valid assumption in most of the ice mass. The depth at which firn becomes ice (density of 839 kg/m 3 ) varies depending on the accumulation rate and the temperature, but this depth is typically 60m to 70m [START_REF] Paterson | The Physics of Glaciers[END_REF], to be compared to the total ice thickness that reaches 3000m.

Boundary condition

We define here the kinematic boundary conditions for the glacier lower and upper surface.

Surface kinematic boundary condition

The evolution of the surface of a glacier is defined by:

∂s ∂t + u s ∂s ∂x + v s ∂s ∂y -w s = Ṁs (2.5)
where Ṁs is the surface mass balance in ice equivalent and u s , v s , w s the components of velocity at the surface.

Lower kinematic boundary condition

The evolution of the lower surface of a glacier is defined by:

∂b ∂t + u b ∂b ∂x + v b ∂b ∂y -w b = Ṁb (2.6)
where Ṁb is the basal mass balance in ice equivalent and u b , v b , w b the components of velocity on the lower surface.

Depth-integrated mass balance equation

Glacier evolution is dictated by geometrical changes. These geometrical changes are due to mass conservation, which includes dynamic thinning/thicknening and external forcings.

The conservation of mass states that the temporal change in ice thickness is the difference between the net mass balance of the glacier (i.e., surface mass balance plus basal mass balance) and the volume flux divergence.

To establish the mass conservation equation, we integrate the incompressibility equation from the bed to the surface:

s(x,y) b(x,y) ∂u ∂x + ∂v ∂y + ∂w ∂z dz = 0 (2.7)
We define the depth-averaged horizontal velocity v = (ū, v) T , such that:

ū = 1 H s(x,y) b(x,y) udz v = 1 H s(x,y) b(x,y) vdz (2.8) Equation (2.7) is equivalent to: s(x,y) b(x,y) ∂u ∂x + ∂v ∂y dz + w (x, y, s(x, y)) -w (x, y, b(x, y)) = 0 (2.9)
We write w s = w (x, y, s(x, y)) and w b = w (x, y, b(x, y)).

To integrate the two other terms in the integral, we use Leibniz rule of integration 1 which gives: Using the bed and surface kinematic boundary conditions and the definition of the depthaveraged velocity gives:

∂ ∂x s(x,
∂ ∂x (H ū) + ∂ ∂y (H v) + ∂s ∂t -Ṁs - ∂b ∂t + Ṁb = 0 (2.12)
As the ice thickness is H = s -b and introducing the bidimensional divergence operator ∇•, the depth-averaged mass conservation equation is:

∂H ∂t = -∇ • (vH) + Ṁs -Ṁb (2.13)
1 Leibniz rule of integration:

∂ ∂x b(x) a(x) f (x, z) dz = b(x) a(x) ∂ ∂x f (x, z) dz + ∂b ∂x f (x, b (x)) - ∂a ∂x f (x, a (x)) (2.10)

Momentum balance

Momentum balance equation

Law of conservation of Momentum

Newton's second law of motion states that the sum of all forces is equal to the rate-of-change of linear momentum. Ice is considered as a deformable continuum body so the equilibrium of linear momentum is the same as for any solid or fluid subjected to gravity (Coriolis effect is negligible):

∇ • σ + ρg = ρ dv dt (2.14)
where σ is the Cauchy stress tensor, ρ the ice density, v the ice velocity, dv/dt the material derivative and g the acceleration due to gravity.

Ice momentum balance

Although recent observations show an acceleration of many glaciers in Greenland and Antarctica [Rignot et al., 2002b;[START_REF] Joughin | Timing of recent accelerations of Pine Island Glacier, Antarctica[END_REF][START_REF] Rignot | Changes in the velocity structure of the Greenland ice sheet[END_REF], these accelerations are small and the terms associated to these accelerations are negligible compared to the acceleration due to gravity in particular [START_REF] Reist | Mathematical analysis and numerical simulation of the motion of a glacier[END_REF]. Ice is therefore generally modeled with quasi-static models [START_REF] Paterson | The Physics of Glaciers[END_REF], p258):

∇ • σ + ρg = 0 (2.15)

Ice behavior law

Incompressibility

Because ice is considered incompressible, its behavior law only involves the deviatoric stress tensor. The isotropic pressure does not contribute to its deformation [START_REF] Hooke | Principles of Glacier Mechanics[END_REF], p13 and [START_REF] Paterson | The Physics of Glaciers[END_REF]). The stress tensor is decomposed into a deviatoric stress and pressure as follows:

σ = σ - 1 3 T r (σ) [I] = σ + p [I] (2.16)
where σ is the deviatoric stress tensor, p the isotropic pressure and [I] the identity matrix.

Let (x, y, z) be a cartesian coordinate system, with z the vertical axis pointing upward.

In terms of deviatoric stress components and pressure, the momentum balance equations (2.14) (2.17)

Isotropy

Ice is often assumed to be an isotropic material, but this is not exact for thick ice masses [START_REF] Gow | Rheological implications of internal structure and crystal fabrics of west antarctic ice sheet as revealed by deep core drilling at byrd-station[END_REF][START_REF] Azuma | Formation processes of ice fabric pattern in ice sheets[END_REF]. Studies on Greenland and Antarctic ice cores [START_REF] Thorsteinsson | Textures and fabrics in theGRIP ice core[END_REF][START_REF] Azuma | Textures and fabrics in the Dome F (Antarctica) ice core[END_REF][START_REF] Bargmann | On a thermodynamically consistent flow model for induced anisotropy in polar ice[END_REF] show that ice is isotropic only near the upper surface, as the crystallographic axes are randomly distributed. Preferred directions start to develop deeper into the ice, from a random orientation distribution at the upper surface to a vertical single maximum fabric along the direction close to the vertical next to the glacier bed. Ice can only be considered isotropic in the first hundreds of meters before a girdle starts to develop (sometimes observed even in the lower part of the firn) and the fabric is completely anisotropic after a couple of thousands of meters. However, the primary purpose of these drillings is to understand the climatic changes in the past. Therefore drillings are performed in areas with very little or no movement. Ice properties in ice streams and other fast moving areas remain poorly understood.

Models that include a description of ice fabric exist. They have been used mainly to study ice divides and compare modeled ice anisotropy and fabrics to ice core data [Gagliardini and[START_REF] Gagliardini | Simulation of anisotropic ice flow and fabric evolution along the GRIP-GISP2 flowline, central Greenland[END_REF] and are often based on orthotropic linear flow laws [START_REF] Gillet-Chaulet | A user-friendly anisotropic flow law for ice-sheet modelling[END_REF]. We will here consider ice as an isotropic material for the sake of simplicity.

This assumption is common to most large scale ice flow models even if theoretical studies start quantifying the influence of anisotropy on the ice flow through the use of enhancement factors [START_REF] Ma | Enhancement factors for grounded ice and ice shelves inferred from an anisotropic ice-flow model[END_REF]].

Glen's flow law

As ice is considered viscous and incompressible, its behavior law only involves the deviatoric stress tensor:

σ = 2µ ε (2.18)
where µ is the viscosity. The viscosity is scalar as ice is considered isotropic.

The most common behavior law used to model ice is viscous power law whose coefficients are based on John W. Glen's experiments [START_REF] Glen | The creep of polycrystalline ice[END_REF]. This law is known as Glen's flow law in glaciology and stands that:

εe = σ e B n (2.19)
where B is the ice hardness, which increases as the ice becomes stiffer, and n is a parameter empirically determined (most studies have found that n 3, see [START_REF] Hooke | Principles of Glacier Mechanics[END_REF], p15 and Paterson [1994], p86). This parameter is called flow law exponent. εe and σ e are the second invariants of the strain rate and deviatoric stress tensors, defined as:

εe = 1 √ 2   i,j=1..3 ε2 ij   1/2 σ e = 1 √ 2   i,j=1..3 σ 2 ij   1/2
(2.20)

The effective viscosity can be identified as:

µ = B 2 εe n-1 n (2.21)

Viscosity temperature dependence

The viscosity parameter B depends on many parameters (ice temperature, fabric, etc). As it is mostly affected by temperature, it can be described with a simple Arrhenius relationship:

A = A 0 exp - Q RT h (2.22)
where A = B -1/n is the creep parameter, T h the Kelvin temperature adjusted for melting point depression, A 0 a prefactor, R the universal gas constant and Q the activation energy for creep.

Field measurements of glacier ice give a value around 60 kJ.mol -1 for the creep activation energy [START_REF] Paterson | The Physics of Glaciers[END_REF][START_REF] Cuffey | The Physics of Glaciers[END_REF]. Ideally the prefactor A 0 should be constant, observations show a variability that could be due to ice fabric or ice purity.

Above -10 • C, ice softens more than predicted by the Arrhenius relationship and tabulated values based on field measurements such as those provided in table 2.1 are often used [START_REF] Paterson | The Physics of Glaciers[END_REF][START_REF] Cuffey | The Physics of Glaciers[END_REF].

Mechanical boundary conditions

Ice/Air interface

The upper surface of ice sheets and glaciers is in contact with the atmosphere. This interface is considered as a free surface as the air pressure is negligible compared to ice pressure. The mechanical boundary condition at the ice/air interface is therefore:

Temperature ( Table 2.1: Recommended factor for the creep parameter A for a flow law exponent n = 3 [START_REF] Paterson | The Physics of Glaciers[END_REF] σn 0

(2.23)

Ice/Ocean interface

For ice shelves and glacier floating termini, the lower surface of ice is in contact with the ocean as well as some parts of the ice front. At the ice/ocean interface, a Neumann boundary condition is therefore applied.

The force applied is equal to the water pressure:

σn = -p w n = ρ w gb n (2.24)
where ρ w is the ice density and b the ice lower surface elevation with respect to sea level.

At the ice front, ice can either be in contact with air or in contact with both air and ocean (ocean on the lower part and air on the upper part). A Neumann boundary condition is therefore applied. If ice is only in contact with air at the ice front, the boundary condition is similar to the ice/air interface:

σn 0 (2.25)
Where ice is in contact with the ocean, the force applied is equal to the water pressure:

σn = ρ w gzn (2.26)
where z is the elevation of the point considered. This condition can therefore be summarized as: σn = -p i n = ρ w g min (0, z) n (2.27) so p i = p w for z < 0 (points below see level) but is equal to zero for areas above sea level.

Ice/Bedrock interface

At the ice/bedrock interface, we specify the normal velocity and the tangential force. A viscous friction law is used to model basal friction:

τ b = -α 2 u τ (2.28)
where τ b is the tangential components of external forces, u τ = u -(u • n) n the tangential velocity at the ice/bedrock interface and α a "friction-like" coefficient. This viscous friction law is commonly used in ice sheet modeling [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment -Theory and application to Ice Stream-B, Antarctica[END_REF][START_REF] Cuffey | The Physics of Glaciers[END_REF].

We ensure that no ice penetrates into the bedrock using a non-interpenetration condition at this ice/bedrock interface. Assuming that the bedrock does not move and there is no melting or freezing at the ice/bedrock interface, the velocity normal to the ice/bedrock interface is set to zero and equation (2.6) is reduced to:

w b = u b ∂b ∂x + v b ∂b ∂y (2.29)
where (u b , v b , w b ) are the components of velocity on the lower surface,

Full-Stokes equations

If we introduce the behavior law described above (section 2.3.2) and add the incompressibility equation in system (2.17), we find the Stokes equations, commonly referred to as full-Stokes equations in glaciology. Let (u, v, w) be the velocity components in a the Cartesian coordinate system, the equations read:

                                   ∂ ∂x

Simplified mechanical models

These three simplified models are all derived from the full-Stokes equations. In this section, we just present the equations, their rigorous derivation will be presented in chapter 3.

Higher-order model

The three dimensional higher-order model was derived by [START_REF] Blatter | Velocity and stress-fields in grounded glaciers: A simple algorithm for including deviatoric stress gradients[END_REF] and [START_REF] Pattyn | A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes[END_REF].

The assumptions made are:

                                   ∂ ∂x

Shelfy-stream approximation

The bidimensional shelfy-stream approximation was pioneered by [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment -Theory and application to Ice Stream-B, Antarctica[END_REF] and [START_REF] Morland | Plane and radial ice-shelf flow with prescribed temperature profile[END_REF]. The assumptions made in this model are:

• ∂w ∂x << ∂u ∂z • ∂w ∂y << ∂v ∂z • ˙ xz = 0 • ˙ yz = 0
So the full-Stokes equations, after depth-average integration are reduced to the following bidimensional equations:

                                   ∂ ∂x

Shallow ice approximation

The shallow-ice approximations is a semi-analytical three dimensional model presented by [START_REF] Hutter | Theoretical glaciology: material science of ice and the mechanics of glaciers and ice sheets[END_REF]. The only components of the stress tensor not neglected in this model are σ xz and σ yz . Additionally the horizontal derivative of vertical velocity are neglected compared to the vertical derivatives of horizontal velocities so the full-Stokes equations are reduced to:

                                 ∂ ∂z

Other models

Other models exist to describe the flow of ice sheets (see e.g. [START_REF] Hindmarsh | A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling[END_REF] for a more complete list of existing models). Combinations of different models are also possible. [START_REF] Pollard | Modelling West Antarctica ice sheet growth and collapse through the past five million years[END_REF] model the evolution of the West Antarctic during the last five million years with a combination of models. Their model uses the shallow ice approximation (SIA) for grounded ice [Hutter , 1982a] and the shallow-shelf approximation (SSA) for floating ice [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment -Theory and application to Ice Stream-B, Antarctica[END_REF]. As the SIA is an analytical model, it is not influenced by the results computed on the floating part with the SSA, the two models can therefore be computed separately. In order to capture the grounding line effects, a mass-flux condition is used at the transition between the two models [START_REF] Schoof | Ice sheet grounding line dynamics: Steady states, stability, and hysteresis[END_REF].

Another hybrid modeling approach, that combines SSA and SIA, is presented by Bueler and Brown [2009]. Both SSA and SIA are computed separately, and the ice velocity is taken as a weighted sum of the SSA and the non-sliding SIA. The SSA is used as a sliding law for the SIA, and a parameter is used to balance the amount of displacement due to SSA and SIA, ranging from SSA only on floating parts to SIA only when ice is frozen to the bedrock. This approach does not consist in model coupling, as the ice flow model is taken as linear combination of the two contributions, assumed to be independent.

As these models took the first steps in integrating hybrid models into the field of ice sheet modeling, here we propose to go further and introduce a strong coupling between the main ice flow approximations. A rigorous coupling of ice flow models is described in chapter 5.

Energy balance 2.4.1 Energy balance equation

Law of conservation of Energy

The energy balance equation states that for an isolated system, the total energy is constant over time. If we consider that the heat flux vector follows Fourier's law, the conservation of energy is:

d dt (ρcT ) = ∇ • (k ∇ T ) + Φ (2.34)
where T is the ice temperature, c is the ice heat capacity, k the ice heat conductivity and Φ the internal source of energy, such as viscous heating due to deformation.

The material derivative is the sum of the local derivative and the convective derivative:

d dt (ρcT ) = ∂ ∂t (ρcT ) + v • ∇ (ρcT ) (2.35)
So the conservation of energy imposes:

∂ ∂t (ρcT ) = -v • ∇ (ρcT ) + ∇ • (k ∇ T ) + Φ (2.36)
The deformational heating or viscous heating due to internal deformation of ice [START_REF] Paterson | The Physics of Glaciers[END_REF] is:

Φ = T r (σ ε) (2.37)

Ice thermal model

The temperature of the ice affects the flow through the parameter B, and also affects the melting rate, which in turn affects the geometry and sliding of the glacier. A thermal model is therefore essential to have a realistic model of ice flow. The temperature distribution is determined by solving the thermal equation that only include one unknown.

We use several assumptions in the thermal model. The heat capacity and conductivity of the ice are considered constant: we neglect their thermal and spatial dependence. These assumptions are common and described by [START_REF] Hooke | Principles of Glacier Mechanics[END_REF], p117, and lead to a linear heat equation:

∂T ∂t = -v • ∇T + k ρc ∆T + Φ ρc (2.38)
The local heat transfer is thus a result of advection, conduction, and internal deformation heating. The ice velocity is one parameter of the thermal equation. The thermal and mechanical problems cannot be solved independently and we therefore have a couple thermomechanical problem. Additionally, if ice temperature reaches the pressure melting point2 , ice starts to melt and the energy is used to transform ice into water. The temperature must therefore verify the assumption T < T pmp where T pmp is the pressure melting point. This constraint is discussed in section 2.4.3.

Thermal boundary conditions

Ice/Air interface

The temperature on the upper surface is considered to be equal to the air temperature, we therefore impose it to be equal to the mean annual temperature:

T s = T air (2.39)
Using the mean annual temperature as a boundary condition for the thermal model is a common way to impose surface tempature (i.e., [START_REF] Hulbe | A new numerical model of coupled inland ice sheet, ice stream, and ice shelf flow and its application to the West Antarctic Ice Sheet[END_REF]; [START_REF] Pattyn | A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes[END_REF]). If this approximation is not true in general, it remains a good assumption for cold and dry sites where the maximum air temperature rarely rises to 0 • C. As the temperature variability decreases with depth, the temperature of the air differs at most by 2 • C compared to the firn temperature at 10 meters, with an average difference around 0.7 • C (e.g., [START_REF] Cuffey | The Physics of Glaciers[END_REF] p. 405, [START_REF] Brandt | Temperature measurements and heat transfer in near-surface snow at the South Pole[END_REF]).

Ice/Ocean interface

At the ice/ocean interface, a heat flux is imposed. It mainly depends on the temperature difference between the ocean and the ice. Simple parameterizations can be used to specify this flux. [START_REF] Holland | Modeling thermodynamic ice-ocean interactions at the base of an ice shelf[END_REF] propose the following relation:

k ∇ T | b • n = -ρ w c pM γ (T b -T pmp ) (2.40)
where T b the ice temperature on the lower surface, T pmp the pressure metling point, k the ice thermal conductivity, n the outward pointing normal vector, ρ w the water density, c pM the mixed layer specific heat capacity and γ the thermal exchange velocity.

Ice/bedrock interface

At the ice/bedrock interface, heat flux comes from both the geothermal heat flux and heat due to basal friction. When ice temperature reaches the pressure melting point at the base, ice starts to melt and energy is dissipated with this melting. The thermal boundary condition is therefore:

k ∇ T | b • n = G -τ b • u τ (2.41)
where G is the geothermal heat flux and the u τ the tangential velocity at the ice/bedrock interface.

Thermal problem

There are two different ice thermal regimes:

1. Cold ice: all the ice is below the melting point and there is no liquid water 2. Temperate ice: ice is at pressure melting point everywhere; solid ice and liquid water coexist.

These two regimes sometimes coexist and ice is then defined as polythermal : some areas contain cold ice, others temperate ice.

The heat budget equation presented previously in this chapter does not account for phase changes and does not include the constraint that ice temperature cannot exceed pressure melting point and is not appropriate to model temperate ice. If the heat provided is sufficient to increase the ice temperature to the pressure melting point, the excess heat induces ice to melt into water. Several methods exist to tackle this problem of phase change. We briefly present here some of the techniques used to represent ice temperature.

Ice temperature in Greenland and Antarctica

Ice reaches melting point not only in mountain glaciers [Aschwanden and Blatter , 2009] but also in some parts of the Greenland ice sheet [Greve, 1997a, b] and the Antarctic ice sheet [START_REF] Calov | Simulation of the antarctic ice sheet with a three-dimensional polythermal ice-sheet model, in support of the epica project[END_REF][START_REF] Savvin | Simulation of the Antarctic ice sheet with a three-dimensional polythermal ice-sheet model, in support of the EPICA project. II. Nested high-resolution treatment of Dronning Maud Land, Antarctica[END_REF][START_REF] Siegert | Spectral roughness of subglacial topography and implications for former ice-sheet dynamics in East Antarctica[END_REF]. Greve [1997a] showed that the basal temperature reaches melting point in almost half of the Greenland ice sheet, these areas being mainly located along the coast. Aschwanden et al. [submitted] find similar results using an enthalpy method (see section 2.4.3). A comparison of basal temperature and melting rate is provided in figures 2.4 and 2.5 for a cold ice model (see section 2.4.3) and a model using an enthalpy formulation.

For the Antarctic ice sheet, [START_REF] Siegert | Spectral roughness of subglacial topography and implications for former ice-sheet dynamics in East Antarctica[END_REF] showed that the ice basal temperature is at the pressure melting point in many areas and that the basal water thickness reaches 0.2 m in some areas, where the subglacial temperatures are warm and the bedrock is smooth, which induced strong basal sliding. so on. We have identified which fields are needed as boundary conditions for an enthalpy formulation, possibly provided through such coupling.

The enthalpy formulation can be used to simulate both fully-cold and fully-temperate glaciers. Neither prior knowledge of the thermal structure nor a parameterization of the cold-temperate transition surface (CTS) is required.

Our application to Greenland shows that an enthalpy formulation can be used for continent-scale ice flow problems. The magnitude and distribution of basal melt rate, a thermodynamical quantity critical to modeling fast ice dynamics in a changing climate, can be expected to be more realistic in an energy-conserving enthalpy formulation than in most existing ice sheet models.

Cold ice model

Cold ice models treat ice as an incompressible viscous fluid that conducts heat and contains only one constituent: solid ice. To prevent ice temperatures from exceeding the pressure melting point, a Dirichlet condition is generally employed to impose the temperature to be at the pressure melting point at these locations [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF]:

T pmp (p) = T 0 -c (p -p 0 ) (2.42)
where T pmp is the ice pressure melting point at pressure p, T 0 the reference pressure melting point (0.01 • C) at pressure p 0 (630 Pa).

The excess heat is artificially transformed into basal mass balance, for example in Rutt et al.

[2009]:

Ṁb = k ρL ∂T * ∂z - ∂T ∂z (2.43)
where Ṁb is the basal mass balance (melting rate), k the ice thermal conductivity, ρ the ice density, L the specific latent heat of fusion, T * the temperature computed (before imposing the temperature to be equal to the pressure melting point) and T the final temperature (equal to the pressure melting point).

Greve [1997b] and Aschwanden et al. [submitted] showed that using this model instead of a polythermal model creates only small differences on the Greenland ice sheet. Greve [1997b] the ice immediately above is cold. Cold ice models do not account for the latent heat stored as liquid water within temperate ice, such schemes are therefore not conserving energy when temperate ice is present. However, a cold ice model is a reasonable and simple way of representing ice temperature in the ice, as long as the pressure melting point is only reached on basal points.

Temperate ice model

Temperate ice is at the pressure melting point everywhere and ice is composed of both solid ice and liquid water. Temperate ice models are based on the classical mixture theories and include a moisture content [START_REF] Müller | Thermodynamics[END_REF]. Each point is considered to be simultaneously occupied by all the constituents of the mixture, here ice and water.

There is a hierarchy of mixture models [START_REF] Müller | Thermodynamics[END_REF] depending on the amount of the components into the mixture. The first category is appropriate for modeling diffusive motion of a particular substance contained as a contaminant of impurity in another substance. In this case, the balance laws of momentum and energy are only formulated for the mixture as a whole. The second category models the interaction forces between the constituents. The energy balance is formulated for the mixture as a whole but both momentum and mass balances are formulated for all the constituents. The last category is a full-mixture thermodynamics where all balance equations are formulated for all the components.

As only a small fraction (less than 5%) of the volume is composed of water [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF]; Greve [1997a], the water can be considered as a tracer component whose motion is relative to the mixture model and described by Fickian diffusion [START_REF] Müller | Thermodynamics[END_REF]. In this model, only one momentum and energy balance equations are formulated for the mixture as a whole but two mass balances are used, one for each component of the mixture [START_REF] Müller | Thermodynamics[END_REF];

Hutter [1993].

Müller [1985] introduces the volume fraction of water in the mixture ω v and the mass fraction of water in the mixture ω m :

ω v = V w V and ω m = m w m (2.44)
where V w is the volume of water, V the total volume of the mixture, m w the mass of water and m the total mass of the mixture. If we write ρ i the ice density and ρ w the water density, the density of the mixture ρ and its velocity v are:

ρ = (1 -ω v ) ρ i + ω v ρ w ρv = (1 -ω v ) ρ i v i + ω v ρ w v w (2.45)
with v i the ice velocity and v w the water velocity.

As the amount of water into the ice is small, the density of temperate ice stays close to the cold ice density (variation of less than 0.5% if there is 5% water in the mixture [Greve, 1997a]), so the mixture is considered incompressible and its mass balance is

∇ • v = 0.
The energy balance of the system is:

∂ ∂t = -∇ • ( v + D) + Φ (2.46)
where is the internal energy of the system, D the diffusivity of the mixture and Φ the internal deformation heat.

As the ice is temperate, all constituents of the mixture are at pressure melting point. The energy of the mixture is therefore:

= c T pmp + L ω m (2.47)
where c is the ice heat capacity and L the latent heat of fusion.

If we consider the diffusivity of the mixture to follow a Fourier law D = -k∇T pmp (because the mixture temperature is the pressure melting point), the thermal diffusivity of temperate ice to be k and the motion of liquid water relative to the mixture motion to be described by a Fickian diffusion law ω m (v w -v) = -ν∇ω m , the energy balance of the mixture is after transformations:

d dt (cT pmp + Lω m ) = Lν∆ω m + ∇ • (k∇T pmp ) + Φ (2.48)

Polythermal ice model

Ice is called polythermal when it contains both cold ice and temperate ice. Polythermal ice models include both types of ice, these two areas being separated by the cold-temperate transition surface (CTS) [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF]. Polythermal ice models have been developed by [START_REF] Fowler | Flow of polythermal glaciers. 1. Model and preliminary-analysis[END_REF]; [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF] and Greve [1997a] to correctly include phase changes.

The main addition of these models is therefore the inclusion of the cold-temperate transition surface to separate the cold ice from the temperate ice. The transition between cold and temperate ice is not smooth but constitutes a phase-change interface and therefore happens at a singular interface where Stefan-type matching conditions are applied.

Polythermal models track liquid water produced in temperate ice areas by dissipation heating. These models therefore have the ability to better compute the ice temperature as well as basal melt rates, which matters in fast flow areas that are controlled by the presence of pressurized water in the ice bed [START_REF] Schoof | Ice-sheet acceleration driven by melt supply variability[END_REF]. This model has been successfully applied, in the case of Greenland in particular [Greve, 1997a, b].

However, energy jump conditions at the cold-temperate surface transition are difficult to apply and make polythermal models very complicated to develop and implement. An alternative solution is to employ enthalpy formulations.

Enthalpy method

An alternative to polythermal models as described in the previous section is the enthalpy method. Although frequently used in fluid dynamics, this method is rarely employed in ice sheet modeling. Mentioned by [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF], a simplified version was used in a flow line model by Calvo et al. [1999]. [START_REF] Aschwanden | Mechanics and thermodynamics of polythermal glaciers[END_REF] derived a mathematical model based on enthalpy before proposing a general enthalpy formulation for glaciology [Aschwanden et al., submitted].

An enthalpy formulation can be used for both cold and temperate ice: a small change in enthalpy reflects a small change in temperature for cold ice and a small change in the fraction of liquid water for temperate ice. Both cold and temperate ice are therefore treated in a single formulation with an enthalpy equation. This method does not require tracking the cold-temperate ice surface and is consequently computationally advantageous. Furthermore, the energy fluxes in supraglacial runoff and subglacial aquifers can be treated in the same framework. This framework unifies the treatment of the energy conservation for en-glacial, sub-glacial and supra-glacial liquid water [Aschwanden et al., submitted].

The specific enthalpy H and internal energy are linked as follows [Aschwanden and Blatter , 2009]:

H = + p ρ (2.49)
where is the specific internal energy and p the pressure. If a material is heated under constant pressure (isobaric), the specific enthalpy is directly dependent on the specific internal energy: ∆H = ∆ . We only consider this case here.

For cold ice, we suppose that the diffusivity follows a Fourier law with a diffusivity coefficient k i . As a change of enthalpy is only associated to a change in temperature, the enthalpy equation for cold ice is:

ρ ∂H ∂t + v • ∇H = ∇ • k i C i ∇H + Q (2.50)
In temperate ice, the ice temperature is at the pressure melting point and a change of enthalpy is associated to a change in the water fraction. The thermal conductivity of the mixture is

k = (1 -ω m ) k i + ω m k w .
As for the case of temperate ice, a Fick-type diffusion is used to express the diffusive water flux [Hutter , 1982b]:

ρ w ω v L (v w -v) = - k 0 L ∇H (2.51)
where k 0 is a small positive constant.

The enthalpy equation for temperate ice is:

ρ ∂H ∂t + ρv • ∇H = ∇ • k 0 L ∇H + ∇ • (k∇T pmp ) + Φ (2.52)
An enthalpy equation valid for both cold and temperate ice can be written as:

ρ ∂H ∂t + v • ∇H = Φ +        ∇ • k i C i ∇H if H < H s ∇ • k∇T pmp + k 0 L ∇H if H s < H < H l (2.53)
where H s is the enthalpy of pure ice and H l the enthalpy of pure liquid water at the pressure melting point T pmp .

This approach is appropriate for most polythermal glaciers except for fully temperate glacier where a momentum balance equation should be specified for each constituent [START_REF] Hutter | Theoretical glaciology: material science of ice and the mechanics of glaciers and ice sheets[END_REF]Aschwanden et al., submitted], as we make the assumption that liquid water is present in relatively small quantities (less than 5%) in temperate ice. In the case of Antarctica and Greenland were most of the ice is cold except for a temperate layer at the ice/bedrock interface, this model is appropriate.

In glaciers and ice sheets, enthalpy is a function of temperature, water content and pressure. The temperature and water fraction are retrieved from the enthalpy (see figure 2.6) as follows:

• For cold ice H < H s :

T = H -H s c + T pmp ω = 0 • For temperate ice H s < H < H l :
February 10, 2012

T = T pmp ω = H -H s L
where H s is the enthalpy of pure ice and H l the enthalpy of pure liquid water at the pressure melting point T pmp . The case of pure water where enthalpy is above the enthalpy of pure liquid water H l is not considered here.

Water fraction ω Temperature T

Enthalpy H 1 0 T pmp H s H l

Cold ice

Liquid water Temperate ice

Methodology adopted in this thesis

The work presented in this thesis focuses primarily on the modeling of the Greenland and Antarctic Ice Sheet. For these continents, ice mainly reaches the pressure melting point close to the bedrock. We decided to use a cold ice model for our work. However, as ice reaches the melting point mostly at the bedrock, we consider that ice is cold everywhere except at the ice/bedrock interface. The boundary condition on the lower surface must be changed to account for ice melting in areas where ice reaches the pressure melting point.

The treatment of this condition is detailed in a forthcoming chapter, section 4.3.

General ice flow problem

The general thermo-dynamic problem describing ice flow is modeled by the following local equations:

∀M ∈ Ω, ∀t ∈ [0, t f ] , Find (u, v, w, p, T ) :                                                          ∂ ∂x 2µ ∂u ∂x + ∂ ∂y µ ∂u ∂y + µ ∂v ∂x + ∂ ∂z µ ∂u ∂z + µ ∂w ∂x - ∂p ∂x = 0 ∂ ∂x µ ∂u ∂y + µ ∂v ∂x + ∂ ∂y 2µ ∂v ∂y + ∂ ∂z µ ∂v ∂z + µ ∂w ∂y - ∂p ∂y = 0 ∂ ∂x µ ∂u ∂z + µ ∂w ∂x + ∂ ∂y µ ∂v ∂z + µ ∂w ∂y + ∂ ∂z 2µ ∂w ∂z - ∂p ∂z -ρg = 0 ∂u ∂x + ∂v ∂y + ∂w ∂z = 0 ∂T ∂t = -v • ∇T + k ρc ∆T + Φ ρc ∂H ∂t = -∇ • (vT ) + Ṁs -Ṁb (2.54)
The problem is completed by the mechanical boundary conditions:

                                       Upper surface Γ s : σn = 0 Ice/Water interface Γ w : σn = p w n Ice Front Γ i : σn = p i n Ice sheet base Γ b : τ b = -α 2 u τ w b = u b ∂b ∂x + v b ∂b ∂y
Other lateral borders Γ u : u = u obs (2.55) and the thermal boundary conditions:

                   Upper surface Γ s : T = T air Icesheet base Γ b : k∇ T | b • n = G -τ b • u τ Iceshelf base Γ w : k ∇ T | b • n = -ρ w c pM γ (T b -T pmp ) Other borders Γ u ∪ Γ i : ∇ T • n = 0 (2.56)
This problem is described using the full-Stokes equations and a cold-ice model, but these components can be replaced by other mechanical and/or thermal models.

The mechanical and thermal problems are coupled. However the mechanical and thermal parts of the problem are generally solved independently. Details on the different mechanical models employed in ice sheet modeling are provided in chapter 3. The procedure used to combine thermal, mechanical and mass balance problems and compute steady-state and evolutive models is then described in chapter 4.

Multi-model methods

Several sets of equations with different degrees of sophistication are often available to model the physics of the same system; this is the case for ice sheet flow, as several approximations of the full-Stokes equations have been developed to reduce computational costs (see chapter 3). Developing new techniques for coupling models and scales of complex mechanical or physical problems is an active research area in the fields of computational mechanics and physics.

In this section we briefly list some common methods that have been developed to couple different models on a partition of the domain occupied by a physical system. We particularly emphasize the Schwarz method and the Arlequin framework [START_REF] Ben Dhia | Problèmes mécaniques multi-échelles : la méthode Arlequin (Multiscale mechanical problems: the Arlequin method)[END_REF], 1999] (see also e.g. Ben Dhia and [START_REF] Ben Dhia | Mathematical analysis of the mixed arlequin method[END_REF]; [START_REF] Rateau | Méthode arlequin pour les problèmes mécaniques multi-échelles, applications a des problèmes de jonction et de fissuration de structures élancées[END_REF]; Ben Dhia and [START_REF] Ben Dhia | The arlequin method as a flexible engineering design tool[END_REF];

Ben Dhia [2008]) since the method developed in this thesis is essentially related to these two approaches.

Problem definition

To describe the coupling methods, we will utilize the framework of continuum mechanics variational formulations. We consider as an example a domain Ω occupied by a material body. This domain is composed of two subdomains Ω 1 and Ω 2 that cover all the domain Ω.

We want to use two different physical or mechanical models on these two subdomains. The partitioning of the domain Ω is such that Ω 1 and Ω 2 overlap over a superposition zone or do not overlap at all (see e.g. figure 2.7).

Domain decomposition without superposition of subdomains is employed by most coupling techniques. In continuum mechanics, coupling of different models in such cases is based on the transmission of surface conditions at the interface. These transmissions are either exact or approximated.

The penalty method is an approximate method commonly used for its ease of use (see e.g., [START_REF] Courant | Variational methods for the solution of problems of equilibrium and vibrations[END_REF]). A penalty term, which can be interpreted as a spring with a very large stiffness, is added to the energy so that when a constraint is not fulfilled, the energy increases drastically. This penalized energy enforces the constraint in an exact manner when the spring stiffness goes to infinity. However the practical choice of a non-infinite parameter is difficult to parameterize. Moreover, the resulting numerical systems are ill-conditioned and therefore decrease the precision of the results.

Lagrangian methods are used to transmit the exact surface conditions. These approches are the basic elements for methods of domain decomposition, either primal or dual, and [START_REF] Farhat | A Method of Finite-element Tearing and Interconnecting and its Parallel Solution Algorithm[END_REF] introduced the FETI (Finite Element Tearing and Interconnecting) method that is based on a domain decomposition in which the subdomains are completely disconnected and Lagrange multipliers are introduced to ensure the compatibility at the interface between the subdomains (see also e.g. Farhat et al. [1994b, a]).

Ω 1 Ω 2
The LATIN (LArge Time INcrement) method was pioneered by [START_REF] Ladeveze | Sur une famille d'algorithmes en méchanique des structures[END_REF] for nonlinear time dependent problems and is based on a non-incremental iterative approach and further developed by Ladeveze and colleagues (see e.g. [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]; [START_REF] Dureisseix | A multi-level and mixed domain decomposition approach for structural analysis[END_REF]; [START_REF] Ladeveze | A micro/macro approach for parallel computing of heterogeneous structures[END_REF]). The algorithm provides an approximation of the solution for the entire domain and the entire time interval at each iteration.

Schwarz methods

This type of coupling technique was directly inspired by a purely mathematical work initiated by H. Schwarz [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF]. All the techniques derived from this work are referred to as Schwarz methods. These methods classically consider the union of two or more overlapping subdomains (see figure 2.8) and solve the equations of the two subdomains independently, using the latest values of the approximate solution as the boundary conditions.

These methods are either additive or multiplicative [START_REF] Xu | Iterative methods by space decomposition and subspace correction[END_REF][START_REF] Griebel | On The Abstract Theory Of Additive And Multiplicative Schwarz Algorithms[END_REF] and have been extensively used for parallel computations. Chimera methods [START_REF] Steger | On the use of composite grid schemes in computational aerodynamics[END_REF] are adapted from the same principles to deal with problems that are based

Ω 1 Ω 2 Γ 1 Γ 2 Figure 2
.8: Domain decomposition with overlapping subdomains on hybrid meshes. Recently, this methodology has begun to be used to handle multiscale problems [START_REF] Lions | Domain decomposition methods for CAD[END_REF][START_REF] Glowinski | Approximation of multi-scale elliptic problems using patches of finite elements[END_REF].

All Schwarz methods derive from the same basic idea described by H. Schwarz: an iterative algorithm is used to compute the unknown field on the two subdomains simultaneously assuming that Dirichlet conditions are known at the boundary of these subdomains; the Dirichlet conditions are then updated with the computed fields until we reach a convergence criterion. We write respectively Γ 1 and Γ 2 the internal boundaries of the two subdomains Ω 1 and Ω 2 (see figure 2.8). With these notations and for a general linear elastic problem where we are interested in finding the displacements u 1 and u 2 on Ω 1 and Ω 2 respectively, this algorithm is:

• Initialization of u 0 1 and u 0 2 • Compute u i 1 using u i-1 2 as a Dirichlet condition on Γ 1 • Compute u i 2 using u i-1 1 as a Dirichlet condition on Γ 2 • Iterate until both u i 1 and u i 2 converge
Many versions of this method exist, so we suggest that interested readers consult the literature [START_REF] Lions | On the Schwarz alternating method[END_REF][START_REF] Dolean | Optimized interface conditions for domain decomposition methods in fluid dynamics[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Dubois | Optimized Schwarz Methods for the Advection-Diffusion Equation and for Problems with Discontinuous Coefficients[END_REF].

Arlequin method

Finally, we describe the multiscale and multimodel Arlequin framework that has been pioneered by H. Ben Dhia [START_REF] Ben Dhia | Problèmes mécaniques multi-échelles : la méthode Arlequin (Multiscale mechanical problems: the Arlequin method)[END_REF], 1999] and further developed by H. Ben Dhia and colleagues (see e.g. Ben Dhia and [START_REF] Ben Dhia | Mathematical analysis of the mixed arlequin method[END_REF][START_REF] Ben Dhia | The arlequin method as a flexible engineering design tool[END_REF]; Ben Dhia [2008]; Ben Dhia et al.

[2008]), and also by other researchers in the field of computational mechanics (see e.g. [START_REF] Hu | Multi-scale modelling of sandwich structures using the Arlequin method Part I: Linear modelling[END_REF][START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF]; [START_REF] Bauman | Adaptive multiscale modeling of polymeric materials with Arlequin coupling and Goals algorithms[END_REF]). This approach is based on the superposition and the linking of different numerical models. Its main properties are to combine a superposition of mechanical models, a distribution of the energy between these models and the use of weak junctions to link them. Since this strategy has been used to derive the simplified method developed in this thesis (see chapter 5), we give here some details about this technique.

We consider here a standard mechanical problem of linear elasticity and we focus on the mathematical aspect of the Arlequin Framework. The theoretical justifications of this method are detailed in Ben Dhia and [START_REF] Ben Dhia | Mathematical analysis of the mixed arlequin method[END_REF]; [START_REF] Rateau | Méthode arlequin pour les problèmes mécaniques multi-échelles, applications a des problèmes de jonction et de fissuration de structures élancées[END_REF]; Ben Dhia and [START_REF] Ben Dhia | The arlequin method as a flexible engineering design tool[END_REF]; Ben Dhia [2008]. Interested readers are urged to consult these references as this section only presents the main ideas of the method.

Let Ω 1 and Ω 2 be two distinct bounded regular domains of R d , with d = 2 or 3 in practice.

Let Ω be the union of these two domains, and Ω S their intersection, assumed to be non zero measured and regular.

model superpostion

Ω 2 Ω 1 Ω S Figure 2
.9: Arlequin domains superposition Let V be the kinematically admissible velocity field and Γ u the part of the domain border where Dirichelt conditions are applied (for the sake of simplicity, we consider homogeneous Dirichet conditions). The linear elasticity problem reads:

Find u ∈ V = {u ∈ H 1 (Ω), u = 0 on Γ u }, ∀v ∈ V, a (u, v) = l (v) (2.57) with: ∀ (u, v) ∈ V × V a (u, v) = Ω σ (u) : ε (v) dΩ ∀v ∈ V l (v) = Ω f .vdΩ (2.58)
The dual volume-based coupling Arlequin formulation of the classical problem (2.57) reads:

February 10, 2012

Find (u 1 , u 2 , Φ d ) ∈ V 1 × V 2 × M d , ∀ (v 1 , v 2 , Ψ d ) ∈ V 1 × V 2 × M d : a (u 1 , v 1 ) + C d (Φ d , v 1 ) = l 1 (v 1 ) a (u 2 , v 2 ) + C d (Φ d , v 2 ) = l 2 (v 2 ) C d (Ψ d , u 1 -u 2 ) = 0 (2.59)
where it is assumed that the domain Ω is partioned into two overlapping subdomains Ω 1 and Ω 2 ; Ω S being the superposition and gluing zone (see figure 2.9). The spaces of kinematically admissible fields in each sub-domain are denoted by V 1 and V 2 . The weighted internal and external virtual works are defined as follows:

∀ (u 1 , v 1 ) ∈ V 1 × V 1 a 1 (u 1 , v 1 ) = Ω 1 α 1 σ (u 1 ) : ε (v 1 ) dΩ 1 ∀ (u 2 , v 2 ) ∈ V 2 × V 2 a 2 (u 2 , v 2 ) = Ω 2 α 2 σ (u 2 ) : ε (v 2 ) dΩ 2 ∀v 1 ∈ V 1 l (v 1 ) = Ω 1 β 1 f .v 1 dΩ 1 ∀v 2 ∈ V 2 l (v 2 ) = Ω 2 β 2 f .v 2 dΩ 2 (2.60)
where weight parameter functions (α 1 , α 2 ) and (β 1 , β 2 ) are defined respectively on Ω 1 and Ω 2 and satisfy (see e.g. figure 2.10, a partition ensuring the partition of energy):

α i ≥ 0 in Ω i α i = 1 in Ω i \ Ω S α 1 + α 2 = 1 in Ω S (2.61) β i ≥ 0 in Ω i β i = 1 in Ω i \ Ω S β 1 + β 2 = 1 in Ω S (2.62) M d is the dual mediator space. It is the dual space of M = V 1 | Ω S = V 2 | Ω S = H 1 (Ω S
) and the volume coupling operator, denoted by C d is [START_REF] Ben Dhia | Numerical modelling of multiscale problems: the Arlequin method[END_REF], 2008]:

∀ (Ψ, v) ∈ M d × M, C d (Ψ, v) = Ψ, v M d ,M (2.63)
where ., . M d ,M stands for the duality bracket. The dual volume coupling is a natural mechanical coupling operator, in the sense that, interpreting the Lagrange multiplier Ψ d as a density of forces, it has to be in the dual space M d of the space of the displacements in Ω S . We refer here to the references mentioned above for other variants of couplig or gluing operators.

0 1 α 2 α 1 Ω 2 Ω 1 Ω S Figure 2
.10: weight parameter functions in 1d

Requirements for coupling ice flow models

We mentioned in the previous section several techniques to couple different mechanical or physical models.

The penalty method is easy to implement, however the appropriate choice of the penalty parameter and the conditioning of the stiffness matrix are the main drawbacks. Schwarz algorithms could be used for the coupling of different mechanical or physical models (and for ice models, in particular). However, the convergence of the scheme could be a severe issue. Finally, the Arlequin framework is very general and allows coupling of any kind of models at different scales, as mentioned above. This framework (which has inspired other methods) seems to be mandatory for the coupling of complex mechanical systems (such as atomistic and continuum). But as it is a very general method, its implementation requires a lot of modifications in the structure of an existing code. Let us however observe that it has been very recently implemented within the Abaqus code [START_REF] Qiao | Implementation of the Arlequin method into ABAQUS: Basic formulations and applications[END_REF].

For these reasons, we present in chapter 5 a new method derived from the Arlequin framework and adapted to coupling ice flow models as it is:

• less general than Arlequin framework but sufficient for ice flow modeling

• able to couple models with different degrees of freedom

• easy to implement

February 10, 2012

Chapter summary

In this chapter we presented a brief overview of the physics and numerics of ice sheet models. We described the ice sheet systems and their interactions with other components of the Earth system. We then focused on the main three aspects of the modeling: mass conservation, momentum conservation and energy conservation. We established the local equations and detailed the boundary conditions applied for these three quantities, and derived the general system of equations that govern ice sheet thermo-dynamics.

Several physical or mechanical models generally exist to describe the same system and can be advantageously combined. In a last section of this chapter, we listed a few methods that exist to couple such models. The next chapter presents several ice flow models of varying order of complexity all derived from the full-Stokes equations. Solving the full-Stokes equations over the entire ice sheet is a numerical challenge especially for long paleoclimatic runs [START_REF] Deconto | Rapid Cenozoic glaciation of Antarctica induced by declining atmospheric CO2[END_REF][START_REF] Pollard | Modelling West Antarctica ice sheet growth and collapse through the past five million years[END_REF]. Several approximations of the momentum balance have been derived from the full-Stokes equations, taking advantage of the shallow aspect of the ice sheets and glaciers. The horizontal dimensions of these ice masses is generally several orders of magnitude larger than its thickness. The approximations consist in neglecting terms from the momentum balance equations and simplifying the strain rate tensor. We describe here the assumptions used to derive three of the most common approximations: the higher-order model, the shallow-shelf approximation or shelfy stream approximation and the shallow ice approximation. We then summarize the domain of validity for each of these approximations.

The weak formulation, based on the principle of virtual displacements or principle of virtual work, is a more general and more rigorous formulation of a given continuum mechanics problem compared to other formulations like the local equations (if singularities are present for example). For this reason, the purpose of this chapter is to derive simplified models of the full-Stokes equations using the weak formulation instead of deriving them from the local equations and boundary conditions as it is usually done [START_REF] Hutter | Theoretical glaciology: material science of ice and the mechanics of glaciers and ice sheets[END_REF][START_REF] Blatter | Velocity and stress-fields in grounded glaciers: A simple algorithm for including deviatoric stress gradients[END_REF][START_REF] Pattyn | A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes[END_REF]. In this chapter we use this approach to derive the higher-order model, the shelfy-stream approximation and the shallow ice approximation from the full-Stokes weak formulation. Previous work was done to derive the simplified models from the full-Stokes equations, based on the Hamilton-type principles or the principle of least action [START_REF] Bassis | Hamilton-type principles applied to ice-sheet dynamics: new approximations for large-scale ice-sheet flow[END_REF][START_REF] Dukowicz | Consistent approximations and boundary conditions for ice-sheet dynamics from a principle of least action[END_REF][START_REF] Dukowicz | Incorporating arbitrary basal topography in the variational formulation of ice-sheet models[END_REF] and this is to our knowledge the first time simplified models are derived using a weak formulation.

Full-Stokes variational equations

We use the local equations (2.30) and boundary conditions (2.55) defined in the previous chapter to establish the variational formulation of full-Stokes equations. We use the classical notations for Hilbert spaces: H 1 (Ω) is the Sobolev space that contains the square integrable functions whose first derivative is also square integrable over the domain, Ω. L 2 (Ω) is the Lebesgue space of square integrable functions over the domain, Ω. Let V u obs be the space of kinematically admissible velocity fields, P the space of admissible pressure fields and x ∈ Ω a point of coordinates (x, y, z):

V u obs = Φ ∈ H 1 (Ω) 3 , Φ = u obs on Γ u and Φ • n = 0 on Γ b P = q ∈ L 2 (Ω) (3.1)
Observe here that in order to simplify the presentation, without loss of generality, we will consider the case u obs = 0 and denote by V the space of velocity solution and kinematically admissible velocity fields. In practive, observed velocities are used to constrain the velocity on Γ u ; we use the classical lifting of Dirichlet boundary conditions.

For any kinematically admissible field (Φ, q) T such that Φ = (φ x , φ y , φ z ) T ∈ V and q ∈ P, the weak formulation of the full-

+ Γ b α 2 u τ • Φ τ dΓ = Ω ρg • ΦdΩ + Γw p w Φ • ndΓ + Γ i p i Φ • ndΓ (3.2)
where α is a friction coefficient, p i the pressure at the ice front, p w the water pressure applied at the ice/ocean interface Γ w , u τ the velocity at the ice/bedrock interface, tangential to the surface normal whose components on a cartesian coordonate system are u τ , v τ and w τ and Φ τ the part of the kinematically admissible field tangential to the ice/bedrock interface Γ b .

To simplify the equations, we apply the weak formulation to four test functions Φ x = (φ x , 0, 0, 0), Φ y = (0, φ y , 0, 0), Φ z = (0, 0, φ z , 0) and q = (0, 0, 0, q). This will simplify the derivation of the simpler models local equations and boundary conditions. With these test functions, the full-

+ Γ b α 2 w τ φ z dΓ = Γw p w n z φ z dΓ - Ω ρgφ z dΩ Ω - ∂u ∂x - ∂v ∂y - ∂w ∂z qdΩ = 0 (3.3)
February 10, 2012

3.2 Higher-order model

Model description

A three dimensional model known as higher-order model was derived from these equations by [START_REF] Blatter | Velocity and stress-fields in grounded glaciers: A simple algorithm for including deviatoric stress gradients[END_REF] and [START_REF] Pattyn | A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes[END_REF]. This model corresponds to an order one approximation of the full-Stokes equations with respect to the aspect ratio H/L where H is the glacier's height and L its length. This is why it is referred to as Higher-order. Two assumptions are necessary to derive this model. The first one is that the gradients of the vertical velocity are small compared to the vertical gradients of the horizontal velocity. The second one is that the bridging effect is negligible [van der Veen and [START_REF] Van Der Veen | Force budget: I. theory and numerical methods[END_REF]. We can summarize these assumptions as:

• ∂w ∂x << ∂u ∂z • ∂w ∂y << ∂v ∂z • ∂σ xz ∂x << ∂σ zz ∂z • ∂σ yz ∂y << ∂σ zz ∂z
An additional assumption on the friction law used will be necessary to fully decouple the system between horizontal and vertical equations.

Determination of variational and local equations associated to the higher-order model

Variational formulation

Let us now use these assumptions in the weak formulation of the full-Stokes equations (3.3) to establish the local equations of this model.

If we include the first two assumptions of the higher-order model (∂w/∂x << ∂u/∂z and ∂w/∂y << ∂v/∂z), for any kinematically admissible field (Φ, q) T such as:

Φ ∈ V q ∈ P (3.
+ Γ b α 2 v τ φ y dΓ = Γw p w n y φ y dΓ + Γ i p i n y φ y dΓ (3.5)
These equations couple the horizontal velocity fields and the ice pressure.

Considering admissible regular fields Φ z = (0, 0, Φ z , 0) with Φ z = 0 on the boundary of Ω and admissible fields q, and after integracting by parts the third system of the global full-Stokes system (3.3), we obtain:

Ω ∂ ∂x µ ∂u ∂z + µ ∂w ∂x φ z + ∂ ∂y µ ∂v ∂z + µ ∂w ∂y φ z + ∂ ∂z 2µ ∂w ∂z φ z - ∂p ∂z φ z dΩ = Ω ρgφ z dΩ Ω - ∂u ∂x - ∂v ∂y - ∂w ∂z qdΩ = 0 (3.6)
We now include the last assumption in the last equation (∂σ xz /∂x << ∂σ zz /∂z and ∂σ yz /∂y << ∂σ zz /∂z). These two assumptions are equivalent to:

• ∂ ∂x µ ∂u ∂z + µ ∂w ∂x << ∂ ∂z 2µ ∂w ∂z -p • ∂ ∂y µ ∂v ∂z + µ ∂w ∂y << ∂ ∂z 2µ ∂w ∂z -p
The first equation of the system (3.6) is therefore reduced to:

Ω ∂ ∂z 2µ ∂w ∂z φ z - ∂p ∂z φ z dΩ = Ω ρgφ z dΩ (3.7)
Thus, we have the following system, for the last two equations:

February 10, 2012            Ω ∂ ∂z 2µ ∂w ∂z - ∂p ∂z -ρg φ z dΩ = 0 Ω - ∂u ∂x - ∂v ∂y - ∂w ∂z qdΩ = 0 (3.8)
Local equation of the pressure field

From equations (3.8), we can deduce the following local equations of the higher-order model, coupling w and p :

           ∂ ∂z 2µ ∂w ∂z - ∂p ∂z -ρg = 0 - ∂u ∂x - ∂v ∂y - ∂w ∂z = 0 (3.9)
Due to the last assumption (bridging effect negligible), we saw that the third equation of the momentum balance (see equation 3.9) is reduced to:

∂ ∂z 2µ ∂w ∂z -p = ρg (3.10)
We integrate this equation between a point at elevation z and the corresponding point on the upper surface s, along the z-direction:

s z ∂ ∂z 2µ ∂w ∂z -p dz = s z ρgdz (3.11)
This is equivalent to:

2µ ∂w ∂z s -2µ ∂w ∂z z -p (s) + p (z) = ρg (s -z) (3.12)
We need here an additional assumption to determine the pressure field everywhere in the domain. We consider the air pressure to be negligible, so the boundary condition at the ice/air interface is a free surface. Thus one can easily show that the pressure field p is related to the vertical component of the velocity field w by:

p (z) = 2µ ∂w ∂z z + ρg (s -z) (3.13)
Now, the local incompressibility equation (see equation 3.9) states that in the domain Ω:

∂w ∂z = - ∂u ∂x - ∂v ∂y (3.14)
So the pressure is equal to:

p (z) = -2µ ∂u ∂x + ∂v ∂y + ρg (s -z) (3.
+ Γ b α 2 u τ φ x dΓ = Ω ρg (s -z) ∂φ x ∂x dΩ + Γw p w n x φ x dΓ + Γ i p i n x φ x dΓ Ω µ ∂u

Determination of the boundary conditions for the higher-order model

General equation for the boundary conditions

The last equation of the weak formulation (incompressibility) does not include any boundary condition and is therefore not considered in this paragraph.

We start from equation (3.16) to establish the boundary conditions. After integration by parts and considering that (u, v) is solution of the problem so they verify equation (3.17), we have:

∂Ω 2µ ∂u ∂x n x + µ ∂u ∂y + µ ∂v ∂x n y + µ ∂u ∂z n z φ x dΓ = ∂Ω ρg (s -z) φ x n x dΓ Γ b -α 2 u τ φ x dΓ + Γw p w n x φ x dΓ + Γ i p i n x φ x dΓ ∂Ω µ ∂u ∂y + µ ∂v ∂x n x + 2µ ∂v ∂y n y + µ ∂v ∂z n z -pn y φ y dΓ = Γ b -α 2 v τ φ y dΓ + Γw p w n y φ y dΓ + Γ i p i n y φ y dΓ ∂Ω µ ∂u ∂z + µ ∂w ∂x n x + µ ∂v ∂z + µ ∂w ∂y n y + 2µ ∂w ∂z n z φ z dΓ = ∂Ω ρg (s -z) φ y n y dΓ Γ b -α 2 w τ φ z dΓ Γw p w n z φ z dΓ (3.18)
We decompose ∂Ω as: ∂Ω = Γ u ∪ Γ b ∪ Γ w ∪ Γ i to establish the boundary conditions on each part of the border.

Ice/Ocean boundary condition for the higher-order model

We first consider kinematically admissible fields that are zero on all the domain border ∂Ω except at the ice/water interface Γ w , as the elevation is z = b, we have:

Γw 2µ ∂u ∂x n x + µ ∂u ∂y + µ ∂v ∂x n y + µ ∂u ∂z n z -ρgHn x -p w n x φ x dΓ = 0 Γw µ ∂u ∂y + µ ∂v ∂x n x + 2µ ∂v ∂y n y + µ ∂v ∂z n z -ρgHn y -p w n y φ y dΓ = 0 = 0 (3.19)
These integrals are zero for any kinematically admissible field whose value is zero on the boundary of the domain except for Γ w . The functions under the integral are therefore zero everywhere on Γ w and the boundary condition at the ice/water interface Γ w is:

σn = p w n (3.20)
Other boundary conditions for the higher-order model

We do the same analysis to derive the following boundary conditions on Γ b , the ice bedrock interface, Γ i , the ice front and Γ s , the ice-air interface, respectively: These static boundary conditions are completed with the kinematical ones applied on Γ u .

Summary

The higher-order model allows to decouple horizontal and vertical velocities using the following local equations: The boundary condition of this model at the ice/air interface Γ s (upper surface) is a free surface.

                                   ∂ ∂x
At the ice/water interface Γ w it is:

σn = p w n (3.25)
For the ice front Γ i , we have:

σn = p i n (3.26)
For the ice/bedrock interface Γ b :

(σn) || = -α 2 u τ (3.27)
February 10, 2012

The kinematic boundary conditions on the rest of the border Γ u are:

v = u obs (3.28)

Model validity

These equations allow the displacement to be a combination of sliding and shear deformation, the model is therefore valid in areas where the movement is mainly due to sliding but also in areas where it is only due to shear deformation (where ice is frozen at the bedrock for example).

The boundary condition at the ice/bedrock interface is:

(σn) || = -α 2 u τ on Γ b (3.29)
In the general case the ice/bedrock interface has no particular property and this condition cannot be decoupled similarly to the higher-order equations into horizontal and vertical components. Indeed the three component of u τ include part of u, v and w. To decouple this boundary condition, we must make the additional assumption that the slope of the bedrock is limite and assumed to be equal to zero (so the normal to the surface is coincident with the z axis). This assumption will be used in the next subsections.

Shelfy-stream approximation

Model assumptions

Another approximation was derived from the full-Stokes equations by [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment -Theory and application to Ice Stream-B, Antarctica[END_REF] and [START_REF] Morland | Plane and radial ice-shelf flow with prescribed temperature profile[END_REF]. This model is known as the shallow-shelf or shelfy-stream approximation and is a vertically-integrated bidimensional ice flow model.

Several assumptions are made in this model including the ones made in the higher-order model. The main one is that the basal drag associated with deforming basal sediments does not affect significantly the vertical gradients of horizontal velocities. All vertical shear is neglected so the horizontal velocity does not vary with depth. Similarly to the higherorder model, the bridging effect is negligible [van der Veen and [START_REF] Van Der Veen | Force budget: I. theory and numerical methods[END_REF]]. These assumptions are summarized as:

• ∂w ∂x << ∂u ∂z • ∂w ∂y << ∂v ∂z • ˙ xz = 0 • ˙ yz = 0
We derive the local equations and boundary conditions of this model by including the previous assumptions in the full-Stokes weak formulation.

Determination of variational and local equations associated to the shelfy-stream approximation

Variational formulation

Let us now use these assumptions in the weak formulation of the full-Stokes equations (3.3) to establish the local equations of this model.

If we include the shelfy-stream assumptions in the full-Stokes weak formulation, for any kinematically admissible field (Φ, q) T such as:

Φ ∈ V q ∈ P (3.
+ Γ b α 2 v τ φ y dΓ = Γw p w n y φ y dΓ + Γ i p i n y φ y dΓ (3.31)
These equations couple the horizontal velocity fields with the pressure.

Considering admissible regular fields Φ z = (0, 0, Φ z , 0) with Φ z = 0 on the boundary of Ω and admissible fields q, using the shelfy-stream assumptions and after integrating by parts the third system of the global full-Stokes system (3.3), we obtain the following system for the last two equations:

           Ω ∂ ∂z 2µ ∂w ∂z - ∂p ∂z -ρg φ z dΩ = 0 Ω - ∂u ∂x - ∂v ∂y - ∂w ∂z qdΩ = 0 (3.32)
February 10, 2012

Local equation of the pressure field

From equations (3.32), we can deduce the following local equations of the shelfy-stream approximation, coupling w and p :

           ∂ ∂z 2µ ∂w ∂z - ∂p ∂z -ρg = 0 - ∂u ∂x - ∂v ∂y - ∂w ∂z = 0 (3.33)
Due to the last assumption (bridging effect negligible), we saw that the third equation of the momentum balance (see equation 3.33) is reduced to:

∂ ∂z 2µ ∂w ∂z -p = ρg (3.34)
As done previously for the higher-order model (assuming the air pressure is negligible), we can deduce that the pressure is equal to:

p (z) = -

Derivation of a bidimensional model for the horizontal velocity

Depth-averaged viscosity

One consequence of the shelfy-stream approximations is that the horizontal velocity components, u and v, do not vary with depth. We can therefore reduce the space of kinematically admissible fields (φ x , φ y ) to functions that do not vary with z. From now on in this paragraph, we consider that φ x and φ y do not vary with z. However, µ varies with z as the temperature is depth dependent. So we decompose the integration over the volume Ω into one integral on the mean section ω and a second one over the thickness H. The previous equations can be modified to: As we have H = s -b, the equation gives:

Ω ρg ∂s ∂x φ x dΩ = ω ρgH ∂s ∂x φ x dω (3.42)
The same integration is performed for the y component.

Depth-integration of ice front terms

We need an additional assumption to depth-integrate the ice front terms. We suppose that the air or water pressure at the ice front Γ i for a point at elevation z is:

p i = ρ w g min (z, 0) (3.43)
For any point below the water level, p i is the water pressure. For points whose elevation is above the water pressure, p i = 0.

As we consider the ice front to be vertical cliffs, it gives:

Γ i p i n x φ x dΓ = γ i 0 b ρ w gzn x φ x dzdl (3.44) if we note γ i = Γ i ∩ ∂ω.
The outward pointing normal does not vary with depth as the ice front is a vertical cliff, so:

Γ i p i n x φ x dΓ = γ i - 1 2 ρ w gb 2 n x φ x dl (3.45)
The same integration is performed for the y component.

Using these two equalities, the first two equations of the full-

Other boundary terms

We now decompose the integral over ∂Ω as

∂Ω = Γ u ∪ Γ w ∪ Γ b ∪ Γ s ∪ Γ i : ∂Ω ρg (s -z) φ x n x dΓ = Γ b ∪Γw ρg (s -z) φ x n x dΓ (3.47)
On Γ s , z = s so the part of the integral over Γ s is zero. φ x and φ y are zero on Γ u as Dirichlet conditions are used on Γ u so the fields φ x and φ y have to be zero on Γ u in order to be part of the kinematically admissible fields, this part of the integral is also zero. Let us now consider the part of the integral over the ice front Γ i :

Γ i ρg (s -z) φ x n x dΓ = γ i ρg s b (s -z) dz φ x n x dl (3.48)
as the normal of the ice front does not vary with depth. It gives:

Γ i ρg (s -z) φ x n x dΓ = γ i ρg H 2 2 φ x n x dl (3.49) Γ i ρg (s -z) φ y n y dΓ = γ i ρg H 2 2 φ y n y dl (3.50)
On Γ b and Γ w , as z = b, this integral is:

Γ b ∪Γw ρg (s -z) φ x n x dΓ = Γ b ∪Γw ρgHφ x n x dΓ (3.51) Γ b ∪Γw ρg (s -z) φ y n y dΓ = Γ b ∪Γw ρgHφ y n y dΓ (3.52)
So the momentum balance equations we obtain after using these relations are:

Parametric surface

To obtain the local equations of the shelfy-stream approximation, all the integrals must be expressed on the same domain. We use the expression of a parametric surface integral. Let f (x, y) be a parameterization of a surface S, g a function defined on S and P the projection of the surface S on the plane z = 0. The integral of g over S is:

S g (x, y, z) dS = P g (x, y, f (x, y)) 1 + ∂f ∂x 2 + ∂f ∂y 2 1/2 dA (3.54)
As the terms in the integrals over Γ b do not vary with z, if we use this relation to change the integral over Γ b into an integral over ω, it gives:

Γ b ∪Γw ρgHφ x n x dΓ = ω ρgHφ x n b x 1 + ∂b ∂x 2 + ∂b ∂y 2 1/2 1 + ∂f ∂x 2 + ∂f ∂y 2 1/2 dΓ (3.55)
where n b x and n b y are components of the outward pointing normal of the lower surface b and f = (s -b) /2 to define the z coordinate of the mean section. To simplify the notation, we introduce J as:

J = 1 + ∂b ∂x 2 + ∂b ∂y 2 1/2 1 + ∂f ∂x 2 + ∂f ∂y 2 1/2
(3.56)

If we note ω b the part of the domain on the ice sheet and ω w the part of the domain on the ice shelf, the other integrals over Γ b and Γ w are similarly rewritten as:

Γ b ∪Γw ρgHφ y n y dΓ = ω ρgHφ y n b y JdΓ (3.57) Γ b α 2 u τ φ x dΓ = ω b α 2 u τ φ x JdΓ (3.58) Γ b α 2 v τ φ y dΓ = ω b α 2 v τ φ y JdΓ (3.

Determination of the local equations for the bidimensional shelfystream approximation

The weak formulation for the shelfy-stream approximation is therefore:

ω ∂ ∂x
       σn • e x = 1 2 ρgH 2 -ρ w gb 2 n x σn • e y = 1 2 ρgH 2 -ρ w gb 2 n y
(3.66)

Summary

The shelfy-stream approximation is a bidimensional model that allows to compute the horizontal and vertical velocity successively. The local equations of this model are described by equation (3.63) and the boundary conditions ar the ice front by equation (3.66). Dirichlet boundary conditions are used on the rest of the border.

Simplification

As for the higher-order model, these equations are not completely decoupled yet as the velocity tangential to the ice/bedrock interface u τ and v τ include in the general the vertical velocity w. As mentioned above for the higher-order model, we assume that the bedrock surface gradient is limited to decouple the friction equation.

The local equations of the shelfy-stream approximation are often written as [START_REF] Macayeal | Large-scale ice flow over a viscous basal sediment -Theory and application to Ice Stream-B, Antarctica[END_REF][START_REF] Macayeal | A tutorial on the use of control methods in ice-sheet modeling[END_REF][START_REF] Hindmarsh | A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling[END_REF][START_REF] Gudmundsson | Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation[END_REF][START_REF] Sergienko | Ice stream basal conditions from block-wise surface data inversion and simple regression models of ice stream flow: Application to bindschadler ice stream[END_REF][START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF]: This form is valid if the slope of the ice lower elevation is limited, so that n b x and n b y can be neglected. It also supposes that the area of the ice lower elevation and mean section are almost identical, so that the coefficient J we introduced is close to 1. This later assumption is true if in addition to small slopes in the ice lower surface, upper surface also experience limited slopes for the grounded part.

Model validity

This model can be used for ice shelves and ice streams sliding on the bed as shear deformation is negligible in these cases. In areas where there is no sliding (frozen bed), the flow is mainly due to shear deformation and the approximations made here are not valid anymore. This model cannot be used in these areas and we need a higher-order model, which includes vertical shear.

Shallow ice approximation

Model assumptions

The shallow-ice approximation is a three dimensional model widely used in ice sheet modeling that was presented by [START_REF] Hutter | Theoretical glaciology: material science of ice and the mechanics of glaciers and ice sheets[END_REF]. While derived from the full-Stokes equations, similarly to the shelfy-stream approximation and higher-order model, the assumptions made in this model are completely different from the previous two models. The only terms that are not neglected in the deviatoric stress tensor are σ xz and σ yz . The horizontal gradients of vertical velocity are also neglected compared to the vertical gradients of horizontal velocities. We can summarize these assumptions as:

• ∂σ xx ∂x ,

Variational formulation

As for the derivation of the previous simplified models, let us use these assumptions in the weak formulation of the full-Stokes equations (3.3) to establish the local equations of this model.

If we include the last two assumptions of the higher-order model (∂w/∂x << ∂u/∂z and ∂w/∂y << ∂v/∂z), for any kinematically admissible field (Φ, q) T such as:

Φ ∈ V q ∈ P (3.
+ Γ b α 2 v τ φ y dΓ = Γw p w n y φ y dΓ + Γ i p i n y φ y dΓ (3.69)
These equations couple the horizontal velocity fields and the ice pressure.

In terms of velocity components, the first three assumptions of the shallow ice approximation are: To include these conditions in the full-Stokes weak formulation, we first need to integrate by part the first terms of the first three equations.

• ∂ ∂x
Considering admissible regular fields Φ x = (Φ x , 0, 0, 0), Φ y = (0, Φ y , 0, 0) and Φ z = (0, 0, Φ z , 0), with Φ x = Φ y = Φ z = 0 on the boundary of Ω and admissible fields q, and after integrating by parts the first three systems of the global full-Stokes system (3.3), we obtain:

                                   Ω ∂p ∂x - ∂ ∂z µ ∂u ∂z φ x dΩ = 0 Ω ∂p ∂y - ∂ ∂z µ ∂v ∂z φ y dΩ = 0 Ω ∂p ∂z + ρg φ z dΩ = 0 Ω - ∂u ∂x - ∂v ∂y - ∂w ∂z qdΩ = 0 (3.70)

Local equation of pressure field

The third equation states that ∂p ∂z = -ρg (3.71)

If we integrate it between the point at altitude z and the upper surface, and use the additional assumption that the air pressure at the ice/air interface is neglected, the pressure at elevation z is such that:

p (z) = ρg (s -z) (3.72)
where s is the upper surface elevation.

Local equations of the horizontal velocity

Using equation (3.72) into the variational formulation of the horizontal velocity, we obtain:

          

Boundary conditions for the shallow ice approximation

The shallow ice equations (3.73) are integrated along vertical columns between the bedrock and a given elevation z. A Dirichlet condition is needed for each column, usually applied at the base of the ice, generally taken as zero if the ice is below the pressure melting point, or derived from a simplified friction law to account for basal sliding.

Summary

The shallow ice is a simple approximation of the full-Stokes equations that allows efficient computing of ice velocity. Its local equations are:

                                 ∂ ∂z µ ∂u ∂z -ρg ∂s ∂x = 0 ∂ ∂z µ ∂v ∂z -ρg ∂s ∂y = 0 ∂p ∂z + ρg = 0 ∂u ∂x + ∂v ∂y + ∂w ∂z = 0
(3.74)

Model validity

The shallow ice is a model valid for ice sheets where movement is only due to vertical shear deformation. This model is not appropriate to model ice shelves or fast ice streams where most of the motion is caused by sliding over the bedrock or lateral spreading.

Domains of validity

We derived four ice flow models from the full-Stokes equations. Simpler models reduce the computational cost while sophisticated ones allow a better representation of the actual physics. The The first approximation widely used in ice sheet modeling was the shallow ice approximation. Its simplicity and low computational cost allowed it to be used in large-scale models [START_REF] Ritz | Sensitivity of a greenland ice sheet model to ice flow and ablation parameters: Consequences for the evolution through the last climatic cycle[END_REF][START_REF] Deconto | Rapid Cenozoic glaciation of Antarctica induced by declining atmospheric CO2[END_REF]. This approximation takes advantage of the thickness-length ratio and the main consequence is that ice deformation is only due to vertical shear stress. For this reason, this approximation is not valid for wavelengths smaller than a few ice thickness. Displacement from basal sliding is also often neglected in this approximation, so results are more accurate for areas where ice is frozen at the bedrock/ice interface.

The second approximation, shallow-shelf or shelfy-stream approximation, was developed to model ice shelves and ice streams. This approximation neglects vertical shear but allows to include basal sliding. As ice shelves movement is mainly due to spreading and ice streams displacement primarily slides over soft beds, so this approximation should advantageously be used for this type of regions. The computational cost is also very low as it is a bidimensional model. This approximation has been extensively used to infer basal parameters of Antarctic and Greenland ice stream [START_REF] Macayeal | The basal stress distribution of Ice Stream E, Antarctica, Inferred by Control Methods[END_REF][START_REF] Joughin | Observation and analysis of ice flow in the large Greenland ice stream[END_REF][START_REF] Vieli | Application of control methods for modelling the flow of Pine Island Glacier, West Antarctica[END_REF][START_REF] Joughin | Integrating satellite observations with modelling: basal shear stress of the Filcher-Ronne ice streams, Antarctica[END_REF][START_REF] Joughin | Basal conditions for Pine Island and Twaites Glaciers, West Antarctica, determined using satellite and airborne data[END_REF] as it allows to include basal sliding and properties of ice shelves [START_REF] Rommelaere | Large-scale rheology of the Ross Ice Shelf, Antarctica, computed by a control method[END_REF][START_REF] Larour | Modélisation numérique du comportement des banquises flottantes, validée par imagerie satellitaire[END_REF][START_REF] Khazendar | Larsen B Ice Shelf rheology preceding its disintegration inferred by a control method[END_REF][START_REF] Khazendar | Roles of marine ice, rheology, and fracture in the flow and stability of the Brunt/Stancomb-Wills Ice Shelf[END_REF].

The higher-order model was introduced later and includes some terms of the vertical shear. Therefore, this model should not only allow to reproduce ice shelves and ice streams but should be able to capture the evolution of ice streams within the ice sheet. The computational cost associated to this model is much more intensive than that of the previous ones and this model has rarely been used.

To better understand the limitations of these approximations, some theoretical studies have been performed by [START_REF] Hindmarsh | A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling[END_REF] and [START_REF] Gudmundsson | Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation[END_REF]. They compare results of several approximations to the full-Stokes solution by analyzing the influence of small perturbations to the ice flow. [START_REF] Hindmarsh | A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling[END_REF] uses an infinite plane with a small inclination and compare several approximations among other are the shallow ice approximation, the shelfy-stream approximation and the higher-order model. Model computations of these approximations are compared to analytical solutions and to the full-Stokes solution.

Hindmarsh [2004] shows that for a linear rheology, the solution with the shallow ice model deviates from the full-Stokes when the wavelength of the perturbation is smaller than 50 ice thickness and the solution of higher-order model deviates for wavelength smaller that 5 to 10 ice thickness. The value at which the higher-order solution starts to deviate from the full-Stokes one depends on the slip ratio. For non-linear rheology, [START_REF] Hindmarsh | A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling[END_REF] shows that the four models we presented are stable. However, the shelfy-stream approximation gives inaccurate results except for high slip ratio, which is not unexpected for an ice shelf and ice stream model. Both shallow ice and shelfy-stream approximations fail for wavelength smaller than 5 to 10 ice thickness, but the shallow ice approximation fails in a more physical way, as it tends to attenuate the errors whereas the shallow-shelf makes them persist. The higher-order model is able to correctly reproduce the full-Stokes solution for perturbations with wavelength below one ice thickness, depending on the slip ratio or the coefficient used in the friction law. Inclusion of longitudinal stresses increase accuracy for short wavelengths but still does not allow to reproduce ice behavior for wavelengths shorter than one ice thickness.

Results from [START_REF] Gudmundsson | Analytical solutions for the surface response to small amplitude perturbations in boundary data in the shallow-ice-stream approximation[END_REF] compare the shallow ice and shelfy-stream approximation to a full-Stokes model for a linear rheology and a linear friction law. According to his experiments, transmission of basal topography perturbations to the surface topography are correctly reproduced for perturbations with wavelengths higher than 300 ice thickness for the shallow ice model and 5 to 8 ice thickness for the shelfy-stream approximation. Transmission of perturbations in the basal sliding are better transferred by the shallow ice model, as results from shelfy-stream model slightly differ from the full-Stokes solution for all wavelength. Finally, transmission of basal sliding perturbations to the surface velocity are correctly reproduced with the shelfy-stream approximation for the horizontal velocity (with a small offset) but not for the vertical velocity, as the solution deviates for wavelength under 500 ice thickness.

All these comparisons have been made for theoretical cases. A study of Pine Island Glacier, West Antarctica, by [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] indirectly compares results from the shelfystream approximation, higher-order and full-Stokes models by reproducing the surface velocity using control methods. They find that the inferred patterns of basal drag found with the three models are similar except in the vicinity of the ground line. These results are used and discussed in chapter 6.

Chapter summary

In this chapter we described three common approximations used to model ice sheet flow, all deriving from the full-Stokes equations. Under classical assumptions, but also under a geometrical one (that seems not to be always sufficiently discussed in the litterature), we established the local equations and boundary conditions for these three model using the variational formulation of full-Stokes as it is a more general and more rigorous formulation for a given continuum mechanical problem. We then detailed the domain of validity for these three models based on the literature. In the following chapter we detail how these formulations are implemented using the finite element method.

In the previous chapters, we detailed the thermal, mechanical and mass conservation equations that govern the physics of ice sheet systems independently. The general thermomechanical evolutive equations that couple these three aspects and that must be solved for are:

∀x ∈ Ω, find (u, v, w, p, T, H) :

                                                         ∂ ∂x
∂T ∂t = -v • ∇T + k ρc ∆T + Φ ρc ∂H ∂t = -∇ • (vH) + Ṁs -Ṁb (4.1)
In this chapter, we focus on the numerical aspects of the problem. After describing the time and space discretizations we adopted, we derive the weak formulations for these two problems and detail the algorithm used to solve the non-linearity introduced by the behavior law in the mechanical problem and the non-linearity introduced by the phase change of ice at pressure melting point in the thermal problem. We then explain how the thermal, mechanical and mass balance problems are combined to perform either steady-state or evolutive simulations. We finally present the model implementation and the main characteristics of the software and show some simple experiments commonly used to verify the validity of the solutions.

Discretization

Space discretization

Glaciers have a very small aspect ratio so the vertical and horizontal discretizations must be different. Here we use triangles to mesh the bidimensional space with the Triangle library [START_REF] Shewchuk | Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator[END_REF]. We use a mesher inspired from BAMG and YAMS [START_REF] Frey | a fully automatic adaptive isotropic surface remeshing procedure[END_REF][START_REF] Hecht | BAMG: Bi-dimensional anisotropic mesh generator[END_REF]) to generate anisotropic triangular meshes whose metric is based on the velocity hessian matrix. This technique distorts the mesh and reduces the discretization errors, while minimizing the number of elements. Interpolation-based a-priori error estimates of a finite element solution based on piecewise linear P 1 elements depends only on its Hessian, provided that the solution is regular enough [START_REF] Habashi | Anisotropic mesh adaptation: towards user-independent, mesh-independent and solverindependent CFD. Part I: general principles[END_REF]. For each element E of the mesh, the error between the P 1 interpolated field u h and the exact field u is bounded as follows:

|u (x) -u h (x)| ≤ Ch 2 sup x∈E |H u (x)| (4.2)
where C is a constant depending on the space dimension, h the characteristic length of the element, H u (x) the Hessian matrix of u (x) and |H u (x)| its spectral norm. Adaptive mesh refinement takes advantage of this property to generate meshes whose element sizes are such that the interpolation error is constant for the entire mesh. Let us observe that other error indicators are available in the literature.

Vertical planes passing through the sides of the triangular edges divide the domain into vertical columns of ice. These vertical columns are limited by the ice bottom on their lower surface and the ice surface on the upper part. We then divide the ice thickness into a set of layers of horizontally variable thickness. Each vertical column is therefore divided into a given number of triangular prisms whose upper and lower surfaces are not necessarily parallel. The three-dimensional mesh is therefore composed of prismatic elements.

The same mesh is used to compute the thermal, mechanical and prognostic models even though it has been adapted only for the mechanical model. Since the variations in velocities are larger and often related to variations in other fields, a mesh adapted to velocity changes is generally also suitable for thermal and prognostic models.

Time discretization

We use finite differences in time for the time discretization. The time discretization is necessary for transient models that represent the temporal evolution of a glacier or ice sheet (evolution of thickness and temperature).

∂f ∂t f i+1 -f i ∆t (4.3)
where ∆t is the time step between t i and t i+1 . We employ the Courant-Friedrichs-Lewy (CFL) condition for stability [START_REF] Courant | Über die partiellen differenzengleichungen der mathematischen physik[END_REF][START_REF] Courant | Variational methods for the solution of problems of equilibrium and vibrations[END_REF].

Mechanical models implementation

We detail here the implementation of the full-Stokes model, the weak formulation and the Galerkin approximation. The derivation of the other ice flow models implementation is similar and is therefore not detailed here. The fundamental problem description, which includes the local equations and boundary conditions, is equivalent to its weak formulation.

Weak formulation

We use the local equations (2.30) and boundary conditions (section 2.3.3) defined in chapter 2 to derive the weak formulation of the full-Stokes equations.

Let V u obs be the space of kinematically admissible velocity fields and P be the space of admissible pressure fields.

V u obs = Φ ∈ H 1 (Ω) 3 , Φ = u obs on Γ u and Φ • n = 0 on Γ b P = q ∈ L 2 (Ω) (4.4)
As done in section 3.1, we remind that in order to simplify the presentation, without loss of generality, we will consider the case u obs = 0 and equate the space of velocity solution and kinematically admissible velocity fields by setting V = V u obs . In practice, non homogeneous Dirichlet conditions based on observed velocities are used and included using lifting.

For any kinematically admissible field (Φ, q) ∈ (V × P), the weak formulation of the full- After integration by parts and using the boundary conditions in section 2.3.3, we obtain the weak formulation of the full-Stokes model as presented in equation (3.2).

For a given, strictly positive viscosity coefficient µ and regular data, this mixed problem can be proven to be well-posed by applying basically theories of Babuska and Brezzi (see e.g. [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]).

Galerkin approximation

The finite element method is based on the discretization of the solution space. We use V h and P h two sub-spaces of respectively V and P of finite dimension. We suppose that {Φ 1 , ..., Φ n } and {q 1 , ..., q m } are two bases of these spaces. We want the weak formulation established in the previous paragraph to be valid for any test functions, this being true in particular for the basis functions. ∀ {i, j} ∈ {

j dΩ + Γ b α 2 u τ • Φ i dΓ = Γw p w n • Φ i dΓ + Γ i p i n • Φ i dΓ - Ω ρgφ iz dΩ (4.6)
We decompose the solution of the problem on these two bases:

v (x) = n k=1 v k Φ k (x) p (x) = m l=1 p l q l (x) (4.7)
Using this decomposition, the weak formulation is:

∀ {i, j} ∈ {[1 : n] × [1 : m]} n k=1 v k Ω 2µ
∂z q j dΩ + Γ b α 2 Φ kτ • Φ i dΓ - m l=1 p l Ω q l ∂φ ix ∂x + ∂φ iy ∂y + ∂φ iz ∂z dΩ = Γw p w n • Φ i dΓ + Γ i p i n • Φ i dΓ - Ω ρgφ iz dΩ (4.8)
where Φ kτ is part of the basis function Φ k tangential to the bedrock surface.

This formulation is valid for any basis function Φ i ∈ V and q j ∈ P. It therefore gives a system of n + m equations of n + m unknowns. In matrix form, this problem is equivalent to finding the vectors (U, P ) solution of the following linear system of equations:

A B B T 0 U P = F 0 (4.9)
with:

A k,i = Γ b α 2 Φ kτ • Φ i dΓ + Ω 2µ
F i = Γw p w n • Φ i dΓ + Γ i p i n • Φ i dΓ - Ω ρgφ iz dΩ (4.12)

Finite element method

The weak formulation of the full-Stokes problem is a mixed formulation as both velocity and pressure are solved. As observed before, the continuous problem has a unique solution, however this is not true for any discretization of this continuous problem and a particular attention must be paid to the choice of discretization spaces [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]. The discretization space adopted defines prismatic element for the three dimensional models.

The basis functions used in the finite element method must ensure that the problem will be numerically stable. It is for example not possible to choose polynomials that have the same order for both velocity and pressure. The choice of finite element must fulfill, in particular, the Ladyzhenskaya-Babuška-Brezzi compatibility condition (LBB) that restricts the approximation spaces in order to ensure the existence of a unique solution to the discrete problem.

Let V and P be two Hilbert spaces, || || V and || || P the associated norms and b a linear form continuous on V × P. The form b satisfies the LBB condition if there is a scalar β > 0 such that:

inf p∈P sup v∈V b (v, p) ||v|| V ||p|| P ≥ β (4.13)
To satisfy the LBB condition we employ the Residual Free Bubble Function (RFBF) stabilization technique [Brezzi et al., 1996b, a]. This is one of the methods that ensures avoidance of spurious pressure modes and numerical instabilities associated with inadequate spatial resolution. This approach has the advantage of providing control over the accuracy of FE schemes and has been used in Earth Science modeling, ocean modeling in particular [START_REF] Nechaev | A diagnostic stabilized finite-element ocean circulation model[END_REF][START_REF] Guillen-Gonzalez | Bubble finite elements for the primitive equations of the ocean[END_REF]. The major idea of the RFBF stabilization is to enrich the approximation function space of velocities with "bubble functions" which have zero values on the boundaries of the elements but contribute to the system matrix through their projections on the original basis functions.

In our case, we use P + 1 P 1 (MINI) elements. An additional node is added in the center of the prism with a quadratic nodal function for the velocity, whose value is zero on the regular prism nodes (known as bubble function in the literature). The nodal functions used for the other nodes (the six vertices of the prism) are linear for the pressure and the velocity. To avoid including these additional nodes and the degrees of freedom associated with these nodes, we use static condensation (see e.g. [START_REF] Franca | Bubble functions prompt unusual stabilized finite-element methods[END_REF]).

The weak formulations and Galerkin approximations of the Shelfy-stream approximation and Higher-order models are derived the same way. These two problems only involve degrees of freedom for velocity and are therefore not mixed problems. Piecewise linear triangular P 1 elements are used for the discretization of the shelfy-stream approximation and P 1 prismatic elements for the higher-order model.

Non linear behavior law

In chapter 2, we introduced Glen's flow law as the deformation law for ice flow. This is a non-linear flow law, and therefore leads to non-linear partial differential equations for the velocity. The non-linearity is solved iteratively using the Picard's iterative method, a fixed point-like strategy. [START_REF] Reist | Mathematical analysis and numerical simulation of the motion of a glacier[END_REF] (p.39) showed that this method can be used for ice flow equations and converges if we start from an initial value "close enough" to the solution.

To prove the convergence of the continuous problem, [START_REF] Reist | Mathematical analysis and numerical simulation of the motion of a glacier[END_REF] used the calculus of variations. He defines the application that for a given u h and its associated viscosity µ (u h ) returns the velocity v h solution of the linear problem described above (see equation (3.2)). Note that a solution of the problem is a fixed-point of this application. [START_REF] Reist | Mathematical analysis and numerical simulation of the motion of a glacier[END_REF] then considers the Frechet-Derivative of this application and observed that its spectral radius is inferior to 1. As this funtion is also the iteration function in the algorithm applied (see below), if we start from an initial guess not too far from the solution, the suite of functions solution of this application converges toward u h .

We solve a linear problem for each iteration until convergence. We start with an estimate of the velocity v 0 and do the following steps until the convergence criterion is met:

• Calculate new viscosity µ k = µ (V k ) • Compute velocity V k+1 such that K (V k ) V k+1 = F When the residual term ||K (V k ) V k -F ||/||F ||
is below the convergence criterion, we stop the iterations.

As the shelfy-stream and higher-order models decouple the horizontal and vertical equations, an additional step is needed to compute the vertical velocity as a by-product of the horizontal velocity contrary to full-Stokes model. This is done using the incompressibility equation with the finite element method, assuming that there is no penetration of ice in the bedrock so the boundary condition at the ice/bedrock interface is a non-penetration condition. Figure 4.2 shows the algorithm used for both the shelfy-stream or higher-order models and for full-Stokes equations, with u h representing the horizontal components of velocity.

Thermal model implementation

The heat budget equation established in the first chapter ( 2) is:

∂T ∂t + v • ∇T = k ρc ∆T + Φ ρc (4.14)
However the ice temperature cannot be above pressure melting point as we saw in section 2.4.3, and appropriate model must be used to correctly model ice temperature. We use here a cold ice model and assume that the areas of temperate ice are limited in space and located on the ice/bedrock interface. This assumption is realistic for most parts of the Antarctic ice sheet but not always true for Greenland (see section 2.4.3). To avoid temperatures that are above the pressure melting point, we use the penalty method.

Including phase change at the ice/bedrock interface

We include here the solid-liquid phase change processes at the ice/bedrock interface. If the ice temperature reaches the pressure melting point, it starts melting and the heat budget equation has to be modified to include the phase change. The temperature/melting condition is described in scheme 4. The ice temperature needs therefore to respect the following condition: T ≤ T pmp . If the temperature reaches the pressure melting point, the heat equation has to include a melting term that represents the heat loss. We consider here that ice is below the metling point everywhere except at the ice/bedrock interface. The boundary condition on the lower surface Γ b must therefore be changed to account for ice melting in areas where ice reaches the pressure melting point:

k ∇ T | b • n = G -τ b • u τ + ρL Ṁb (4.15)
where L is the ice latent heat of fusion and Ṁb the ice melting rate in m/s, positive or zero.

The general thermal boundary condition at the ice/bedrock interface Γ b is:

February 10, 2012 k ∇ T | b • n = G -τ b • u τ +    ρL Ṁb if T = T pmp 0 if T < T pmp (4.16)

Signorini formalism

To include this melting condition, we use a formalism similar to Signorini's contact law or Kuhn and Tucker conditions (see e.g. [START_REF] Duvaut | Inequalities in Thermoelasticity and Magnetohydrodynamics[END_REF]).

In our problem, the conditions to fulfill are:

T ≤ T pmp ⇔ Temperature below or at the pressure melting point Ṁb ≥ 0 ⇔ Positive or zero melting (T -T pmp ) Ṁb = 0 ⇔ Condition of complementarity (4.17) Inequalities are not easy to introduce in weak formulations. We therefore use a penalty method (e.g., [START_REF] Courant | Variational methods for the solution of problems of equilibrium and vibrations[END_REF]) and introduce a Sign field S (see Ben Dhia [1988]; Ben Dhia and Zammali [2007]) to approximate them by equations:

ρL Ṁb = κ S × (T -T pmp ) (4.18)
where the sign field function S is equal to 1 R + (T -T pmp ), 1 R + being the indicator function of the sem-axis of positive reals, defined by:

1 R + (x) = 1 if x ≥ 0 0 else (4.19)
If κ tends towards infinity, we recover Signorini's formalism (see figure 4.4). In our problem, the unit of κ is W m -2 K -1 , so κS (T -T pmp ) is equivalent to a heat flux. κ is equivalent to a heat transfer coefficient in a Newton's cooling law. Indeed in a Newton's cooling law, the energy, E, is solution of the differential equation such as:

dE dt = -hA∆T (t) (4.20)
where h is the heat transfer coefficient, A the surface area of the heat being transfered and ∆T the time-dependant thermal gradient between an object and its environment.

As for contact problems, one can also use a Lagrangian formalism that is more rigorous (e.g. [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]; [START_REF] Wriggers | Computational Contact Mechanics[END_REF]; Ben Dhia and Zarroug [2002]; Ben Dhia and Zammali [2007]).

Weak formulation

The heat budget equation is, find T and S such that:

∂T ∂t + v • ∇T = k ρc ∆T + Φ ρc (4.21)
The boundary conditions of the model are:

                   Upper surface Γ s : T = T air Icesheet base Γ b : k∇ T | b • n = G -τ b • u τ + κ S × (T -T pmp ) Iceshelf base Γ w : k ∇ T | b • n = -ρ w c pM γ (T b -T pmp ) Other borders Γ u ∪ Γ i : ∇ T • n = 0 (4.22)
where Γ s is the upper surface, Γ b the ice/bedrock interface, Γ w the ice/water interface, T air the air temperature at the surface, T b the ice temperature on the lower surface, T pmp the pressure melting point, G the geothermal flux, τ b the friction stress at the ice/bedrock interface, k the ice thermal conductivity, u τ the ice velocity tangential to the ice/bedrock interface, n the outward pointing normal vector, ρ w the water density, c pM the mixed layer specific heat capacity [START_REF] Holland | Modeling thermodynamic ice-ocean interactions at the base of an ice shelf[END_REF] and γ the thermal exchange velocity.

Similarly to the previous section, we establish the weak formulation of this model. Let V T air be the space of kinematically admissible temperature fields.

V T air = θ ∈ H 1 (Ω) , θ = T air on Γ s (4.23)
Similarly to the mechanical problem, we consider here that homogeneous Dirichlet conditions are applied. For any admissible temperature fields θ ∈ V 0 , the weak formulation of the thermal equation is:

Ω ∂T ∂t θ + k ρc ∇T • ∇θ + θv • ∇T dΩ + Γw ρ w c pM γ ρc T θdS = Ω Φ ρc θdΩ + Γw ρ w c pM γ ρc T pmp θdS + Γ b 1 ρc (G -τ b • u τ + κ S × (T -T pmp )) θdS (4.24)
We use finite differences in time and an implicit time scheme. The Galerkin approximation is derived similarly for the thermal equation as for the mechanical equations. The weak formulation becomes:

Find T ∈ V, ∀ θ ∈ V a (T, θ) = l (θ) (4.25)
with:

a (T, θ) = Ω 1 ∆t T θ + θv • ∇T + k ρc ∇T • ∇θ ∂Ω + Γw ρ w c pM γ ρc T θdS - Γ b 1 ρc (κ ST ) θdΓ (4.26) l (θ) = Ω 1 ∆t T -θdΩ + Φ ρc θ dΩ + Γw ρ w c pM γ ρc T pmp θdΓ + Γ b 1 ρc (G -τ b • u τ -κ ST pmp ) θdΓ (4.27) S = 1 R + (T -T pmp ) (4.28)
where T -is the ice temperature at the previous time step. The non-linearity induced by the constraint that ice temperature must remain below the pressure melting point is discussed in the next paragraph.

Thermal non-linearity

Ice melting introduces a non-linearity in the problem as areas where ice is at the pressure melting point are not a priori known. This thermal non-linearity is raised up by the unknown sign field S and we use the following algorithm to solve this nonlinear equation for one particular time step (see also figure 4.5):

1. Initial field S = S 1 2. Compute solution T 1 and update

S = S 2 = 1 R + (T 1 -T pmp ) 3. If S 2 = S 1 stop, else iterate until S i = S i-1
The melting is computed as a by-product:

Ṁb = κ ρL S × (T -T pmp ) (4.29)

Stabilized formulation

If the convective coefficient is larger than the diffusive coefficient, the Galerkin method leads to oscillations. The stabilized formulation includes an additional term so the problem we solve for is:

Find T ∈ V, ∀θ ∈ V A (T, θ) = L (θ) (4.30) with: A (T, θ) = a (T, θ) + b (T, θ) L (θ) = l (θ) + m (θ) (4.31) and b (T, θ) = Ω τ - k ρc T + v • ∇T v • ∇θ m (θ) = Ω τ Φ ρc v • ∇θ (4.32)
with τ the stabilization parameter that is detailed in the next section.

Stabilization parameter

The parameter τ that appears in the stabilized formulation depends on the mesh properties, the velocity and thermal parameters. It is expressed as:

τ = h K 2||v|| min 1, h K ||v||m K 2κ (4.33)
with the parameters:

• h K is the element diameter, the minimal length of the element edges

• m K is the stabilization parameter • ||v|| is the velocity norm u 2 + v 2 + w 2 1/2 • κ = k ρc
is the diffusivity coefficient of the thermal equation

The stabilization parameter depends on the type of elements and is equal to 1/3 for wedges elements.

To avoid dividing by a zero if the velocity is zero, we use:

• τ = h K 2||v|| if h K ||v||m K 2κ > 1 • τ = h 2 K m K 4κ if h K ||v||m K 2κ < 1
The velocity norm is computed on each gauss point and the parameter must therefore be computed at each gauss point.

Combining thermal and mechanical models 4.4.1 Steady-state solution

The thermal and mechanical equations are coupled by the viscosity parameter (temperature dependent), the advection and deformational heating. In order to have a thermo-mechanical steady state, we need to reach convergence of both velocity and temperature.

Here, we use a fixed point scheme (Picard scheme), and we solve alternatively for the velocity and the temperature until both fields have converged according to a convergence criterion.

We start with a priori estimates of velocities and temperature. Since solving temperature and velocity is also a non-linear problem and requires iterative schemes, a double loop is needed to solve the thermo-mechanical steady-state equations.

To reach the steady-state equilibrium, we use the convergence scheme described in figure 4.6. For each step, we compute the velocity and the temperature as explained above. We stop the iterations when the relative variations of both temperature and velocity between two iterations are below a specified tolerance.

Evolutive models

This section details the steps performed for the resolution of transient models. All the equations are resolved independently to reduce the numerical cost. Changes in ice thickness, H, will change the forces applied to the ice sheet system. For example a change in surface slope modifies the driving stress, which will change the ice velocity.

We start with a given configuration of the glacier where the ice geometry (thickness, surface, bedrock) and temperatures are know. The velocity can be estimated using observations. For each time step, we start by computing the new temperature. The ice viscosity parameter, B, is updated accordingly. The new velocity is then calculated. We use the mass conservation to calculate the thickness for the next time step and change the geometry of the glacier accordingly. The mesh is finally updated to match the new configuration (horizontal layers follow the new surface and bed elevations).

The time steps are either computed based on the CFL condition or prescribed (in this case we ensure that the prescribed time step is compatible with CFL condition). Scheme 4.7 summarizes these steps.

We use the Arbitrary Eulerian-Lagrangian method [START_REF] Donea | Advances in computational mechanics[END_REF] to update the mesh at each time step. It is convenient for our problem to use a mesh that is horizontally fixed (Eulerian approach) but that evolves vertically as the lower and upper surfaces change with time. The only equation that is affected by the mesh velocity is the thermal model, for which the advection needs to be modified to take into account the vertical velocity of the mesh.

Initial state: v 0 , H 0 and T 0

Energy balance

Final state: v f , H f and T f

T i+1 = f e (T i , v i , H i ) t < t f Momentum balance v i+1 = f m (v i , H i , T i+1 )
Mass balance

H i+1 = f ma (H i , v i , T i+1 )
Update mesh

Model implementation

All these algorithms have been implemented in the Ice Sheet System Model (ISSM) developed in collaboration between the Jet Propulsion Laboratory, University of California Irvine and the MSSMat Laboratory at École Centrale Paris [START_REF] Larour | Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM)[END_REF]. ISSM is a thermo-dynamic finite element model massively parallelized. We describe here its main properties.

Languages

ISSM uses the C language for the implementation of the finite elements. Objects management relies on the C++ language so that polymorphic capacities can be used. The advantage of this solution is to be both flexible and scalable.

For ease of use, a Matlab interface has been developed. It allows easy pre and postprocessing. The C/C++ core is interfaced to the Matlab environment using Matlab's external API (Application Program Interface). This results in a seamless integration of ISSM in Matlab.

Parallel architecture

ISSM can be used in a serial mode within Matlab but is above all a parallel architecture designed to run models on large clusters. When running on these large clusters, the C/C++ core is compiled as a standalone executable. Parallelism is achieved through the Message Passing Interface (MPI) so the software can be used on both distributed or shared memory clusters.

The PETSc library is used to define numerical objects such as vectors, matrices and solvers. These objects are directly used in ISSM and are abstracted to hide the difference between the serial and parallel implementation. PETSc also provides access to a large array of solvers.

We mainly rely on MUMPS (Multifrontal Massively Parallel Sparse direct solver) to compute the models solutions. This solver is poorly scalable but does not suffer from convergence issues as it relies on a direct solving method. The solving phase represents more than 90% of the computational time when using this direct solver but still allows us to tackle problems with several millions of degrees of freedom. Further work will be necessary to use iterative solvers, especially to solve the full-Stokes equations that require appropriate preconditioning.

Mesh partitioning

To partition the mesh elements and nodes across a cluster, ISSM relies on METIS, a library for partitioning and ordering matrices and graphs. This partitioning scheme results in partitions that have equal numbers of elements on each cluster node which ensures a good computational load balance, as resulting stiffness matrices are well-partitioned. Figure 4.8 presents an example of partitioning of Pine Island Glacier, West Antarctica, into partitions that have the same number of elements.

Figure 4.8: Partitioning of the elements of Pine Island Glacier, West Antarctica, in 20 subdomains. Each domain is represented by a different color, the black lines are the elements edges. The model domain of this example is chosen to coincide with the geographic limits of Pine Island basin and limited to areas where observed velocities are available.

Model validation

Verification

We developed a suite of tools and a set of tests to automatically carry out nightly runs. These nightly runs consist in approximately a hundred test cases based on a variety of configurations that range from synthetic square ice shelves and ice sheets to simplified version of Pine Island and Nioghalvfjerdsfjorden glaciers. Tests include computation of velocity, temperature, thickness for steady-state or transient cases and about two thousand output fields are tested against archive files. An HTML report is automatically generated and sent by email every night. This ensures that no error is introduced in the software while developing new capacities.

Comparison with ISMIP-HOM benchmark

Model results of ISSM were compared against existing benchmarks such as ISMIP-HOM (Ice Sheet Model Intercomparison Project for Higher-Order Models) benchmark [Pattyn et al., 2008]. This benchmark targets validation for three dimensional ice flow models. Our results show an excellent agreement with the participating models, for all A to F tests. Here, we present some of these results, specifically for tests A, C and F. For these tests, we rely on regular triangular elements for the bidimensional mesh with 50 intervals in both horizontal directions. The bidimensional meshes are vertically extruded into 20 vertical layers to form prismatic elements. These tests are performed with the higher-order and the full-Stokes models.

Test A is a static experiment that involves an ice slab flowing over a sinusoidal bumpy-bed (see figure 4.9) with zero velocity at the base. We use the penalty method to impose periodic boundary conditions on the sides. The test is performed for six domain lengths (L) ranging from 5 km to 160 km. Figure 4.10 shows that our modeled surface velocity agrees well with the benchmark velocity computed by other software in the benchmark. Results from both higher-order and full-Stokes models are within the interval described by the benchmark models [Pattyn et al., 2008]. Experiment C of the ISMIP-HOM benchmark is similar to A, except for a flat bedrock and a bed that is not frozen as shown on figure 4.11. The basal friction coefficient is prescribed by a sinusoidal periodic law, so periodic boundary conditions are also applied here. Our results using both full-Stokes and higher-order models agree quite well with the benchmark results (see figure 4.12), and are also within the intervals described by other models. We note that the results for both models are very similar at all domain lengths.

The last ISMIP-HOM benchmark presented here is experiment F, which is a transient experiment for a slab of ice that flows over a sloping bed. The initial bedrock and surface are parallel and ice thickness is prescribed at 1,000 m. A gaussian bump is introduced at the center of the basal topography. As in the previous two experiments, penalties are applied to simulate the periodic boundary conditions. They are used for both velocity and thickness. The free surface and velocity evolve until a steady-state solution is reached. This experiment is run with basal sliding (slip bed) and without basal sliding (no-slip bed). Results are presented in figure 4.13 for both models. The steady-state surface found by both models and for both basal conditions agrees very well with the results of the ISMIP-HOM benchmarks. The steady-state surface velocity found for the non-slip bed is slightly higher than the results of other models. Steady-state surface velocity solved by ISSM is higher for the higher-order solution and lower for the full-Stokes solution by about 1m/yr compared to the case of the slip bed experiment. Only a handful of ice sheet models have performed these tests. In addition, this test exhibits challenges in terms of mass conservation, especially for finite element based treatments of mass transport [START_REF] Gagliardini | The ISMIP-HOM benchmark experiments performed using the Finite-Element code Elmer[END_REF].

In order to check whether our implementation conserves mass throughout the transient run, we computed the total mass of the ice slab at each step of the run. The mass is constant throughout the run, at a precision of 8 digits. This demonstrates that our numerical implementation efficiently conserves mass within single float precision. Chapter 5

A tiling method with applications to multi-modeling of ice sheet problems

Several ice flow models of increasing complexity have been developed to represent the dynamics of ice sheet systems. Simpler models require less computational time and ressources whereas sophisticated models capture the ice dynamics more accurately in all types of area. Hybrid models that combine several approximations have the potential to accurately capture ice dynamics while maintaining reasonable computational time.

In this chapter, we first describe a new coupling technique, the Tiling method, adapted specifically for our problem. We then apply it in order to combine the shelfy-stream approximation, the higher-order model and the full-Stokes equations.

Tiling method

The Tiling method is inspired by the Arlequin framework previously presented in section 2.6.3. It is based on a superposition zone where two different models are strongly blended, hence the name Tiling method. It mainly differs from the Arlequin framework in that there is only one mesh, which implies that this technique is not adapted to multiscale problems. The Tiling method does not require a domain decomposition and the domain is not "physically" split into two subdomains. Let us outline here the main ideas of this approach by considering a generic problem.

New formulation of a continuous problem

Let Ω be the model domain and V the set of kinematically admissible fields for Ω. We write a and l the bilinear and linear forms associated to the variational formulations and a Ω and l Ω their restriction to Ω. These two forms respectively represent the virtual work of internal and external forces. The generic problem to solve for is:

Find u ∈ V, ∀v ∈ V a Ω (u, v) = l Ω (v) (5.1)
Instead of solving this mono-model problem with the finite element method, the idea is to first reformulate this problem using a domain decomposition. We consider two subdomains of Ω denoted by Ω 1 and Ω 2 such that:

Ω = Ω 1 ∪ Ω 2 (5.2)
These two subdomains are overlapping in a transition zone. The overlapping region, denoted by Ω t , is called the transition zone or superposition zone. It is defined by: (see Fig. 5.1)

Ω t = Ω 1 ∩ Ω 2 = ∅ (5.3)
Let V 1 and V 2 be the sub-spaces of kinematically admissible fields for the two subdomains, such that

V = V 1 + V 2 .
If we use the bilinearity of a, the linearity of l, and if we decompose the test function, v, on the two spaces, the problem formulation using the domain decomposition becomes:

Find u = u 1 | Ω 1 + u 2 | Ω 2 ∈ (V t1 + V t2 ) , ∀v = v 1 | Ω 1 + v 2 | Ω 2 ∈ (V t1 + V t2 ) : a Ω 1 u 1 | Ω 1 , v 1 | Ω 1 model 1 + a Ω 2 u 2 | Ω 2 , v 2 | Ω 2 model 2 + a Ωt u 1 | Ω t , v 2 | Ω t + a Ωt u 2 | Ω t , v 1 | Ω t model coupling = l v 1 | Ω 1 model 1 + l v 2 | Ω 2 model 2
(5.7)

The coupling terms allow the coupling of u 1 and u 2 on the superposition zone.

This continuous problem is singular because of the existence of the volume superposition zone where the two subdomains are used simultaneously leading to an infinite redundant set of equations. We take advantage of the problem discretization to overcome this redundancy problem.

Domain discretization

In order to tackle the redundancy issue for finite element discretization of the continuous problem (5.7), our basic idea consists of reducing the superposition zone to only "one layer" of elements so that all nodes of Ω t are on the border of this zone ∂Ω t . Each node is associated with only one model except on ∂Ω t where one of the two models is constrained. Therefore, both u 1 and u 2 are never simultaneously computed on the same node in the mesh. Under these conditions, the discretized problem is regular.

To build this superposition zone, one possibility is to define a line (or a surface for three dimensional domains) that represents the border between the two subdomains and coincides with the elements edges. Then, find all the elements containing at least one of these points, keeping only the elements in one of the two subdomains (if all are kept, there will be two layers of elements, one inside the first subdomain and the other one inside the second subdomain). Figure 5.3 exhibits an example for a bidimensional domain decomposition.

Let us observe here that in this discrete setting, coupling with the tiling method may recall the one used by [START_REF] Belytschko | A coupled finite element -Element-free Galerkin method[END_REF] to couple finite and discrete elements. Let us also observe that the same formulation (5.7) has also been derived by Ben Dhia and Jamond to introduce a fine local model with a defect in a safe and global one; the transition zone being made of a one layer of coarse global elements (Ben Dhia and Jamond, 2008, private unpublished communication).

We hereafter adapt this method to the coupling of different ice flow models and show that for certain couplings, essential new technical tools must be designed in order for the method to be relevant. This is believed to be an original contribution, developed in this thesis.

We see that the final formulation is the sum of the two mono-model formulations and there are two additional coupling terms on the left hand side: a 1-2,Ωt (u 1 , v 2 ) and a 2-1,Ωt (u 2 , v 1 ). These terms, which must be defined for each model coupling, are non-zero only for the elements inside the blending superposition zone, where both u 1 and u 2 are simultaneously non-zero. Elements outside of this zone are treated as mono-model.

This method is easily parallelizable and requires limited modifications of an existing code, as the additional coupling terms must be taken into account only in the superposition zone.

Coupling shelfy-stream and higher-order models

In this section, we detail the coupling between the bidimensional shelfy-stream approximation (SSA) and the three dimensional higher-order model (HO) (see sections 3.3 and 3.2 for more details on these two models) using the Tiling method. These two models are derived from the full-Stokes equations and the horizontal and vertical velocities are computed separately.

We use the tiling method, as described in section 5.1 to couple these two models. SSA and HO solve the same system of decoupled equations with the same degrees of freedom, the latter with a three dimensional model and the former with a bidimensional depth-integrated model. Coupling them with the Tiling method is relatively straightforward, even if some adjustments must be made in the superposition zone.

Notations

Let u M and u P be the horizontal velocity associated respectively to the shelfy-stream approximation (SSA) and the higher-order (HO) model. The assumptions made for these two models lead to the following kinematics:

u M = u M (x, y) e x + v M (x, y) e y
u P = u P (x, y, z) e x + v P (x, y, z) e y (5.9) Therefore the strain rates ε associated to these velocities are:

ε (u M ) =             ∂u M ∂x 1 2 ∂u M ∂y + ∂v M ∂x 0 1 2 ∂u M ∂y + ∂v M ∂x ∂v M ∂y 0 0 0 - ∂u M ∂x - ∂v M ∂y             (5.10) ε (u P ) =             ∂u P ∂x 1
+ Γ b α 2 v P φ M y dΓ (5.20)
The finite element discretization and the mesh in particular must of course follow the conditions added in section 5.1 to avoid redundancy: the transition zone must be composed of only "one layer" of elements.

Treatment of ice viscosity

Two different models coexist in the transition zone and special care must be given to the treatment of ice viscosity, µ, in this area. Indeed, the viscosity depends on the effective strain rate and this strain rate is not constant with depth in the transition zone because of the contribution of HO. The strain rate used in the transition zone is:

ε (u) = ε (u M + u P ) = ε (u M ) + ε (u P ) (5.21)
ε (u) therefore varies with depth.

We use the most general viscosity law, the one used in HO, for both models in the superposition zone. As this viscosity varies throughout the thickness, SSA is no longer a bidimensional model and cannot be computed using a 2d surface. We must integrate the formulation over the volume and use a HO type viscosity, µ, instead of the depth-averaged viscosity, μ, usually used for SSA.

Coupling higher-order and full-Stokes models

In this section, we detail the coupling between the higher-order (HO) and full-Stokes (FS) models using the tiling method (see sections 3.2 and 3.1 for details on these models). The degrees of freedom solved for in these two models are different: in HO, the equations are decoupled, the horizontal components of velocity (u and v) and the vertical component (w) are solve sequentially, whereas in FS all the three components of velocity as well as the pressure are solved simultaneously.

Coupling FS to HO horizontal equations ensures the continuity of u and v but not w. To ensure the continuity of the vertical velocity, we must include a vertical velocity component in HO that is coupled to FS in the superposition zone. We use an iterative algorithm with an a priori estimate of the HO vertical velocity, which is updated at each iteration using the incompressibility equation so that vertical and horizontal velocities of HO remain consistent.

The pressure is a Lagrange multiplier and therefore does not have to be coupled between the two models.

Notations

Let v p and v s be the velocity fields associated respectively to the higher-order (HO) and the full-Stokes (FS) models, (u p , v p , w p ) and (u s , v s , w s ) their components in a cartesian coordinate system. The subscripts p and s refer respectively to HO and FS.

We use the superscript h to indicate that we refer only to the horizontal components of a vector and the superscript v to indicate that we only refer to the vertical component of this same vector. A tilde indicates that these components are extended with zeros in order to have a three dimensional vector. For example:

v p =   u p v p w p   , v h p = u p v p , ṽh p =   u p v p 0   , ṽv p =   0 0 w p   (5.22)

Domain decomposition

We want to find the velocity on a domain Ω where HO is applied everywhere except on a region where FS is used (see Fig. 5.5). We divide the domain Ω into two subdomains Ω p and Ω s corresponding to the areas where HO and FS are respectively applied. These two subdomains overlap on a superposition zone that is not empty: Ω ps = Ω p ∩ Ω s = ∅. We also define Γ p and Γ s the border of the transition zone such that: (5.23) As mentioned in section 5.1, the discretization of the domain must be such that the superposition zone is composed of only one layer of elements to avoid redundancy.

Γ p = ∂Ω p ∩ ∂Ω ps Γ s = ∂Ω s ∩ ∂Ω ps

Kinematically admissible fields

The solution of HO horizontal equations in Ω is composed of the two horizontal components of velocity. The space of kinematically admissible fields V h p,Ω is such that:

V h p,Ω = Φ ∈ H 1 (Ω) 2 , Φ = 0 on Γ u (5.24)
The solution of FS is the combination of velocity v s and pressure p s fields. The spaces of kinematically admissible fields for the velocity and pressure fields are respectively V s,Ω and P such that:

V s,Ω = Φ ∈ H 1 (Ω) 3 , Φ = 0 on Γ u and Φ • n = 0 on Γ b P = q ∈ L 2 (Ω) (5.25)
We are seeking a velocity field v = v p + v s and a pressure field p = p s . We use additional boundary conditions for the velocity fields:

v s = 0 on Γ s v p = 0 on Γ p (5.26)
The space of kinematically admissible fields for HO horizontal equations and FS for the hybrid model are therefore reduced to V h tp and V ts such that:

V h tp = Φ ∈ V h p,Ω , Φ = 0 on Γ p V ts = {Φ ∈ V s,Ω , Φ = 0 on Γ s } (5.27)
There is no boundary condition on the pressure p as it is a Lagrange multiplier which enforces the incompressibility in the full-Stokes equations. Furthermore, for this same reason, the pressure does not need to be coupled between HO and FS. The pressure is only computed on FS; for HO it is estimated afterwards using the hydrostatic assumption.

Hybrid problem

Before establishing the weak formulation of the hybrid problem, we write the mono-model formulations for HO and FS.

Higher-order weak formulation

The space of kinematically admissible horizontal velocity fields is V h p,Ω such that:

V h p,Ω = Φ ∈ H 1 (Ω) 2 , Φ = 0 on Γ u (5.28)
For any admissible field Φ p h ∈ V h p,Ω , we have:

a p v h p , Φ p h = l p Φ p h (5.29)
where a p and l p are the bilinear and linear forms associated to the virtual work of internal and external forces (see section 3.2).

Full-Stokes weak formulation

The space of kinematically admissible velocity and pressure fields are V s,Ω and P such that:

V s,Ω = Φ s ∈ H 1 (Ω) 3 , Φ = 0 on Γ u P = q ∈ L 2 (Ω)
(5.30) So for any admissible fields Ψ s = (Φ s , q) ∈ V s,Ω × P, we have:

a s (v s , Ψ s ) + b s (p, Ψ s ) = l s (Ψ s ) (5.31)
where a s and b s are two bilinear forms associated to the virtual work of internal forces and l s is a linear form associated to the virtual work of external forces (see section 3.1 for more details on these forms).

Hybrid model weak formulation

Following the description of the tiling method in section 5.1, the weak formulation of the hybrid HO/FS problem is:

Find v h p , v s , p ∈ V h tp × V ts × P , ∀ Φ h p , Φ s , q ∈ V h tp × V ts × P : a p,Ωp v h p | Ωp , Φ h p | Ωp HO + a s,Ωs v s | Ωs , Ψ s | Ωs + b s,Ωs p | Ωs , Ψ s | Ωs FS + a p-s,Ωps v p | Ωps , Ψ s | Ωps + a s-p,Ωps v h s | Ωps , Φ h p | Ωps coupling terms = l p,Ωp Φ h p | Ωp HO + l s,Ωs Ψ s | Ωs

FS

(5.32) with a i,Ω i the restriction of a i to Ω i , Ψ s = (Φ s , q), v p = ṽh p + ṽv p and v h s the horizontal components of v s . ṽv p is an a priori estimate of the vertical velocity component of HO. As this vertical component is not part of the solution we are seeking, terms including it must be on the right hand side of the equation: (5.33) In this example of coupling HO and FS, coupling terms are added both to the left hand side to couple ṽh p to FS terms and to the right hand side to ensure the continuity of ṽv p . if we write (φ x , φ y , φ z , q) the components of Ψ s .

a p,Ωp v h p | Ωp , Φ h p | Ωp HO + a s,Ωs v s | Ωs , Ψ s | Ωs + b s,

Details of the coupling terms

Vertical velocity

An additional constraint is to ensure that the incompressibility equation is verified everywhere in the domain, in the superposition zone Ω ps in particular.

To solve for the vertical velocity, we use the incompressibility equation, applied to the complete velocity v = v p + v s . The incompressibility equation is written in terms of the components of v p and v s :

∂ (w p + w s ) ∂z = - ∂ (u p + u s ) ∂x - ∂ (v p + v s ) ∂y (5.37)
The kinematically admissible fields needed to recover HO vertical velocity on Ω are V v p such that: 5.38) In the hybrid problem, this equation has already been used on Ω s and we only need to solve it on Ω p . As mentionned before, additional homogeneous Dirichlet conditions must be added to ensure the continuity of w = w p + w s . The space of kinematically admissible fields becomes:

V v p = Φ z ∈ H 1 (Ω), Φ z = 0 on Γ u (
V v tp = Φ z ∈ V v p , Φ z = 0 on Γ p (5.39)
The weak formulation associated to this equation is: 5.40) where Γ surf and Γ bed are the lower and upper boundary of the domain (ice/air interface at the surface and ice/bedrock at the bottom).

Ωp w p ∂Φ ∂z dΩ = Ωp ∂w s ∂z ΦdΩ + Ωp ∂ (u p + u s ) ∂x + ∂ (v p + v s ) ∂y ΦdΩ + Γ surf w p Φn z dΓ + Γ bed w p Φn z dΓ (
We use a non-penetration boundary condition at the ice/bedrock interface and do not constrain the vertical velocity on the upper surface, so this velocity is unknown and the integral on the upper surface will be part of the stiffness matrix (on the left hand side).

Once again this condition (no penetration of ice in the bedrock) applies to the complete basal velocity v p + v s :

w p (b) = (u p + u s ) ∂b ∂x + (v p + v s ) ∂b ∂y -w s (b) -Ṁb (5.41)
where Ṁb is the basal melting, positive when melting occurs at the bedrock.

If we add the previous boundary condition into the weak formulation, we get:

Ωp w p ∂ (Φ) ∂z dΩ - Γ surf w p Φn z dΓ = Ωp ∂ (u p + u s ) ∂x + ∂ (v p + v s ) ∂y ΦdΩ + Ωp ∂ (w s ) ∂z ΦdΩ + Γ bed (u p + u s ) ∂b ∂x + (v p + v s ) ∂b ∂y -w s -Ṁb Φn z dΓ (5.42)
The modifications introduced affect the vertical solution only on Ω ps , where FS velocity, v s , is not equal to zero. On the rest of HO domain, this formulation is identical to the usual vertical velocity weak formulation.

Iterative scheme

The following scheme (5.6) We start from a given velocity field and compute the associated ice viscosity. We compute v h p and v s using the hybrid model that combines FS and the horizontal equations of HO.

We then update w p using the incompressibility equation on Ω p and update the viscosity. These steps are repeated until the convergence criterion is reached.

This convergence scheme is similar to the one used for HO or SSA except that the vertical equation is computed at each iteration. This adds only a small cost to the computation as the incompressibility equation only involves one scalar field.

Coupling shelfy-stream and full-Stokes models

The shelfy-stream approximation is a bidimensional model very similar to the higher-order model, except for the terms ˙ xz and ˙ yz that are neglected. The model equations are then vertically integrated in order to have a bidimensional model. A method similar to the one described above to couple the higher-order and full-Stokes models can be applied to couple the full-Stokes and shelfy-stream models.

This coupling combines the difficulties encountered in the two previous couplings. In the superposition zone, we must use a viscosity µ that varies with depth even for the shelfystream approximation. The shelfy-stream horizontal equations and full-Stokes equations are solved simultaneously, which ensures continuity of the horizontal velocity. We use an apriori estimate of the vertical velocity for the shelfy-stream approximation, which is updated at each iteration of the viscosity, similar to the way coupling is done between higher-order and full-Stokes models.

Chapter summary

In this chapter we first presented a new method for coupling different mechanical models. This method is mainly derived from the Arlequin framework and Schwarz methods. The tiling method is less general than the Arlequin framework and therefore is not suitable for any combination of meshes and models. However, its implementation is not invasive and requires limited modifications in an existing software. This solution is adapted to the case of coupling different ice flow models, all deriving from the full-Stokes equations. We then explained the application of the tiling method for the coupling of the shelfy-stream, higherorder and full-Stokes models and detailed the adaptations that must be made in this specific case. In the next chapter, we present some results obtained using this technique. We run a diagnostic model to compute the velocity with these three models. Figure 6.4 shows the three components of velocity for each model.

As expected for an ice shelf, there is no difference between the results from the bidimensional and three dimensional models. Indeed ice shelves motion is dominated by lateral spreading and there is almost no vertical shear. Results from the hybrid model are identical to the two mono-models and the transition zone is impossible to detect visually.

Ice sheet

We perform the same experiment for a square ice sheet with a geometry similar to the ice shelf experiment (see figure 6.1), except that the bedrock is 20 m above the hydrostatic equilibrium. A small friction is applied on the base (α=20, see equation 2.28). The repartition of elements between shelfy-stream and higher-order in the hybrid model is identical to February 10, 2012 Shelfy-Stream Higher-Order Superposition zone Figure 6.3: Type of elements used in tiling on a horizontal cross section, the type of approximation used is the same for each element of a column the previous example (ice shelf experiment, see figure 6.3). In this example, homogeneous Dirichlet conditions are applied on all four edges of the domain; there is no ice front.

We compute the velocity for the three models: the two mono-models and the hybrid model. Figure 6.5 shows the three components of surface velocity for the three models and allows comparison of the results.

In this example, the velocities computed with the three models are very similar on this example. Ice motion on ice sheets is a combination of sliding and vertical shear. The friction coefficient applied is small so most of the motion is due to sliding and the computed velocities from the bidimensional and three dimensional models are almost equal. If the friction were higher and most of the displacement was caused by vertical shearing, only the three dimensional model would be able to include this effect. This experiment confirms the feasibility of our approach but does not enable us to conclude about its relevance.

Comparison of tiling and penalty method

In this section, we compare the results obtained for hybrid models with two different coupling technique: the tiling method and the penalty method. Areas where the shelfy-stream approximation and the higher-order models are used are almost identical for the two methods, however the tiling method uses a coupling in the volume while coupling takes place on a surface with the penalty method, therefore resulting in small differences at the transition between the two models. Figure 6.6 exhibits these differences in the type of elements used in the models.

We apply these two techniques on the previous experiments (square ice shelf and square ice Higher-order model 584s 1109s

Tiling Shelfy-stream/Higher-order 165s 397s

Penalties Shelfy-stream/Higher-order 257s 484s Table 6.1: Computational time required to solve the problem for each model

Coupling higher-order and full-Stokes

As in the previous section for the coupling between the shelfy-stream and higher-order models, we now study the coupling between the higher-order and full-Stokes models on square ice shelves and ice sheets and compare the results with the associated mono-models.

Ice shelf

We use the geometry described above on a square ice shelf with homogeneous Dirichlet conditions on three edges and an ice front where water pressure is applied on the fourth edge. We compare the results of models based on higher-order model, full-Stokes model and a hybrid model that combine the previous two models. For the hybrid model, we use the higher-order model on the left half of the ice shelf and full-Stokes model on the right half. The type of model used is the same throughout the thickness. Figure 6.9 shows the type of elements used on a horizontal cross-section of the domain.

Results of the velocity computed with these three models are presented in figure 6.10. The horizontal velocity is the same for the three models contrary to the vertical velocity. This difference comes from different treatment adopted in the two models. In full-Stokes, the vertical velocity is simply computed with the other component, whereas in the higher-order model, it is deduced using the incompressibility equation. However, one needs a boundary condition for this equation and no simple boundary condition can be applied (for ice sheets we use the non-penetration of ice in the bedrock); we therefore suppose that the ice shelf is in steady-state. This simple case in not in steady-state, therefore the vertical velocity in the higher-order is not accurate and differs from the one of full-Stokes. As the vertical velocities of higher-order and full-Stokes models are completely different, the one for the hybrid model is a combination of the two.

The jump in the pressure field at the transition zone is due to the fact that the pressure is not coupled between the models and its continuity is therefore not ensured. This field is computed at the same time as the velocity fields for full-Stokes, using no additional assumption, while it is calculated afterwards as a by-product using the hydrostatic assumption for the higher-order model. This assumption is globally correct but does not conduct to the exact same pressure than with full-Stokes; therefore, the transition between the two subdomains appears in the pressure field.

Higher-Order Full-Stokes Superposition zone Figure 6.9: Type of elements used in tiling on a horizontal cross section

Ice sheet

We now perform an experiment on a square ice sheet with similar geometry and boundary conditions. We compute the velocity and pressure for three models: a higher-order model, a full-Stokes model and a hybrid model that combines the latter two model with the tiling method (see figure 6.9 for the repartition of elements in the hybrid model).

Figure 6.11 shows the results for the three models. The three components of velocity are similar for the three models. As for the square ice shelf experiment, the boundary between the higher-order and full-Stokes model is visible on the hybrid model as the pressure is computed afterwards in the higher-order model using the hydrostatic assumption.

Convergence time Table 6.2 compares the computation time required to compute the ice velocity and pressure for the two previous experiments (square ice shelf and ice sheet) and the three models (higher-order, full-Stokes and hybrid models). Running an experiment with a full-Stokes model is almost ten times longer than doing the same experiment with a higher-order model in this case, even though there are only twice as many degrees of freedom. Using a hybrid model instead reduced this time by at least two in our examples. Additional experiments with varying number of elements should be performed before conclusions can be drawn about the computational gain of this method.

Flow over a rough bed

In this experiment, we test the efficiency of the Tiling method and its ability to combine models. The experiment consists of an ice slab flowing over a sloping bed, the ice surface is a sloping surface with same inclination as the bed. Near the ice front, a very rough bed is introduced, with bumps and hollows whose typical length scale is equal to one ice thickness 6.12. The average ice thickness is 1000 meters. This experiment is useful for testing the Tiling method, as the shelfy-stream, higher-order and full-Stokes models do not yield the same results. We assume the ice sheet to be isothermal so the ice hardness, B, does not vary; and we apply a friction coefficient at the ice/bedrock interface that is uniform. The boundary conditions are similar to the ones applied in the previous section (see figure 6.2). Non-homogeneous Dirichlet conditions are applied on three edges (Dirichlet conditions values are imposed to be equal to zero for the x and z directions and to 100 m/yr for the y direction), and water pressure is applied on the fourth edge.

Three ice flow mechanical models described in this thesis are employed to compute the ice velocity: shelfy-stream, higher-order and full-Stokes models. Ice velocity and pressure computed with these three models are presented in figure 6.13. This is not a realistic case as there would not be such a transition in bedrock roughness, but this experiment is a good test case to test the Tiling method, as the three models do not yield the same results.

The shelfy-stream and higher-order models are almost not affected by the rough bed (the results from these two models are very similar), while, in contrast, the velocity in the full-Stokes model near the ice front clearly influenced by the rough bed. The difference between the velocities of full-Stokes and the simplified models varies between 50 and 100 m/yr at the ice front and decreases upstream.

These results suggest that a full-Stokes model is required on and around the rough bed area to correctly capture the ice flow. We create three hybrid models that couple higher-order

Application to Pine Island Glacier

Pine Island is a major outlet glacier in West Antarctica that has been accelerating, thinning and retreating since the 1970's [START_REF] Rignot | Fast recession of a West Antarctic glacier[END_REF][START_REF] Schmeltz | Sensitivity of Pine Island Glacier, West Antarctica, to changes in ice-shelf and basal conditions: a model study[END_REF][START_REF] Rignot | Changes in West Antarctic ice stream dynamics observed with ALOS PALSAR data[END_REF]. It has been extensively studied [START_REF] Payne | Recent dramatic thinning of largest West Antarctic ice stream triggered by oceans[END_REF][START_REF] Joughin | Basal conditions for Pine Island and Twaites Glaciers, West Antarctica, determined using satellite and airborne data[END_REF][START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] to understand its basal properties. Results from [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] showed the importance of non-hydrostatic effects in the area above the grounding line and concluded that simple models like shelfy-stream or higher-order could not be used in this area as the rising bed exerts a backpressure on the flow and the bridging effect cannot be neglected (see figure 6.16). On the rest of the basin, however, both sophisticated and simpler models result in the same patterns of basal drag, suggesting that shelfy-stream and higher-order can be used on most of the domain for static studies. We are going to apply the Tiling method in order to use shelfy-stream or higher-order on the entire domain except in the vicinity of the grounding line, where FS is required. To reproduce the glacier configuration, we use satellite radar interferometry-derived surface velocities from Rignot et al. [2002b] (see figure 6.17-a), surface air temperatures from [START_REF] Giovinetto | Dependence of Antarctic surface mass balance on temperature, elevation, and distance to open ocean[END_REF], basal heat fluxes from [START_REF] Maule | Heat flux anomalies in Antarctica revealed by satellite magnetic data[END_REF], surface elevation of Antarctica from Bamber et al. [2009a] and ice thickness from [START_REF] Vaughan | New boundary conditions for the West Antarctic Ice Sheet: Subglacial topography beneath Pine Island Glacier[END_REF]. The model domain of Pine Island glacier is similar to the one used in [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] and corresponds to the limits of Pine Island Basin. We use an anisotropic triangular mesh whose metric is based on the hessian matrix of the velocity. This mesh is then vertically extruded into ten layers to create the three-dimensional mesh. The horizontal size of the elements varies between 500 meters in the vicinity of the grounding line and 5 km at the ice divide, the model contains about 150,000 prismatic elements (see mesh on figure 6.18).

The ice hardness depends mainly on temperature so we use a thermal model derived from the energy conservation to estimate B, assuming thermal steady-state:

ρc v • ∇T = k th ∆T + Φ (6.1)
where T is the ice temperature, c the ice heat capacity, k th the ice thermal conductivity, Φ the deformational heating and ∆ the Laplace operator. The surface temperature is imposed and the boundary condition at the ice/bedrock includes geothermal and frictional heat fluxes.

The basal drag coefficient and the ice shelf viscosity cannot be directly measured and are therefore inferred using a control method. We use a partial differential equation constrained optimization algorithm, which consists of a gradient minimization of a cost function that measures the misfit between observed and modeled horizontal surface velocities. [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] showed that a full-Stokes model is required in the vicinity of the grounding line, so we use the full-Stokes solution to infer the unknown parameters. We infer 1) the ice shelf viscosity on the ice shelves and 2) the basal drag coefficient at the ice/bedrock interface on grounded ice. We use the same parameters and boundary conditions for all experiments with lower order and hybrid models. We first run the shelfy-stream and higher-order models using the exact same parameters and then create hybrid models combining full-Stokes with the shelfy-stream or higher-order models and test their validity for this configuration.

Application of an inverse method enables a good fit between the model and observations (see figure 6.17 a and b). The average difference between the observed and modeled velocity is 20 m/yr. If we use the basal drag inferred from the full-Stokes model and apply it as a boundary condition for the ice/bedrock interface with a pure shelfy-stream or pure higher-order model, the modeled velocity deviates significantly from the observed velocity because they both neglect the bridging effect [START_REF] Kamb | Stress-gradient Coupling In Glacier Flow: 4 Effects Of The T Term[END_REF][START_REF] Van Der Veen | Force budget: I. theory and numerical methods[END_REF][START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF]. In the vicinity of the grounding line, the modeled velocity reaches respectively 5,200 m/yr and 5,000 m/yr with the shelfy-stream and higher-order models instead of 3,200 m/yr (see figure 6.17 b-d and figure 6.19).

Upstream of this area, however, the shelfy-stream and higher-order models are in good agreement with the results obtained from full-Stokes and with the observed velocities, with differences not exceeding 200 m/yr. Therefore, the basal drag field inverted by full-Stokes is usable as a boundary condition for the shelfy-stream or higher-order models up to 50 km upstream of the grounding line, which corroborates the conclusions of [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF].

We apply the tiling method on a hybrid model composed of full-Stokes model in the vicinity of the grounding line, from 50 km upstream to 35 km downstream of the grounding line on the fast moving ice stream, and the shelfy-stream of higher-order model on the rest of the domain (see figures 6.18 and 6.17 e-f). A forward model for this configuration shows an excellent fit with the observed velocities using full-Stokes basal drag as a boundary condition (see figure 6.19). The results deviate on average by 31 m/yr from full-Stokes for the hybrid higher-order/full-Stokes model and by 36 m/yr for the hybrid shelfy-stream/full-Stokes model, with a maximum difference around 200 m/yr on the fast ice stream, which represents less than 5% on this fast moving area. This difference is comparable to the difference between observed and modeled velocities (around 20 m/yr).

The computational time for a forward model of Pine Island Glacier is respectively 6.5s, 133s and and 805s for the shelfy-stream, higher-order and full-Stokes models on a 16-cpu machine, which represents a factor of 20 between the shelfy-stream and higher-order models and 6 between higher-order and full-Stokes. On these simulations, the area covered by the patch of full-Stokes elements in the hybrid model represents less than 10% of the glacier's domain but the number of full-Stokes elements is about 30% of the total number of elements. This is because we use an anisotropic mesh so the size of the elements greatly varies. The computational cost needed to run a hybrid higher-order/full-Stokes model reduced by half in this case (391s) and by more than three for a hybrid model shelfy-stream/full-Stokes (267s). It is therefore possible to include the accuracy of full-Stokes at a reasonable cost in critical areas where it is required, thus improving the model accuracy (see table 6.4).

Application to the Greenland Ice Sheet

The objective of this thesis was to develop a technique for simulations of ice sheet flow at the continental scale, using a full-Stokes model in critical areas and simpler models in areas where they are valid. In this section, we present a static model of the Greenland ice sheet using a hybrid model that combines the full-Stokes and higher-order models.

Datasets

This simulation is based on the dataset provided by the SeaRISE assessment (http:// websrv.cs.umt.edu/isis/index.php/SeaRISE_Assessment). The observations include mean annual surface temperature [START_REF] Fausto | Improving surface boundary conditions with focus on coupling snow densification and meltwater retention in large-scale ice-sheet models of Greenland[END_REF], precipitations [START_REF] Burgess | A spatially calibrated model of annual accumulation rate on the Greenland Ice Sheet (1958-2007[END_REF], basal heat flux [START_REF] Shapiro | Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica[END_REF], bedrock topography, ice thickness and surface elevation [START_REF] Bamber | A new ice thickness and bed data set for the Greenland ice sheet 1. Measurement, data reduction, and errors[END_REF]. Observed surface velocities are also provided [START_REF] Joughin | Greenland flow variability from ice-sheet-wide velocity mapping[END_REF] and filled with balance velocities (see e.g. [START_REF] Bamber | An analysis of balance velocities over the Greenland ice sheet and comparison with synthetic aperture radar interferometry[END_REF]) in areas where observations are not available.

Initialization

The objective of the SeaRISE assessment is to run evolutive models of the Greenland ice sheet for 500 years following different scenarios. We only present here the first step of these experiments, which is the model initialization. The initial geometry is defined by the present-day bedrock and surface topography. The extent of the domain is determined by the area where the ice is more than 5 meters deep.

A thermo-mechanical steady-state is run (see section 4.4.1 for more details) to ensure that velocity and temperature fields are consistent, and to initialize the ice hardness. The basal friction at the base of the ice sheet remains poorly understood and we rely on inverse methods to infer this parameter and match the observed surface velocity. We use the procedure described in [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] to perform this assimilation. This initialization is performed using a higher-order model that only includes grounded ice, as the spatial extent of ice shelves in Greenland is limited, to simplify the simulation.

We use an anisotropic bidimensional mesh, with element size ranging from 1 km along the coast to 25 km in the interior of the ice sheet. This mesh is extruded into 10 layers to create a model of about 500,000 elements.

Hybrid model

The discharge of Greenland ice is controlled by a few fast flowing ice streams, so we want to use the most sophisticated model in these critical areas. We construct a hybrid model based on the three-dimensional higher-order model, in which several patches of full-Stokes are introduced (see figure 6.20). Dirichlet conditions are applied for areas where the ice thickness is below 10 m. The velocity computed with this hybrid model is showed on figure 6.21, the transition between the different mechanical models is seamless.

Full-Stokes

Higher-order Superposition zone

Possible improvements

This simulation confirms the feasability and the relevance of our approach of constructing a hybrid model containing several patches of full-Stokes with different extent. This tool is integrated into ISSM so that only the diagnostic computation is affected by this change, and transient models can easily be performed using a hybrid model.

A detailed analysis of this result is beyond the scope of this thesis, but investigations should be made concerning the need for using a full-Stokes model, the repartition and size of full-Stokes patches and their evolution with time. Finally, investigating the differences between this hybrid model and pure higher-order or full-Stokes models would be of great interest.

Chapter summary

Simplification of the ice flow equations reduces the computational time required to run numerical models, which is useful for large-scale modeling, but it has been shown that these In this chapter we presented applications of the tiling method on both synthetic cases and real glaciers. In order to test the validity of the method, we first tested it on examples where both coupled models behave similarly. We then analyzed the case of a square ice sheet with a bump on the bed, which revealed that the full-Stokes model behaves differently than the other two models in such cases. We finally applied it on Pine Island Glacier, West Antarctica, as we know that a full-Stokes is needed in the vicinity of the grounding line for this glacier. All results from hybrid models, models that are combinations of different models, confirmed our expectations as their computational time is reduced compared to simple full-Stokes or higher-order models. This shows that ice flow models can be improved by including sophisticated models only in critical areas where such models are needed, therefore resulting in more efficient use of computational resources.

It is now essential to determine where each model is valid, and where sophisticated models are required. The next step is now to improve our knowledge and understanding of these ice flow models, their domain of validity in particular, for both idealized and real geometries. Criteria should be devised to determine what approximations to use. The regions where each approximation is employed should then automatically evolve with time based on these physical criteria.

Chapter 7

Conclusion

In this concluding chapter, we present the main developments achieved in this thesis, discuss the results and present the limitations of the present work. Suggestions for the application of the tiling method and for further developments of the numerical model are then presented before the general conclusion of this manuscript.

Main achievements and limitations

The main purpose of this work was to develop and apply a technique to couple different mechanical models of ice sheet flow of varying complexity. Before working on this problem, the first step was to build a parallel finite element platform inspired from Cielo, an inhouse finite element code, with the ISSM team. We then implemented all these models into the Ice Sheet System Model. The shelfy-stream approximation was already part of Cielo, whereas other models had to be included in this framework. Mathieu Morlighem and I started by adding the higher-order, full-Stokes models and shallow ice approximation into ISSM; the latter model is rarely used as it is a semi-analytical solution. Other models such as thermal, prognostic or balance thickness solutions were also added. A fair amount of time was devoted to implementing these solutions, creating a reliable platform and carrying out tests to validate the results (comparison of our results with analytical solutions and benchmarks experiments such as those presented in chapter 4). Finally, all the capabilities implemented in ISSM led us to participate in the SeaRISE assessment (http://websrv. cs.umt.edu/isis/index.php/SeaRISE_Assessment), an effort organized by the National Aeronautics and Space Administration (NASA) to estimate the upper bound of ice sheet contribution to sea-level rise by performing sensitivity experiments over the next 500 years and involving a dozen ice sheet models.

We then investigated several model coupling techniques. Our first approach was to employ simple techniques, such as the penalty method or Rigid Body motion. If this approach was adapted to couple the shelfy-stream and higher-order models, it was difficult to extend it to a full-Stokes model. We then considered iterative methods to combine full-Stokes to one of the two other simpler models, but the convergence of this method was very slow and required many iterations. Initialization with a good approximation of the solution helped reduce the computational time but remained inefficient for areas where results from a full-Stokes model were different from those from simpler models. We then considered using a more general method, such as the Arlequin framework. However, implementing such a general framework in an existing parallel software requires a lot of modifications and is not really necessary for our problem (no real need to include multi-scale effects, coupling of models relatively similar, ...). After looking at how to simplify this solution, we realized that a method such as the tiling method, that combines models by local superposition and linking was an appropriate solution to combine ice flow models.

After analyzing its conditions of applications, the tiling method was first applied to monomodel systems to ensure its validity. We then applied this techinque to the coupling of the shelfy-stream approximation with a higher-order model; this implementation was straightforward as the only specificity of this case is to ensure a compatible viscosity between the two models. To couple full-Stokes to one of the simpler models, an additional step was required: the vertical velocity of the models must be coupled even though it is computed separately for the simpler models. This was achieved by using an iterative scheme in which the vertical velocity of the simpler model is updated at each step of the convergence. After testing on idealized geometries, for cases where the two models behave either similarly or differently, the method was applied to real glaciers. One major application was the analysis of Pine Island Glacier, West Antarctica. Previous studies [START_REF] Morlighem | Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica[END_REF] showed that a full-Stokes model was needed in the vicinity of the grounding line but that simpler models could be used in the rest of the domain. A hybrid model that combines full-Stokes model around the grounding line, from 50 km upstream to 20 km downstream of the grounding line, and either shelfy-stream or higher-order model in the rest of the basin was used to compute the ice flow, which supported their results.

Combining ice flow models of varying complexity has the potential to improve models' physical accuracy while limiting their computational costs. Even though this technique improves results, it is clear that more work needs to be done in determining the domains of validity of each model, particularly when modeling real glaciers. Several studies showed that a full-Stokes model is required in the vicinity of the grounding line [START_REF] Nowicki | Conditions for a steady ice sheet-ice shelf junction[END_REF][START_REF] Durand | Full Stokes modeling of marine ice sheets: influence of the grid size[END_REF] but the extent of this area is poorly know. Furthermore, it is not clear if other areas should also be modeled with full-Stokes models, for example steep regions, as the assumptions made for the friction law might not be valid. Criteria should be defined to determine a priori the regions where each approximation can be employed, based on topography and observed velocity. A strategy similar to the one proposed by Ben Dhia et al. [2011] for atomic to continuum coupling could be used to define the size of the higher-order or full-Stokes patches.

Another limitation of this model is that the domain decomposition between the different models is static in our implementation, as these areas cannot evolve with time. Once criteria defining the domains of validity for each mechanical model are available, an important step will be to let the domain decomposition vary within a transient run, to better consider the evolution of the glacier.

The tiling method can be applied to many different systems and models and should not be restrained to ice sheet modeling. This method can be directly applied, for example, to beam models where more sophisticated models can be used close to the clamped edge of the beam, and simpler models can be used far from this edge.

Another important study achieved during my PhD but that is not developed in this thesis is the analysis of data consistency, though it is not included in this thesis. Ice flux divergence is an important quantity that determines the rate of thickness change of a glacier. However, large variations are detected when combining high-resolution surface velocity datasets with low-resolution measurements of ice thickness, which is inconsistent with what we know of glaciers surface mass balance, basal mass balance and thinning rate. We studied this problem on Nioghalvfjerdsfjorden, in North Greenland, and demonstrated that these inconsistencies, common to all glaciers, come from the interpolation of ice thickness measured on thickness tracks onto regular grids using schemes that do not conserve mass [START_REF] Seroussi | Ice flux divergence anomalies on 79north Glacier, Greenland[END_REF].

Perspectives

Providing realistic projections of ice sheet response in a warming climate and estimating their contribution to sea level rise is a major challenge for the glaciology community. Ice sheet numerical models are a fantastic and unique tools to better understand the physical processes controlling ice sheet flow and to project their behavior in a changing environment. Numerical modeling has gained a lot of attention in the last decades and the Ice Sheet System Model is part of this new generation of models, focused on developing high-resolution, higher-order models or data assimilation. Despite great progress, a lot of work remains.

Here I describe suggested improvements for ISSM.

The thermal model based on penalties sometimes takes a long time to converge and could be replaced by an enthalpy model capable of modeling both cold and temperate ice. Such a model also has the advantage of being easily coupled with a cryo-hydrology model. Most thermal models neglect the heat transport due to water flow through the ice [START_REF] Fowler | Flow of polythermal glaciers. 1. Model and preliminary-analysis[END_REF]. This heat transport is thought to influence the thermodynamics of glaciers, especially in Greenland where large amount of melted ice flows throughout the cryo-hydrologic system during the melt season [START_REF] Fountain | Water flow through temperate glaciers[END_REF]].

Phillips et al.

[2010] created a parameterization that allows incorporation of cryo-hydrologic warming into thermal ice sheet models. This cryo-hydrologic system can affect the thermal regimes of ice sheet within years to decade whereas the conduction of heat takes at least centuries to affect the thermal response of ice sheets. This mechanism has the potential to explain the rapid response of ice sheets to climate warming. This parameterization is based on the spatial average of temperature over representative areas and is controlled by the density and spacing of the cryo-hydrologic features. It does not require a complete description of the cryo-hydrologic system geometry or the complete resolution of the heat transfer system and is therefore computationally advantageous. [START_REF] Phillips | Cryo-hydrologic warming: A potential mechanism for rapid thermal response of ice sheets[END_REF] adds a term corresponding to heat exchange with the cryo-hydrologic system in the heat budget equation:

∂T ∂t + v • ∇T = k∆T + Φ ρc + k ρc T ch -T R 2 (7.1)
where T ch is the water temperature in the cryo-hydrologic system and R the average spacing between the cryo-hydrologic features. The temperature of the water in the cryo-hydrologic system is at the pressure melting point during the melt season. During winter, the water can refreeze if its enthalpy is reduced enough by the contact with the ice. This method was developed for a column model and has to be extended to three dimensional models.

Including such a model combined with an enthalpic treatment of the energy budget will allow to better represent the thermal aspect of our model.

Ice shelves and floating glacier termini are strongly affected by modifications in climate as demonstrated by the recent, spectacular collapses of several ice shelves in the Antarctic Peninsula. Ice shelves play a major role in the stability of their tributary glaciers as buttresses for the inland ice and control the discharge of ice to the ocean. The interactions between floating and grounded ice are still poorly understood, which limits the ability to predict the future evolution of major outlet glaciers and ice streams. A precise description of the grounding line is necessary to capture these effects however the exact location of the grounding line and its rate of advance or retreat is numerically challenging.

In ISSM, as in most ice sheet models, the position of the grounding line is determined by a hydrostatic equilibrium criterion. Recently, [START_REF] Nowicki | Conditions for a steady ice sheet-ice shelf junction[END_REF] developed a more rigorous grounding line model based on contact mechanics and the solution of the full-Stokes equations. Several full-Stokes models now use this technique to model grounding line dynamics [START_REF] Nowicki | Conditions for a steady ice sheet-ice shelf junction[END_REF][START_REF] Durand | Full Stokes modeling of marine ice sheets: influence of the grid size[END_REF], however they are all flow band models that rely on parameterization to account for the ice shelf buttressing that results principally from contact with embayment walls and ice rises, and they use a contact condition to determine the grounding line position. Furthermore, these models have been used only on idealized geometries and have not been applied to real glaciers or ice streams. Despite their limitations, these models showed that a very high spatial resolution is required in the vicinity of the grounding line and that ice shelves do not comply with the flotation condition in this region, as they are not in hydrostatic equilibrium in the first kilometers downstream of the grounding line. In the future, with improved computer resources, it might worth investigating these localized processes with a full Arlequin formulation.

Other models have been developed for planview simulations that include the two horizontal dimensions. Such models [START_REF] Goldberg | Grounding line movement and ice shelf buttressing in marine ice sheets[END_REF][START_REF] Joughin | Greenland flow variability from ice-sheet-wide velocity mapping[END_REF] are based on the Shallow Shelf approximation. [START_REF] Goldberg | Grounding line movement and ice shelf buttressing in marine ice sheets[END_REF] show that the extent and profile of an ice shelf affect the dynamics of the grounded ice sheet. Contrary to full-Stokes models, these models cannot employ a contact condition to determine the position of the grounding line, and therefore rely on the hydrostatic equilibrium condition, which has been demonstrated not to be accurate in the first kilometers downstream of the grounding line. Accurately modeling the grounding line position and its evolution requires three-dimensional full-Stokes models.

Another key feature missing in ISSM is the treatment of the ice front motion. Currently, ice fronts are fixed and no calving law is applied, which is equivalent to applying a freeflux boundary condition, and therefore having a stable fixed front. The position of ice shelves calving front is of major importance for the stability of the Antarctic Ice Sheet, as observed after the collapse of Larsen A or Larsen B [START_REF] Rott | Northern Larsen Ice Shelf, Antarctica: further retreat after collapse[END_REF][START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF]. The position of an ice front affects the behavior of the glaciers as it influences the stress field through lateral drag in an ice shelf embayment, the flow of their tributary glaciers, and therefore mass flux across the grounding line. A comprehensive calving law describing the advance and retreat of this front remains a challenge as a full calving model would require very small time steps and consideration of features with scales around 1 meter, such as iceberg formation. Therefore, calving parameterization is often employed rather than a calving model [START_REF] Amundson | A unifying framework for iceberg-calving models[END_REF]. Even though these parameterizations do not consider all the physical processes involved, including a moving ice front and a simple calving parameterization will help improve the accuracy of model results. 6.7 Comparison of the three components of velocity for a hybrid model combining the shelfy-stream approximation with the higher-order model with tiling method (left) and the penalty method (right) on a square ice shelf . . . . .

6.8

Comparison of the three components of velocity for a hybrid model combining the shelfy-stream approximation with the higher-order model with tiling method (left) and the penalty method (right) on a square ice sheet . . . . . 6.9 Type of elements used in tiling on a horizontal cross section . . . . . . . . .

  of glaciers and ice caps 1.1.1 The past 65 million years

Figure 2 .

 2 Figure 2.1: Description of ice sheet systems, http://en.wikipedia.org/wiki/Ice_shelf

Figure 2 . 2 :

 22 Figure 2.2: Description of the modeled ice sheet system (cross-section)

  y, b(x, y))+ w s -w b = 0(2.11)

Figure 6 .

 6 Figure 6. Pressure-adjusted temperature at the base for the control run (left) and the cold-mode run (right). Hatched area indicates where the ice is temperate. Values are in degrees Celsius and contour interval is 2 • C. The dashed line is the cold-temperate transition surface.

Figure 2 . 4 :

 24 Figure 2.4: Pressure adjusted temperature at the base of the Greenland ice sheet using an enthalpy formulation (left) and a cold ice model (right). Hatched areas indicate where the ice is temperate (ice at the pressure melting point). The dashed line is the cold-temperate transition surface. Courtesy of A. Aschwanden[Aschwanden et al., submitted]

Figure 8 .

 8 Figure 8. Basal melt rate for the control run (left) and the cold-mode run (right). Values are in millimeters per year. The dashed line is the cold-temperate transition surface.

Figure 2 . 5 :

 25 Figure 2.5: Basal melt rate of the Greenland Ice Sheet for the enthalpy formulation (left) and the cold ice model (right). Values are in millimiters. The dashed line is the cold-temperate transition surface. Courtesy of A. Aschwanden[Aschwanden et al., submitted]

Figure 2 . 6 :

 26 Figure 2.6: Temperature/Enthalpy (black line) and Temperature/Water fraction (blue line) functions. H s and H l are the enthalpy of pure ice and pure water at the pressure melting point T pmp and mark the limit between the cold ice, temperate ice and liquid water domains.

Figure 2 . 7 :

 27 Figure 2.7: Domain decomposition without (left) or with (right) superposition of the two subdomains

52February

  || = -α 2 u τ on Γ b (3.21) σn = p i n on Γ i (3.22) σn = 0 on Γ s (3.23)where τ b = (σn) || .

  59) Γw p w n x φ x dΓ = ωw p w n b x φ x JdΓ (3.60) Γw p w n y φ y dΓ = ωw p w n b y φ y JdΓ (3.61)

Figure 4 .Figure 4 . 1 :

 441 Figure 4.1: Meshes of Jakoshavn Isbrae, West Grenland and InSAR surface velocities overlaid on a MODIS mosaic of Greenland. a) Uniform mesh b) Anisotropic mesh whose metric is based on the velocity hessian matrix. Both meshes comprise around 1400 triangular elements. Adapted from Larour et al. [2012]

Figure 4 . 3 :

 43 Figure 4.3: Description of the melting-rate/temperature condition. T is the ice temperature, T pmp the ice pressure melting point and Ṁb the melting rate

T -T pmp Ṁb 0 κFigure 4 . 4 :

 044 Figure 4.4: Description of the melting-rate/temperature condition. T is the ice temperature, T pmp the ice pressure melting point, Ṁb the melting rate and κ the penalty coefficient

Figure 4

 4 Figure 4.7: Algorithm used for evolutive models

Figure 4 . 9 :

 49 Figure 4.9: Geometry of ISMIP-HOM benchmark test A, representing the upper and lower surfaces of the model domain.

Figure 4 .

 4 Figure 4.10: Results of ISMIP-HOM benchmark test A. Surface velocity (m/yr) across the bump at y=L/4 for different length scales, with L ranging from 5 to 160 km. Values computed with ISSM higher-order (NFS ISSM) and full-Stokes (FS ISSM) models are compared to the values found by other models (NFS and FS) in the original ISMIP-HOM benchmark [Pattyn et al., 2008].

Figure 4 .

 4 Figure 4.12: Results of ISMIP-HOM benchmark test C. Surface velocity (m/yr) across the bump at y=L/4 for different length scales, with L ranging from 5 to 160 km. Values computed with ISSM higher-order (NFS ISSM) and full-Stokes (FS ISSM) models are compared to the values found by other models (NFS and FS) in the original ISMIP-HOM benchmark[Pattyn et al., 2008].

Figure 5 .FebruaryFigure 5 . 5 :

 555 Figure 5.5 shows the decomposition of the domain on a horizontal cross section.Let ∂Ω p be the border of HO domain, ∂Ω s the border of FS domain, Γ u the part of the border of Ω where a Dirichlet condition is used, Γ σ the part of the border of Ω where a Neumann boundary condition is used (seeFig.5.5). For the sake of simplicity, we consider that all Dirichlet conditions are homogeneous Dirichlet conditions.

Figure 6 Figure 6 . 2 :

 662 Figure 6.1: Thickness on a flow line

Figure 6 . 4 :Figure 6 . 7 :

 6467 Figure 6.4: Comparison of the three components of surface velocity for the shelfy-stream approximation (left), the hybrid model (center) and the higher-order model for a square ice shelf

Figure 6 .

 6 Figure 6.10: Comparison of ice velocity and pressure at the base of a square ice shelf for a higher-order model (left), a hybrid model (center) and full-Stokes model (right)

Figure 6 .

 6 Figure 6.12: Geometry of the rough bed experiment

Figure 6 .

 6 Figure 6.16: Basal drag infered on Pine Island Glacier from several ice flow models: A)MacAyeal B)Blatter-Pattyn C)Full-Stokes from Morlighem et al. [2010]

Figure 6 .

 6 Figure 6.18: Mesh and type of elements used on Pine Island Glacier, West Antarctica for the hybrid model (top view). On the superposition zone, the elements are both FS and higher-order or FS and shelfy-stream.

Figure 6 .

 6 Figure 6.20: Modeled surface velocity of the Greenland ice sheet (m/yr)

Figure 6 .February

 6 Figure 6.21: Modeled surface velocity of the Greenland ice sheet (m/yr)

5. 2

 2 Example of solution in 1d . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Example of superposition zone in 2d . . . . . . . . . . . . . . . . . . . . . . 5.4 Decomposition of the domain between SSA and HO . . . . . . . . . . . . . 5.5 Domain decomposition between HO, FS and the transition zone . . . . . . . 5.6 Convergence algorithms used to combine HO and FS with the tiling method 6.1 Thickness on a flow line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Boundary conditions on a horizontal cross section (units are in km) . . . . . 6.3 Type of elements used in tiling on a horizontal cross section, the type of approximation used is the same for each element of a column . . . . . . . . 6.4 Comparison of the three components of surface velocity for the shelfy-stream approximation (left), the hybrid model (center) and the higher-order model for a square ice shelf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Comparison of the three components of surface velocity for the shelfy-stream approximation (left), the hybrid model (center) and the higher-order model for a square ice sheet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Type of elements used in tiling and penalty methods on a horizontal cross section, the type of approximation used is the same for all elements on the same column . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6. 10

 10 Comparison of ice velocity and pressure at the base of a square ice shelf for a higher-order model (left), a hybrid model (center) and full-Stokes model (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.11 Comparison of the ice velocity and pressure on the lower surface (ice/bedrock interface) of a square ice sheet for a higher-order model (left), a hybrid model (center) and a full-Stokes model (right) . . . . . . . . . . . . . . . . . . . . 6.12 Geometry of the rough bed experiment . . . . . . . . . . . . . . . . . . . . . 6.13 Comparison of ice velocity and pressure on the lower surface of the ice sheet for a shelfy-stream (left), a higher-order (center) and a full-Stokes model (right)134 February 10, 2012

  Energy balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Thermal problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 General ice flow problem . . . . . . . . . . . . . . . . . . . . . . . 36

	Chapter 2
	Physics of glaciers and numerical
	methods: a brief state of the art

Contents 2.1 Ice sheet system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.1 Components of ice sheet systems . . . . . . . . . . . . . . . . . . . 14 2.1.2 Three-dimensional geometry and notations . . . . . . . . . . . . . 15 2.1.3 Plane view geometry and notations . . . . . . . . . . . . . . . . . . 16 2.2 Mass balance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.1 Mass balance equation . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2.2 Boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.3 Depth-integrated mass balance equation . . . . . . . . . . . . . . . 18 2.3 Momentum balance . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 Momentum balance equation . . . . . . . . . . . . . . . . . . . . . 20 2.3.2 Ice behavior law . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.3 Mechanical boundary conditions . . . . . . . . . . . . . . . . . . . 22 2.3.4 Full-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.5 Simplified mechanical models . . . . . . . . . . . . . . . . . . . . . 25 2.4 28 2.4.1 Energy balance equation . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.2 Thermal boundary conditions . . . . . . . . . . . . . . . . . . . . . 29 2.4.3 2.6 Multi-model methods . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.6.1 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.6.2 Schwarz methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.6.3 Arlequin method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.6.4 Requirements for coupling ice flow models . . . . . . . . . . . . . . 43 2.7 Chapter summary . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Model validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.4 Shallow ice approximation . . . . . . . . . . . . . . . . . . . . . . Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.4.5 Model validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.5 Domains of validity . . . . . . . . . . . . . . . . . . . . . . . . . . 67

	Chapter 3
	Simplified mechanical models for
	ice sheet flow
	64
	3.4.1 Model assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
	3.4.2 Shallow ice approximation variational and local equations . . . . . 65
	3.4.3 Shallow ice approximation boundary conditions . . . . . . . . . . . 67
	3.4.4

Contents 3.1 Full-Stokes variational equations . . . . . . . . . . . . . . . . . . 46 3.2 Higher-order model . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2.1 Model description . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2.2 Higher-order model variational and local equations . . . . . . . . . 48 3.2.3 Higher-order model boundary conditions . . . . . . . . . . . . . . . 51 3.2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.2.5 Model validity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3 Shelfy-stream approximation . . . . . . . . . . . . . . . . . . . . 54 3.3.1 Model assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.2 Shelfy-stream approximation variational and local equations . . . . 55 3.3.3 Bidimensional model for horizontal velocity . . . . . . . . . . . . . 57 3.3.4 Shelfy-stream approximation bidimensional local equations . . . . 62 3.3.5 Shelfy-stream approximation boundary conditions . . . . . . . . . 62 3.3.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3.7 Simplification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.3.8

  Stokes weak formulation is:

	Ω	2µ	∂u ∂x	∂φ x ∂x	+ µ	∂u ∂y	+ µ	∂v ∂x	∂φ x ∂y	+ µ	∂u ∂z	+ µ	∂w ∂x	∂φ x ∂z	-p	∂φ x ∂x	dΩ
														+
	Γ Ω µ ∂u ∂y + µ ∂v ∂x ∂φ y ∂x + 2µ ∂v ∂y ∂φ y ∂y + µ	∂v ∂z	+ µ	∂w ∂y	∂φ y ∂z	-p	∂φ y ∂y	dΩ
														+			p i n y φ y dΓ
															Γ b		Γ i
	Ω	µ	∂u ∂z	+ µ	∂w ∂x	∂φ z ∂x	+ µ	∂v ∂z	+ µ	∂w ∂y	∂φ z ∂y	+ 2µ	∂w ∂z	∂φ z ∂z	-p	∂φ z ∂z	dΩ

b α 2 u τ φ x dΓ = Γw p w n x φ x dΓ + Γ i p i n x φ x dΓ α 2 v τ φ y dΓ = Γw p w n y φ y dΓ +

  Local equations for the horizontal velocityIf we use equation (3.15) for the pressure in the first two equations of the weak full-Stokes problem, equation (3.3), we obtain the horizontal part of the weak higher-order model:

	Ω	4µ	∂u ∂x	+ 2µ	∂v ∂y	∂φ x ∂x	+ µ	∂u ∂y	+ µ	∂v ∂x	∂φ x ∂y	+ µ	∂u ∂z	∂φ x ∂z	dΩ

15)

  We now integrate by parts the integrals over ω and Ω:

	Depth-integration of the volumic terms
	-	ω	∂ ∂x	4μH	∂u ∂x	+ 2μH	∂v ∂y	+	∂ ∂y	μH	∂u ∂y	+ μH	∂v ∂x	φ x dω
			+	∂ω	μH 4	∂u ∂x	+ 2	∂v ∂y	n x + μH	∂u ∂y	+	∂v ∂x	n y φ x dΓ +	Γ b	α 2 u τ φ x dΓ
					= -	Ω	ρg	∂s ∂x	φ x dΩ +	Γ i	p i n x φ x dΓ
	-	ω	ω ∂ ∂x + ∂ω = -Ω	s μH b μH µdz ∂u ∂y + μH 4 ∂v ∂u ∂x ∂x + Γ b α 2 u τ φ x dΓ = + 2 ∂v ∂y ∂φ x ∂x + ∂ ∂y 2μH Ω + ∂u ∂x ρg (s -z) ∂u ∂y + + 4μH ∂u ∂y + ∂v ∂x n x + μH 2 ∂u ∂x + 4 ∂v ∂y ρg ∂s ∂y φ y dΩ +	∂v ∂x ∂v ∂y n y φ y dΓ + φ y dω ∂φ x ∂y dω ∂φ x ∂x dΩ +	Γ b	α 2 v τ φ y dΓ
				ω		b	s	µdz				µ	∂u ∂y	+ µ	∂v ∂x	∂φ y ∂x	+ 2µ	∂u ∂x	+ 4µ	∂v ∂y	∂φ y ∂y	dω
								+	Γ b	α 2 v τ φ y dΓ =	Ω	ρg (s -z)	∂φ y ∂y	dΩ +	Γ i	p i n y φ y dΓ (3.37)
			We now introduce the depth-averaged viscosity as: Ω ρg ∂s ∂x φ x dΩ = ω s b ρg ∂s ∂x dz φ x dω	(3.41)
																μ =	1 H	b	s	µdz	(3.38)
			So the previous equations are:
				ω	μH				4	∂u ∂x	+ 2	∂v ∂y	∂φ x ∂x	+	∂u ∂y	+	∂v ∂x	∂φ x ∂y	dω
													+	Γ b	α 2 u τ φ x dΓ =	Ω	ρg (s -z)	∂φ x ∂x	dΩ +	Γ i	p i n x φ x dΓ
				ω	μH				∂u ∂y	+	∂v ∂x	∂φ y ∂x	+ 2	∂u ∂x	+ 4	∂v ∂y	∂φ y ∂y	dω
								+
			Γ February 10, 2012

Γw

p w n x φ x dΓ + Γ i p i n x φ x dΓ Γw p w n y φ y dΓ + Γw p w n x φ x dΓ + b α 2 v τ φ y dΓ = Ω ρg (s -z) ∂φ y ∂y dΩ + Γw p w n y φ y dΓ + Γ i p i n y φ y dΓ (3.39) ∂Ω ρg (s -z) φ x n x dΓ + Γw p w n x φ x dΓ + ∂Ω ρg (s -z) φ y n y dΓ + Γw p w n y φ y dΓ + Γ i p i n y φ y dΓ (3.40)

As the surface s only depends on x and y, it is possible to decompose the integral over the volume Ω into an integral over the thickness H and a second one over the mean section ω, which gives:

  Stokes system are:

	ω	∂ ∂x	4μH	∂u ∂x	+ 2μH	∂v ∂y	+	∂ ∂y	μH	∂u ∂y	+ μH	∂v ∂x	φ x dω
			-	∂ω		μH 4	∂u ∂x	+ 2	∂v ∂y	n x + μH	∂u ∂y	+	∂v ∂x	n y φ x dΓ -	Γ b	α 2 u τ φ x dΓ
			=	ω	ρgH	∂s ∂x	φ x dω -	γ i	1 2	ρ w gb 2 n x φ x dl
	ω	∂ ∂x	μH	∂u ∂y	+ μH	∂v ∂x	+	∂ ∂y	2μH	∂u ∂x	+ 4μH	∂v ∂y	φ y dω
			-	∂ω	μH		∂u ∂y	+	∂v ∂x	n x + μH 2	∂u ∂x	+ 4	∂v ∂y	n y φ y dΓ -	Γ b	α 2 v τ φ y dΓ
	=	ω	ρgH	∂s ∂y	φ y dω -	γ i	1 2	ρ

∂Ω ρg (s -z) φ x n x dΓ -Γw p w n x φ x dΓ + ∂Ω ρg (s -z) φ y n y dΓ -Γw p w n y φ y dΓ + w gb 2 n y φ y dl (3.46)

  -ρ w gb 2 n x φ x dlThe only part of the border where there is no Dirichlet boundary condition is the ice front γ i . On this part of the border, the boundary condition is therefore:

				∂ω ρgH 2 ∂ω μH 4 ∂u ∂x + 2 ∂v ∂y n x + μH ∂u ∂y + ∂v ∂x n y φ x dΓ = γ i 1 2 μH ∂u ∂y + ∂v ∂x n x + μH 2 ∂u ∂x + 4 ∂v ∂y n y φ y dΓ = γ i 1 2 ρgH 2 -ρ w gb 2 n y φ y dl
															(3.64)
					4μH	∂u ∂x	+ 2μH	∂v ∂y	+	∂ ∂y	μH	∂u ∂y	+ μH	∂v ∂x	φ x dω
			-		∂ω	μH 4	∂u ∂x	+ 2	∂v ∂y	n x + μH	∂u ∂y	+	∂v ∂x	n y φ x dΓ -	ω b	α 2 u τ φ x JdΓ
	=	ω	ρgH	∂s ∂x	φ x dω-     μH 4 γ i ρg ∂u H 2 2 ∂x + 2 φ x n x dΓ-∂v ∂y n x + μH ω ρgHφ x n b x JdΓ+ ∂u ∂y + ∂v ∂x 	γ i n y = 1 2 ρ w gb 2 n x φ x dl-1 ρgH 2 -ρ w gb 2 n x ωw p w n b x φ x JdΓ 2 (3.65)
	ω	∂ ∂x		μH  ∂u ∂y     μH + μH ∂u ∂v ∂x ∂y +	+ ∂v ∂x	∂ ∂y n x + μH 2 2μH ∂u ∂x	+ 4μH ∂u ∂x + 4	∂y ∂v ∂v ∂y	φ y dω n y = 1 2	ρgH 2 -ρ w gb 2 n y
			-		∂ω	μH		∂u ∂y	+	∂v ∂x	n x + μH 2	∂u ∂x	+ 4	∂v ∂y	n y φ y dΓ -	ω b	α 2 v τ φ y JdΓ
	= ρ ∂ ω ρgH ∂s ∂y φ y dω-γ i ρg H 2 2 φ y n y dΓ-γ i 1 2 ∂x 4μH ∂u ∂x + 2μH ∂v ∂y + ∂ ∂y μH ∂u ∂y + μH ∂v ∂x -α 2 u τ Jδ b = ρgH	∂s ∂x	-ρgHn b x J -p w n b x Jδ w
	∂ ∂x		μH	∂u ∂y	+ μH	∂v ∂x		+	∂ ∂y	2μH	∂u ∂x	+ 4μH	∂v ∂y	-α 2 v τ Jδ b = ρgH	∂s ∂y	-ρgHn b y J -p w n b y Jδ w
															(3.63)
	3.3.5 Determination of boundary conditions for the shelfy-stream ap-
			proximation

ω ρgHφ y n b y JdΓ+ w gb 2 n y φ y dlωw p w n b y φ y JdΓ (3.62) Now considering admissible horizontal and regular fields nul on the boundary of ω, we can easily deduce the local horizontal equations:

We now derive the boundary conditions from the weak formulation (eq. 3.62). If (u, v) is a solution of the problem, it verifies the local equations everywhere in ω. The boundary conditions are therefore:

So the boundary condition at the ice front is:

  table below summarizes the main properties of these models:

	Model	Dim. Unknowns Reference
	Full-Stokes (FS)	3d	4	Stokes [1845]
	Blatter-Pattyn (BP)	3d	2 + 1	Blatter [1995]; Pattyn [2003]
	Shallow shelf (SSA)	2d	2 + 1	MacAyeal [1989]; Morland and Zainuddin [1987]
	Shallow ice (SIA)	1d	1 + 1 + 1 Hutter [1983]

February 10, 2012

  Ωs p | Ωs , Ψ s | Ωs

				FS
	+ a p-s,Ωps	ṽh p | Ωps	, Ψ s | Ωps + a s-p,Ωps v h s | Ωps	, Φ h p | Ωps
				coupling terms
	= l p,Ωp Φ h p | Ωp	+ l s,Ωs Ψ s | Ωs	-a p-s,Ωps	ṽv p | Ωps	, Ψ s | Ωps
		HO		FS		coupling terms

  Ωs v s | Ωs , Ψ s | Ωs , b s,Ωs p | Ωs , Ψ s | Ωs , l p,Ωp Φ h p | Ωpand l s,Ωs Ψ s | Ωs are the usual terms of HO and FS formulations. The other terms are additional terms due to the coupling of the two models. , φ y , φ z , q) the components of Ψ s .

	a p-s,Ωps	ṽh p , Ψ s =	Ωps			2µ	∂u p ∂x		∂φ x ∂x	+ µ	∂u p ∂y	+ µ	∂v p ∂x	∂φ x ∂y	+ µ	∂u p ∂z	∂φ x ∂z	dΩ
					+	Ωps		µ	∂u p ∂y	+ µ	∂v p ∂x	∂φ y ∂x	+ 2µ	∂v p ∂y	∂φ y ∂y	+ µ	∂v p ∂z	∂φ z ∂z	dΩ
					+	Ωps		µ	∂u p ∂z		∂φ z ∂x	+ µ	∂v p ∂z	∂φ z ∂y	dΩ
					+	Ωps		-	∂u p ∂x	-	∂v p ∂y	qdΩ
																		(5.35)
	if we write (φ x a p-s,Ωps ṽv p , Ψ s =		Ωps			µ	∂w p ∂x		∂φ x ∂z	dΩ
					+	Ωps				µ	∂w p ∂y	∂φ z ∂z	dΩ
					+	Ωps				µ	∂w p ∂x	∂φ z ∂x	+ µ	∂w p ∂y	∂φ z ∂y	+ 2µ	∂w p ∂z	∂φ z ∂z	dΩ
	The terms a p,Ωp v h p | Ωp , a s,These terms are: , Φ h p | Ωp + Ωps		-	∂w p ∂z	qdΩ	(5.36)
	a s-p,Ωps v h s , Φ h p	=	Ωps	4µ	∂u s ∂x	+ 2µ		∂v s ∂y		∂φ x ∂x	+ µ	∂u s ∂y	+ µ	∂v s ∂x	∂φ x ∂y	+ µ	∂u s ∂z	∂φ x ∂z	dΩ
			+	Ωps	µ	∂u s ∂y	+ µ	∂v s ∂x		∂φ y ∂x	+ 4µ	∂v s ∂y	+ 2µ	∂u s ∂x	∂φ y ∂y	+ µ	∂v s ∂z	∂φ y ∂z	dΩ
																		(5.34)
	if we write (φ x , φ y ) the components of Φ h p .

  sumarizes how iterations are performed and how the a priori estimate of w p is updated. + a s (v s , Ψ s ) + b s (p, Ψ s )

	Initial velocity v = ṽh p + ṽv p + v s
	Compute v s + ṽh p
	Compute w p Viscosity and w p a p v p h , Φ h = l p Φ h p + l s (Φ s ) -a s p +ṽ v p ṽv p , Φ s convergence ∂w p ∂z = -div v s +ṽ h p p Solution v = v s +ṽ h a p-s ṽh p , Ψ s + a s-p v s h , Φ h p
	Figure 5.6: Convergence algorithms used to combine HO and FS with the tiling method

February 10, 2012

The pressure melting point is a term used in glaciology to define the temperature at which ice melts under a given pressure.February 10, 2012

3.6 Chapter summary . . . . . . . . . . . . . . . . . . . . . . . . . . .

Acknowledgements

I would like to take a moment to personally thank all the people who helped me and supported me during this journey.

ACKNOWLEDGMENTS

The authors thank J. Brown and M. Truffer for stimulating discussions. This work was supported by Swiss National Science Foundation grant no. 200021-107480, by NASA grant no. NNX09AJ38G, and by a grant of HPC resources from the Arctic Region Supercomputing Center as part of the Department of Defense High Performance Computing Modernization Program.

Stabilized finite elements

Limits of the finite element method for convective dominated flows [START_REF] Brooks | Streamline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] noticed that in convection dominated flows, results are sometimes affected by wiggles (spurious node to node oscillations). These anomalies are created by the finite element discretization of the problem described previously.

Several methods have been developed to tackle this problem. Here, we follow [START_REF] Hughes | A New Finite-element Formulation For Computational Fluid-dynamics Part. 8: The Galerkin Least-squares Method For Advectivediffusive Equations[END_REF]; [START_REF] Franca | Revisiting stabilized finite element methods for the advective-diffusive equation[END_REF] and add stabilization terms to the galerkin formulation. These mesh-dependent terms are consistent with the formulation and numerically stabilize the solution.

Chapter summary

In this chapter, we presented the numerical aspects and choices made in this thesis. We first described the spatial discretization adopted and the mesh used in this work as well as the time discretization. We then detailed the implementation of the mechanical problems, established the weak formulation of the full-Stokes model and the associated Galerkin approximation, as well the choices made in the implementation of the finite element model and the resolution of the non linearity of these equations. We also focused on the thermal problem to properly include the phase change when ice reaches pressure melting point; we described the formalism we adopted here and the choices made for the resolution of the non linearity of the thermal problem and in the stabilization of this convection dominated problem. We finally concluded this chapter by describing the main characteristics of the Ice Sheet System Model, the software in which all these formulations have been implemented. [Pattyn et al., 2008]. Analytical results derived with a first-order perturbation analysis of a flow down a uniformly inclined plane [START_REF] Gudmundsson | Transmission of basal variability to a glacier surface[END_REF] are also shown in black for the no-slip bed. The solution u is taken as the sum of the contributions of velocity from the two subdomains:

In order to have a continuous solution, u, we need to impose additional boundary conditions on the boundary of the superposition zone ∂Ω t . All points on ∂Ω t have homogeneous Dirichlet conditions for one of the two contributions (Fig. 5.2). For the sake of simplicity and without limiting the extent of the method, we detail here these additional Dirichlet conditions in the specific case where one subdomain is strictly included in the second one, i.e., no nodes of this first subdomain is on the border of the domain ∂Ω. Under these conditions, the additional Dirichlet conditions are:

The kinematically admissible fields are now:

Multi-model problem

The basic idea of the Tiling method consists of taking advantage of the reformulation of the finite element problem in order to introduce two different mechanical models for the two subdomains Ω 1 and Ω 2 . Schematically, let a 1,Ω 1 , a 2,Ω 2 and l 1,Ω 1 , l 2,Ω 2 be the bilinear and linear forms associated to the variational formulations of the two mechanical models in their respective domains. Let a 1-2,Ωt and a 2-1,Ωt be the transition blending bilinear forms in the transition domain Ω t . We define a weak multi-model formulation associated to problem (5.7) by:

(5.8)

February 10, 2012

Domain decomposition

We divide the domain Ω into two subdomains Ω M and Ω P corresponding to the areas where SSA and HO are respectively applied. These two subdomains overlap into a superposition zone that is not empty Ω M P = Ω M ∩ Ω P = ∅. We also write Γ M and Γ P the intersections of the border of the transition zone ∂Ω M P with respectively Ω M and Ω P :

(5.12)

Kinematically admissible fields

We are seeking a velocity field u that is solution of SSA on Ω M and solution of HO on Ω P . The velocity u is taken as the sum of the contributions of velocity on the two subdomains:

Homogeneous Dirichlet conditions must be added on the border of the transition zone as specified in section 5.1 to avoid discontinuity of the solution. In the case where one subdomain is stricty included into the second one:

106

February 10, 2012

COUPLING SHELFY-STREAM AND HIGHER-ORDER MODELS

If V M and V P are the space of kinematically admissible fields for SSA and HO for the same problem with a mono-model approach, the kinematically admissible fields for SSA and HO, with the hybrid approach are:

(5.15)

Hybrid formulation

We now establish the weak formulation of this problem. Using the tiling formulation (5.8), the problem becomes:

where a (u, Φ) and l (Φ) are:

(5.17)

The elements outside of the transition zone are not affected by the coupling, contrary to the elements in the superposition zone Ω M P because the coupling terms are equal to zero outside of the superposition zone. The bilinear forms associated with the coupling terms that must be added within the superposition zone are:

Chapter 6

Results

Contents

CHAPTER 6. RESULTS

In this chapter we present results of hybrid models coupled with the tiling method for both synthetic cases and real glaciers. We first apply this method to square ice shelves and ice sheets using hybrid models that combine first the shelfy-stream approximation with the higher-order model and then the higher-order and full-Stokes models. We compare these results to others that were computed with different coupling techniques, penalties and iterative methods in particular. We then study the case on an ice sheet flowing over a very rough bed. We finally apply this method to Pine Island Glacier, West Antarctica, where we know that full-Stokes is needed in the vicinity of the grounding line and to the Greenland ice sheet.

Synthetic experiments

The aim of this section is to validate the results obtained with the tiling method on simple geometries. To do so, we use square examples and compare the results obtained with the tiling method to the two mono-models used in the coupling and to other coupling methods when they are available. We start with the coupling between the shelfy-stream and higherorder models, then look at the results for coupling between higher-order and full-Stokes models. For these experiments, we use piecewise linear P 1 elements for the shelfy-stream and higher-order models, and P 1 + P 1 elements for full-Stokes model.

Geometry

We use a very simple geometry, which consists of a square ice sheet or ice shelf. The ice shelf is in hydrostatic equilibrium and the ice sheet base is set slightly above the equilibrium in order to have a non zero friction at the base. The ice density and viscosity are the same everywhere. Figure 6.1 represents the ice thickness on a flow line of the ice shelf.

We constrain the velocity (homogeneous Dirichlet conditions) on three edges and water pressure is applied on the fourth edge to simulate an ice front (see figure 6.2 for a representation of the boundary conditions on a horizontal cross section).

Coupling shelfy-stream and higher-order models

Ice shelf

We use the geometry described above and compare the results when using the shelfy-stream approximation on the whole domain, a higher-order model on the whole domain and when coupling these two models.

For the coupling, we use the shelfy-stream approximation on the right half of the ice shelf and the higher-order model on the left half of the domain. The type of model used is the same throughout the thickness. Figure 6.3 shows the type of element on a horizontal section of the ice shelf.

Computational time

We now investigate the computational time required to run these models. Table 6.1 summarizes the computational time needed to run the square ice shelf and square ice sheet experiments for the four models (shelfy-stream approximation, higher-order model, hybrid model using the tiling method and hybrid model using the penalty method).

Computing the three-dimensional velocity is at least twenty times longer than using a twodimensional model. This is expected as the number of degrees of freedom for this model is much higher: it is multiplied by the number of vertical layers of the model (if there are 10 vertical layers in the model, there will be 10 times more degrees of freedom to solve for in the three dimensional model).

The computational time required to run the hybrid models is situated between the computation needed to run a pure shelfy-stream model and the one of a pure higher-order model.

The two coupling methods lead to similar computational times. .15 shows the computed ice velocity and pressure for the three models. As expected, if the size of the patch of full-Stokes elements is too small, the flow is not slowed enough, and a very large patch that covers most of the domain is not necessary. With patches large enough to cover the influence area of the rough bed, the hybrid models are able to capture the effect of the rough bed, so that the modeled velocity is similar to the one modeled with a pure full-Stokes model.

The computational time of the hybrid models is longer than the time needed to compute the higher-order model but shorter than the time of full-Stokes model and depends on the size of the full-Stokes patch (see table 6.3 for a comparison of the computational time of all the models). The computed velocities all use basal-drag inverted for using full-Stokes model as a boundary condition at the ice/bedrock interface but different forward model: pure full-Stokes (red solid line), higherorder (green dotted line), shelfy-stream (loght blue dotted line), hybrid higher-order/full-Stokes (dark blue dahed line) and hybrid shelfy-stream/full-Stokes (pink dashed line). The vertical line represents the position of the grounding line in 1996 [Rignot et al., 2002a].

Model

Computational time Shelfy-stream approximation 6.