Boundary of spatial domain -

INTRODUCTION

This PhD is an interdisciplinary effort merging control engineering and glaciology. In this introductory chapter we answer several questions about the motivations behind our efforts, as well as establish the main topic.

In this chapter we:

Introduce the topic of this PhD Talk about the motivations behind our work Justify the use of control theory to solve a problem in geophysics

Introduce the main contributions of this work

Present the organization of this manuscript Figure 1.1: Control and Glaciology "fusing" together to solve a problem.

Main Topic

Through seasonal changes, our Earth experiences a diverse array of phenomena. Snowfall is one of those events, and depending on how much snow falls, the temperatures of the oceans and atmosphere, these fresh layers of snow can either melt or freeze into ice. Over time, if there is more freezing than melting, the layers can stack up. However, if more melting occurs, then layers of ice are lost. On earth, this process has been going on for millions of years, covering large areas of the planet with sheets of ice. Three major ice masses are present on the Earth's surface: the Arctic ice cap, and the Antarctic and Greenland ice-sheets. The last two holding more than 90% of the fresh water on Earth.

Once ice-sheets have formed, they are anything but a static environment. Ice is a viscous fluid [START_REF] Hutter | Theoretical Glaciology[END_REF],

meaning that bodies of ice will continuously deform due to forces acting upon them. To visualize this, imagine scooping honey out of a jar and pouring it on a piece of toast. The honey will gather on the surface of the bread and slowly start to spread outwards. Here, gravity is pulling on the honey causing it to deform. The exact same thing happens to ice-sheets, but at a much slower pace of course. However, this is not the only reason driving ice-sheets to move. Due to their immense weight, ice at the bottom of ice-sheets is subjected to high pressures lowering the melting point of ice thus causing it to melt. A layer of water forms, lubricating the rocks and sediments, causing the ice-sheets to slide. These processes are illustrated in Fig. 1.2. Another way to understand some of the dynamics of ice-sheets is to imagine a photographer taking pictures in a studio. If you have ever been in one of those studios, you might remember seeing large dark umbrellas with shiny and reflective undersides. Usually, photographers will place a light source behind an umbrella. This way, light will reflect all over the underside, then diffuse within the studio.

This gives the impression of having more natural lighting. The process of having the light spread out more evenly over space as opposed to coming from a more concentrated area, is called a diffusion process. Ice-sheets undergo a similar process as the bulk of the ice seeks to spread evenly over its underlying surface. However, this does not paint a full picture. Now, remember going to a park as a child and sliding down a slide. The process of you being on the top of the slide then finding yourself, unchanged, at the bottom of it, is an example of an advection process. You were transported from one part of space to another. Ice-sheets undergo a similar thing, where masses of ice will flow from higher places to lower ones around their domain causing visible changes to the surface topography of ice sheets.

The variation over time of the surface elevation h of grounded ice-sheets can be written in the form of a diffusion equation known as the ice thickness equation [START_REF] Huybrechts | The Antarctic ice sheet and environmental change: a three-dimensional modelling study= Der antarktische Eisschild und globale Umweltveränderungen: eine dreidimensionale Modellstudie[END_REF]. We assume that the bedrock elevation b = h -H varies much slower than H, thus its variation over time is omitted:

∂H ∂t = ∇ • D∇h + a (1.1)
where H is the ice thickness, a the mass balance (accounting for accumulation and ablation), and D the diffusion coefficient. This type of equation falls into the family of parabolic differential equations [START_REF] Evans | Partial differential equations[END_REF], the most famous of which is the heat equation. For ice-sheets, D is spatially varying in two dimensions and depends on the ice viscosity, ice thickness and bedrock elevation (and their slopes), temperature, and depending on the sliding law used, a sliding coefficient that we denote as A s [START_REF] Huybrechts | The Antarctic ice sheet and environmental change: a three-dimensional modelling study= Der antarktische Eisschild und globale Umweltveränderungen: eine dreidimensionale Modellstudie[END_REF][START_REF] Greve | Comparison of Numerical Schemes for the Solution of the Ice-Thickness Equation in a Dynamic/Thermodynamic Ice-Sheet Model[END_REF][START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF]. We can rewrite (1.1) to illustrate the effect of having a spatially varying diffusion coefficient: ∂h(t, x, y) ∂t = ∇D(x, y) • ∇h(t, x, y) + D(x, y)∇ 2 h(t, x, y) + a(x, y)

= ∂ ∂ x D ∂h ∂ x + ∂ ∂ y D ∂h ∂ y + a = ∂D ∂ x ∂h ∂ x + D ∂ 2 h ∂ x 2 + ∂D ∂ y ∂h ∂ y + D ∂ 2 h ∂ y 2 + a (1.2)
We see in (1.2) that because D varies through space, changes in h come from both D and its variations. This coefficient dictates how ice tries to evenly spread out over its surface (D and the second order partial derivatives of h), and how parts of the ice is transported along its slopes (the first order partial derivatives of D and h). To describe sliding, one can use a sliding law. In our case, we introduce this phenomenon through a Weertman sliding law [START_REF] Weertman | On the Sliding of Glaciers[END_REF]. This enables us to relate diffusion to basal sliding, which makes these two parameters dependant of each other. Detailed derivation of the ice thickness equation along with the exact laws and assumptions required to do so are given in the next chapter. The topic of this PhD is to estimate the values of these two coefficients for Antarctica.

In the next section, we justify why producing accurate estimates of D and A s is important.

Motivations Behind the Work

When it comes to understanding a phenomena, researchers collect measurements, set up environments, change control variables and observe their effects, and use physical laws to come up with mathematical models that can replicate reality. Simulating a system can offer a lot of insights into the hidden mechanisms of that system.

Let's go back to the honey pouring example. Maybe you have noticed, while pouring honey, the different ways it behaves before spreading over a surface. Next time you have some, try the following: scoop some of it and let it drip. Slowly, move your hand up or down and observe how the end of the honey thread can either fall in random patterns or coil up before collapsing. This kind of behavior is common to all viscous fluids, and was studied and simulated in [START_REF] Ribe | Coiling of viscous jets[END_REF][START_REF] Ribe | Stability of liquid rope coiling[END_REF]. It was found that different flow regimes exist, and that going from one regime to another depends on the fluid's viscosity, density, flow rate, the radius of the outlet, and the height of the fall. One might think that such behavior only occurs on smaller scales, but it is not true. On a much bigger scale, in the Earth's mantle, a similar process can occur [START_REF] Griffiths | Folding of Viscous Plumes Impinging On A Density Or Viscosity Interface[END_REF] when a plume of viscous matter descends through an environment with a different viscosity, causing it to bend or fold.

Scientists have investigated ice-sheet dynamics since the 1950s, but it took a few decades before numerical models emerged. These models depend on a set of parameters that determine the behavior of the model. For example, and without going into any technical details, let's consider a simple model that simulates how a specific ice-sheet evolves through time. Assume that our goal is to replicate available measurements of its surface topography. Let's also assume that the amount of snowfall and melt are known quantities, and that the boundary conditions are well defined, but that we do not know the geometry of the bedrock underneath the ice. This uncertainty in the bed's geometry will highly impact the results. It will not be possible to obtain the same results with a flat bed, an inclined one, or one with a more complex geometry, let alone match real measurements.

Though all ice-sheet models stem from the same constitutive relationships, they differ in their degrees of complexity, physical assumptions and simplifications, incorporation of physical boundary conditions, and numerical schemes. This causes their results to sometimes diverge even when run under similar environmental and numerical conditions. Great effort was spent to systematically compare ice-sheet models. Starting with EISMINT (European ice-sheet Modelling INitiative) [START_REF] Huybrechts | The EISMINT benchmarks for testing ice-sheet models[END_REF], a series of experiments was proposed to highlight the most efficient and reliable numerical techniques, as well as the pros and cons of the tested models. These experiments were set in three main categories. The first consists in fixing modelled processes and parameters and describing boundary conditions with as much detail as possible. The aim is to set up experiments with known analytical solutions, thus assessing the accuracy of the employed numerical schemes. Higher-order and full-Stokes ice-sheet models were intercompared in [START_REF] Pattyn | Benchmark experiments for higher-order and full-Stokes ice sheet models (ISMIP-HOM)[END_REF], while marine ice-sheet models were studied in [START_REF] Pattyn | Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP[END_REF][START_REF] Pattyn | Grounding-line migration in plan-view marine ice-sheet models: results of the ice2sea MISMIP3d intercomparison[END_REF] to evaluate the accuracy in simulating grounding line migration (the line dividing ice in contact with the bedrock and ocean).

The second type of experiments follows from the first, except that each model is allowed to incorporate additional processes deemed essential, as well as setting up its variables to preferred values.

Here, the goal is to assess the complexity of each individual model. Thermomechanical coupling was investigated on ten different ice-sheet models in [START_REF] Payne | Results from the EISMINT model intercomparison: The effects of thermomechanical coupling[END_REF]. And large-scale ice-sheet instabilities were studied on nine ice-sheet models in [START_REF] Calov | Results from the Ice-Sheet Model Intercomparison Project-Heinrich Event Intercomparison (ISMIP HEINO)[END_REF]. The focus in these kind of studies was on the nature of the phenomena rather than the numerical features of the models. In [START_REF] Asay-Davis | Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP+), ISOMIP v. 2 (ISOMIP+) and MIS-OMIP v. 1 (MISOMIP1)[END_REF] the set-up for three interrelated model intercomparison projects was proposed. This set-up tackled the coupling between ice-sheet and ocean models.

The third category of experiments focuses on real ice-sheets, like Antarctica or Greenland. One of the major questions tackled in this category is the contribution large losses in ice masses have on sea-level rise. In its fifth annual report [START_REF] Rajendra | Climate change 2014: synthesis report[END_REF], the IPCC (Intergovernmental Panel on Climate Change)

gave estimated ranges of global sea-level rise under different scenarios. By the year 2100, estimates range from 0.26 m to 0.82 m across all scenarios. However, observations show a global acceleration in mass losses for ice sheets. Because ice-sheet models cannot fully replicate these phenomena, estimates remain conservative.

The interactions between ice sheets, the ocean and the atmosphere, attract much attention because of the impact future predictions of sea level rise have on society. In [START_REF] Bindschadler | Ice-sheet model sensitivities to environmental forcing and their use in projecting future sea level (the SeaRISE project)[END_REF] the sensitivity of ten icesheet models to changes in surface mass balance, sub-ice-shelf melting and basal sliding was studied. The researchers found that the responses, of the tested ice-sheet models, vary linearly relative to the strength (measured as the change in the ice volume above flotation) of the imposed conditions. Researchers also found that the results of different forcings can be linearly combined.

This meant that studies employing simple relationships between climate and ice sheets proved to be adequate. However, refinement of the estimates is still needed. In [START_REF] Nowicki | Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project I: Antarctica[END_REF][START_REF] Nowicki | Insights into spatial sensitivities of ice mass response to environmental change from the SeaRISE ice sheet modeling project II: Greenland[END_REF], several ice-sheet models were subjected to atmospheric, oceanic and subglacial forcing scenarios, for both Antarctica and Greenland. This body of work belongs to a series of research focusing on the sensitivity of models to combinations of changes in the atmosphere, basal sliding underneath the grounded ice sheet, melting under floating ice. One of the researchers' findings, was that basal sliding plays an important role in the thinning of the Antarctic ice sheet. Regions of fast flow thin because of the increased discharge of ice, and even though thinning decreases towards the interior of Antarctica, it can still propagate over time. Improved understanding and estimation of the basal conditions beneath the ice will greatly benefit such research.

The previously mentioned studies have played an important role in qualitatively and quantitatively understanding the sensitivities of ice-sheet models, and in assessing why ice-sheet models can

diverge. Yet, another issue relative to ice sheet models is the impact the initialization has on the end results [START_REF] Seroussi | initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6[END_REF]. Often, two main approaches can be identified. The first method relies on long transient simulations from the last glacial maximum to present-day. This approach is successful in capturing climate history but the produced present-day ice-sheet profiles can greatly differ from the observations.

The second method uses data-assimilation techniques on present-day observations. While successful in reproducing these observations, results from this approach do not reproduce well past climate evolution. Predicting how much Greenland and Antarctica will contribute to sea level rise, ice-sheet models must be able to reproduce present-day observations. We note that in the second method, researchers often rely on techniques based on control engineering. 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] spatial distributions of basal drag were estimated by minimizing a cost function penalizing differences between observed and simulated horizontal surface velocities. And in [START_REF] Morlighem | Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model[END_REF], researchers inverted for basal friction through a variational approach introduced in [START_REF] Douglas | The basal stress distribution of Ice Stream E, Antarctica, inferred by control methods[END_REF]. In [START_REF] Arthern | Flow speed within the Antarctic ice sheet and its controls inferred from satellite observations[END_REF], inversion for the basal drag coefficient and ice stiffness parameter is carried using satellite observations of surface velocity, snow accumulation rate, and rate of change of surface elevation. Here we begin to see that control engineering is not foreign to glaciology.

All the previously mentioned works agree that an improved modelling of the physical processes and a better understanding numerical characteristics of ice-sheet models will be highly beneficial. It is clear that model initialization is a topic that requires extensive and multidisciplinary efforts. Our contributions fall into the inverse problem category because we use measurements of the surface elevation, bedrock topography, and mass balance for Antarctica to infer the basal sliding coefficient.

The choice of this coefficient comes from our problem formulation. We detail our choice and approach in the subsequent chapters. We also estimate the diffusion coefficient because it is closely related to basal sliding, and inferring one of these two parameters can lead to the other.

Control Theory as a Tool

Before answering how control theory can be used to solve a problem in geophysics, it is important to understand what is control engineering. Inherently, control is an interdisciplinary science. Numerous discussions with friends and colleagues have offered very different answers regarding the nature of this science. Yet, the most elegant and simple description can be found in Stuart Benett's book on the early history of control [START_REF] Bennett | A History of Control Engineering[END_REF] : "..engineers and mathematicians came together to create the control engineer". Of course, the resourcefulness and ingenuity of an engineer was the dominant element in the early days of control, mathematics gradually entered the scene as the field progressed and more complex problems emerged. To start putting things into perspective, let's explore how control evolved out of necessity by talking about windmills. The earliest mention of windmills dates back to the 9 th century. Vertical windmills were used in Persia to grind grain [START_REF] Bunch | The History of Science and Technology: A Browser's Guide to the Great Discoveries, Inventions, and the People Who Made Them, from the Dawn of Time to Today[END_REF]. Contrary to their modern cousins, they stood on vertical shafts with sails that rotated parallel to the wind's direction. Of course the design evolved over time, and in the 18 th century windmills were still a crucial element of human life. One major issue was their irregularity of motion as it heavily depends on the speed and direction of the wind. A solution was introduced by an English blacksmith in 1745. Edmund Lee [START_REF] Bennett | A History of Control Engineering[END_REF][START_REF] Bunch | The History of Science and Technology: A Browser's Guide to the Great Discoveries, Inventions, and the People Who Made Them, from the Dawn of Time to Today[END_REF] placed a smaller windmill, known as a fantail or fly, behind and perpendicular to the main sails to keep the windmill facing the incoming wind. However, this did not regulate the speed of rotation of the main sails. Traditionally, main sails were made of cloth. These had to be replaced with shutters that could be held closed using weights. By adjusting these weights, one could choose how strong the wind needs to be in order to blow them open, thus reducing the contact surface of the sails and slowing down their rotation. In 1772, Andrew Meikle [28], a Scottish millwright, improved on that design by replacing the weights with springs for each individual sail. The fantail reacted to the direction of the wind, the shutters to its intensity, and both devices fed back into the windmill the necessary adjustments.

Jumping to the early 20 th century, control started to take its modern form [START_REF] Bennett | A History of Control Engineering[END_REF]. The rush of industry made it necessary for more elaborate methods to increase the efficiency and reliability of production lines. All through the industrial spectrum, factories craved regulation. Production had to remain ongoing and stable; this required regulating the used devices. For example, motor speeds, frequencies, voltages, temperatures, flows, pressures had to be maintained at desired values [START_REF] Maxwell | On Governors[END_REF][START_REF] Routh | A treatise on the stability of a given state of motion: particularly steady motion[END_REF]. If the desired performance was not achieved, engineers would work on redesigning the process or machines used, and of course this led to a lot of trial and error. One major obstacle was the absence of a common language with which to communicate control engineering ideas and standards. Competing companies did innovate but many advancements were kept away from the public. However, and for whatever reason, in the 1930s companies did allow their employees to start publishing technical papers. This was of course a huge leap. And with the onset of World War II, the military relied on the knowledge and expertise of tech companies and their engineers to meet requirements. Notions of feedback, stability, tuning, frequency domain methods, all were introduced in that period [START_REF] Nyquist | Regeneration theory[END_REF][START_REF] Ashby | An introduction to cybernetics[END_REF].

Control theory continued to evolve, and the 1950s-1960s are considered to be the classical years.

During this period a lot of the modern language of control was introduced. Abstract concepts and theorems were established, and applications went beyond the industrial and military fields. In other words, control was not limited in its applications [START_REF] Bennett | A brief history of automatic control[END_REF]. If a concept can be reduced to its mathematical essence, then control can be a viable approach. We previously mentioned that control methods were already employed 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][START_REF] Morlighem | Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model[END_REF][START_REF] Douglas | The basal stress distribution of Ice Stream E, Antarctica, inferred by control methods[END_REF][START_REF] Arthern | Flow speed within the Antarctic ice sheet and its controls inferred from satellite observations[END_REF] to infer basal characteristics under ice sheets. We follow in this tradition, and utilize a branch of control engineering known as Lyapunov theory [START_REF] Lyapunov | The general problem of the stability of motion[END_REF][START_REF] Khalil | Nonlinear Systems[END_REF]. Control is a toolbox offering its user a wide array of strategies to tackle a problem and approach a domain. In our case, the domain is geophysics. With a background in control engineering, we offer a fresh set of eyes. Hopefully our perspective will prove to be useful.

Figure 1.4: Control in its early days, playing with windmills.

Main Contributions

The contributions of this PhD focus on estimating parameters relevant ice-sheet dynamics, in order to replicate measurements of the Antarctic surface topography.

Part I: Estimating Basal Sliding

The first part of this PhD focuses of estimating basal sliding. This two-dimensional (2D) parameter is spatially varying. The method proposed in this thesis is to derive an update or estimation formula for this coefficient, while still ensuring that the simulations of ice-sheet dynamics tend to the available measurements of surface topography. To test this strategy, we focus on a simplified version of the problem. Simplifications came in two parts. First, by reducing the spatial dimensions. This meant that 2D ice-sheets became one-dimensional (1D) flowlines. A flowline is simply the centerline of a glacier along which ice flows from the inner parts towards the edges. Second, by adopting one of the well known approximations in glaciology; the Shallow Ice Approximation (SIA). To understand what it implies, remember that ice is a slow moving viscous fluid. Its dynamics are described by the Stokes flow equations. One can exploit the fact that for large ice-sheets, the thickness to length ratio is small, meaning that ice-sheets are much longer than they are thick. This allows us to neglect some of the forces acting within the ice. This particular assumption works well for the interior of ice-sheets where the ice is much slower than near the edges. Doing so simplified the equations at hand and enabled us to obtain a new update law for the sliding coefficient. This law was tested

against case examples and real data.

Part II: Estimating Diffusion

The second part of this PhD emerged due to some of the difficulties encountered when trying to scale up the previous efforts back to two dimensions. This pushed us to shift our perspective and reformulate the problem. Instead of looking for basal sliding we look for the diffusion coefficient.

This perspective is rather new in glaciology as researchers rarely directly focus on this coefficient as it is of a more abstract nature, focusing instead on variables that have a more explicit physical meaning. However, this does not mean that estimating diffusion is useless. As will be apparent later, the diffusion coefficient can in fact be used to find basal sliding. Rather than coming up with a new update law, we used a pre-established one [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] and adapted it for Antarctica, because Antarctic data needs to be handled with more care. Our method is again tested against 1D and 2D case examples, as well as real data.

Manuscript Organization

This manuscript is divided into five chapters.

Chapter 2 presents some of the theoretical background involved in this PhD, for both glaciology and control. We start in Section 2.1 with concepts related to the physics of glaciers and their dynamics obtained through conservation laws of mass and momentum. We do not dwell on conservation of energy as thermodynamics is not explicitly involved in our work. The different forces acting within a glacier are explained and illustrated. Depending on dealing with ice-sheets or ice-shelves, certain forces dominate others. For ice sheets, the geometry can be exploited, allowing the omission of certain stresses in the balance of forces, which give rise to SIA. In Section 2.2 we expand the balance equations and define the boundary conditions that produce the ice thickness equation. Details of SIA and the modified ice thickness equation are given in Section 2.3. We then move to control.

Convergence of simulations towards measurements and stability of the equations representing the dynamics can be studied using Lyapunov theory. This is done in Section 2.4. We end the chapter with Section 2.5 where we expand on relevant works in glacilogy and control.

Chapter 3 focuses on the estimation of basal sliding in a 1D model. The available measurements of surface topography represent the steady-state of our system. Linearizing around this steady-state allows us to focus on deviations away from it, thus we obtain the linearized 1D flowline model. This process is described in sections 3.1 and 3.2. In Section 3.3 we conduct our Lyapunov analysis and define a new update law that guarantees convergence of the simulated data towards the measurements. In Section 3.4 we show how to iteratively implement this new update formula. In Section 3.5

we test this law in study cases for the linear model. In Section 3.6 we realize the same for the nonlinear model and compare the results with another method already used in ice-sheet models inversion for basal sliding. Then in Section 3.7 we validate our method on real data of the Antarctic ice-sheet.

In Chapter 4 we look for the diffusion coefficient in both 1D and 2D. In Section 4.1 we adapt the method of adaptive distributed parameter systems identification for our 2D problem. Following this method, the Lyapunov analysis is done in this section, we also introduce the notions of observer and update law for diffusion. Due to the stiff nature of the problem we introduce Tikhonov regularization in Section 4.2. This smoothing technique is used in the 1D study case. Numerical simulations and results for 1D are shown in Section 4.3. To gain better insights into the employed method, we study in this section the effects of changing some of the simulation variables. The same is done in Section 4.4 for the 2D study case. And in Section 4.5 we apply the method to real measurements of the Antarctic ice-sheet.

C h a p t e r 2

THEORETICAL BACKGROUND AND STATE OF THE ART

In this chapter, we offer some background in glaciology and control topics pertinent to this PhD. We start with a quick review of the balance equations and illustrate some of the forces acting within ice sheets. We then review ice sheet dynamics and derive the equations describing them. This is followed by a set of simplifications that will give us the main dynamical equation which is at the center of our work. Moving to control, we explore the Lyapunov methods and justify why they fit within our framework. We end the chapter with relevant advances in inverse problems done in glaciology and control.

In this chapter we: constructed our estimation methods. One of our motivations is to include all the technical background necessary for the understanding of this manuscript. Readers will not have to look too far away for explanations facilitating their understanding.

Introduce

Balance Equations

ice sheets are continent-wide masses of ice in contact with the ground. Their thickness can vary from a few hundred to a few thousand meters and they typically cover areas larger than 50, 000 km 2 .

Ice-shelves are outward extensions of ice sheets, float on ocean water, and are fed by the outflow of ice sheets. Due to gravity, these bodies of ice deform and flow as they tend to thin out and spread horizontally. Snowfall can accumulate in the higher parts while the lower areas experience melting and calving. In order to come up with the necessary dynamical equations, we offer here a brief review of relevant principles in fluid mechanics. For a more detailed description, we refer the reader to the following reference books [START_REF] White | Fluid Mechanics. 4th Ed[END_REF][START_REF]Chapter 3 -Equations of Fluid Dynamics[END_REF]. The following review is assembled using material from [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF].

This section is structured as follows:

1. Define the general balance law equation.

2. Define the mass balance (continuity) equation for an incompressible fluid.

3. Define the momentum balance equation.

To simplify further the readability of the chapter keep in mind the following: the mass balance and the boundary conditions (at the top and base of the ice sheet) will be used to find the ice thickness equation. The momentum balance equation and the shallow ice approximation will be used together to modify the ice thickness equation, allowing us to rewrite it in its diffusive form.

The General Balance Law Equation

Consider some material volume ω with a boundary surface ∂ω. This material could be anything as long as it has the same particles for all time t. For some point on ∂ω we consider a unit normal vector n and a velocity field vector v. The orthonormal basis vectors are x, y and z. Let X be the coordinates vector in material space. This setting is shown in Fig. 2.1. All vectors and matrices are denoted with a bold font.

Let G(ω, t) be a physical quantity on the entirety of the material volume ω, and let g(X, t) be its density. Assume that G is additive over subsets of ω. G can represent quantities like mass or momentum. Over time, changes in G can occur due to three processes:

• A flux F(∂ω, t), with density φ(X, t), of G through the boundary ∂ω.

• An internal production P(ω, t), with density p(X, t), of G within ω.

• An external supply S(ω, t), with density s(X, t), of G outside of ω. We assign a positive flux for a flow outward from ω. We can now define changes in G with time as the balance between F, P and S. We use the densities of these quantities to define this balance :

d dt ∫ ω g(X, t) dv = - ∮ ∂ω φ(X, t) • n da + ∫ ω p(X, t) dv + ∫ ω s(X, t) dv (2.1)
Where X is the position vector assigned to each particle in physical space, dv is a local volume change, da is a scalar surface element, and " • " is the dot product.

When it comes to dealing with control volumes, one can define the rate of change of some property of the fluid using Reynolds' transport theorem: Translating the above into an equation, we get:

d dt ∫ ω Γ(X, t) dv = ∫ ω ∂Γ ∂t dv + ∮ ∂ω Γv • n da (2.2)
Before going further, we recall the divergence theorem (Gauss's theorem or Green's theorem):

Theorem 2. The integral of the divergence of a vector field f over a volume ω is equal to the integral of the flux of f through the surface ∂ω:

∫ ω div(f) dv = ∮ ∂ω f • n da (2.3)
Note that the divergence operator for some vector field a is:

div(a) = ∂a x ∂ x + ∂a y ∂ y + ∂a z ∂z
This theorem enables us to replace surface integrals with volume integrals. More precisely, we will replace the first term in the right-hand side of (2.1) and the second term in the right-hand side of (2.2) by their equivalent volume integrals. Note that in (2.2) we replace Γ by g.

∮ ∂ω φ(X, t) • n da = ∫ ω div(φ) dv ∮ ∂ω gv • n da = ∫ ω div(gv) dv (2.4)
We can now use (2.2) and (2.4) to replace terms in (2.1), and gathering all terms under the same integral:

∫ ω ∂g ∂t + div(gv) + div(φ) -p -s dv = 0 (2.5)
As the above expression must remain true, the terms within the integral must be equal to zero. This gives the general balance equation in local form:

∂g ∂t = -div(gv) -div(φ) + p + s (2.6)
What (2.6) implies is that the changes in density g of a chosen material quantity G are balanced with the densities of production P and supply S within ω, the negative divergences of the flux density across the the boundary ∂ω, and the advective flux density in the normal direction n across ∂ω.

The Mass Balance Equation

Now, let us assume that the quantity G is the total mass M of the object with volume ω. We know that mass will not change and that it has a density ρ. We can write an expression stating just that:

d dt ∫ ω ρ dv = 0 (2.7)
If we compare (2.7) to the general balance equation in its integral form, we can identify all the variables in (2.1). We find that g = ρ, and that φ, p and s are zero. This enables us to rewrite the local balance equation thus obtaining the mass balance equation:

∂ ρ ∂t = -div(ρv) (2.8)
In the case of an incompressible fluid, which is a reasonable approximation for ice, the density does not change over time, and (2.8) becomes:

div(v) = 0 (2.9)
From this point onward we only deal with an incompressible fluid.

The Momentum Balance Equation

The quantity G can also represent the total momentum P. We know that momentum is mass multiplied by velocity, thus, its density is mass density times velocity, and P is the integral of that over the whole volume:

g = ρv P = ∫ ω ρv dv (2.10)
Before we proceed to obtain the momentum balance equation, consider Fig. 2.2 where the different forces acting on the material are shown. These forces can either be external forces f acting on volume elements of ω, or internal forces, usually termed stresses, e acting on the boundary surface ∂ω.

Newton's second law of motion states that the sum of all forces F on a body is equal to the rate of change of momentum P of that body. This translates to:

d dt P = F d dt ∫ ω ρv dv = ∮ ∂ω e da + ∫ ω f dv (2.11)
Similar to what we did earlier, we compare (2.11) to the general balance equation in its integral form (2.1) and identify the variables. We obtain g = ρv, p = 0, s = f, and φ = -e. Where the stress e (2.12)

The diagonal elements are known as the normal stresses, while the rest are the shear stresses.

They are all illustrated in Fig. 2.3 on each face of a cube. For each face, the element corresponding to e k k (k = x, y, or z) points outward and dictates the direction of the other two components. For example, take the bottom face, e zz has to point downwards, imposing a direction opposite to z, this implies that e xz must be opposite x, and e yz must be opposite y.

Having defined the terms for the momentum case, we can now plug them in (2.6) to obtain the momentum balance equation (note that ρ is set as constant):

∂(ρv) ∂t = -div(ρv v) + div(e) + f (2.13)
An equivalent form is:

ρ dv dt = div(e) + f (2.14)

Dynamics

Having defined the mass balance (2.9) and momentum balance (2.14) equations, we can now move on to finding the dynamical equation describing the evolution of ice thickness over time. In order to achieve that, in this subsection we will:

1. Develop the mass balance equation (2.9).

2. Define the ice sheet boundary conditions at the top and bottom.

3. Define the ice thickness equation.

Deriving the ice sheet Mass Balance Equation

Recalling (2.9) and developing it into its three components gives:

div(v) = ∂v x ∂ x + ∂v y ∂ y + ∂v z ∂z = 0 (2.15)
From now on, the ice surface (or free surface) elevation is denoted as h(t, x, y), the ice base (or bedrock) elevation as b(t, x, y), and the ice thickness as H(t, x, y). Of course, for a grounded ice sheet (where the ice is in touch with the bedrock), it is evident that H = hb. of the ice sheet (dotted red) has a surface elevation of h(t, x, y), its base in touch with the lithosphere (dotted green) is the bedrock elevation b(t, x, y), while the base submerged by the ocean (dotted black) is z (t, x, y), and the sea level (dotted blue) is z s (t, x, y).

In order to obtain the ice thickness equation, the volume flux Q must appear in (2.16). This flux has two components, one in the x direction, the other in y. Each component is the vertical integral of the horizontal velocity v x or v y . We thus must move the derivatives out of the integrals using the Leibniz integral rule, shown below for partial derivatives in x (the same principle applies for those in y):

∫ h(t,x,y) b(t,x,y) ∂ ∂ x f (t, x, y, z) dz = ∂ ∂ x ∫ h b f dz -f | z=h ∂h ∂ x + f | z=b ∂b ∂ x (2.17)
Using (2.17) for the first two terms in (2.16), and expanding the third gives us:

∂ ∂ x ∫ h b v x dz -v x | z=h ∂h ∂ x + v x | z=b ∂b ∂ x + ∂ ∂ y ∫ h b v y dz -v y | z=h ∂h ∂ y + v y | z=b ∂b ∂ y +v z | z=h -v z | z=b = 0 (2.18)
The integral terms in (2.18) constitute the two components of Q. The rest are to be rewritten using the boundary conditions that we define next.

Boundary Conditions

The boundary conditions for the ice top and ice base are given in this section.

Ice Top

At the surface of an ice sheet, the ice is in contact with the atmosphere. Let's consider a point on this free surface now denoted by F s , V to be the velocity of the free surface, U the velocity of the ice on this surface, and n the unit formal vector to the surface at this point. These vectors are shown in Fig. 2.5 (left). Remembering that z is an elevation above sea level and h(t, x, y) the height of the surface topography of the ice sheet, we can write the free surface as:

F s (t, X) = z -h(t, x, y) = 0 (2.19)
The time derivative of F s following the motion of the free surface is:

dF s dt = ∂F s ∂t + grad(F s ) • V = 0 (2.20)
We can also define a new variable, the ice volume flux through the free surface, also known as the surface mass balance, a s perpendicular to the free surface:

a s = (V -U) • n (2.21)
Notice that accumulation occurs when a s is positive, meaning that the normal component of V is greater than that of U. A loss of ice occurs in the opposite case. We can also remark that the unit normal vector n to the surface can be written as:

n = grad(F) | grad(F) | (2.
∂F s ∂t + grad(F s ) • U = -| grad(F s ) | a s (2.23)
Replacing F s by zh in (2.23) gives:

∂h ∂t + ∂h ∂ x U x + ∂h ∂ y U y -U z = | grad(F s ) | a s (2.24)

Ice Base

A similar reasoning can be made for the ice base. We now consider a surface F b , with the velocity vectors U and V, and the unit normal vector n as shown in Fig. 2.5 (right). Some changes to the equations of the free surface and mass balance are due: Notice that for a b , if the normal component of V is bigger than that of U melt occurs. The kinematic boundary condition at the base is:

F b (t, X) = b(x, y, t) -z = 0 a b = (U -V) • n (2.
∂b ∂t + ∂b ∂ x U x + ∂b ∂ y U y -U z = | grad(F b ) | a b (2.26)

Ice Thickness Equation

Having found the boundary conditions (2.24) and (2.26), we can now use them to replace the nonintegral terms in (2.18):

-v x | z=h ∂h ∂ x -v y | z=h ∂h ∂ y + v z | z=h = ∂h ∂t -| grad(F s ) | a s and v x | z=b ∂b ∂ x + v y | z=b ∂b ∂ y -v z | z=b = - ∂b ∂t + | grad(F b ) | a b (2.27)
Which gives us:

∂ ∂ x ∫ h b v x dz + ∂ ∂ y ∫ h b v y dz + ∂h ∂t -| grad(F s ) | a s - ∂b ∂t + | grad(F b ) | a b = 0 (2.28)
To make (2.28) more compact, we can use the fact that, for a grounded ice sheet, the ice thickness H is the difference between the surface elevation h and the bedrock elevation b (taken to be unchanging in the considered time scale). Also, we define the volume fluxes Q x and Q y to be the vertical integrals of the surface velocities v x and v y , respectively. The accumulation and ablation terms can be projected onto the vertical direction z, simplifying the equations | grad(F s ) | a s = a s and | grad(F b ) | a b = a b , thus we obtain the ice thickness equation:

∂H ∂t = -div(Q) + a s -a b (2.29)

Shallow Ice Approximation

In the previous sections, we derived the mass balance (2.9), momentum balance (2.14), and ice thickness (2.29) equations. These equations are descriptions of real physical phenomena. And though we would like to be as close as possible to reality, one can exploit certain physical characteristics in order to simplify them. One of these approximations is the Shallow Ice Approximation (SIA) [START_REF] Hutter | Theoretical Glaciology[END_REF]. This approximation exploits the fact that a small aspect ratio exists between the thickness of an ice sheet and its length, allowing us to neglect some of the components in the equations we have seen. In a big part this assumption is true as ice sheets tend to be much longer than they are thick. But the SIA fails in marginal regions [START_REF] Baral | Asymptotic Theories of Large-Scale Motion, Temperature, and Moisture Distribution in Land-Based Polythermal Ice Sheets: A Critical Review and New Developments[END_REF] especially when acting forces become of equal importance [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]. We did however base a good deal of our work on SIA as we are mostly interested in the grounded ice sheet where this approximation can remain valid. The simplifications introduced by SIA offer another advantage, as the balance equations can take practical forms that will allow us to derive our estimation law in Chapter 3.

In the previous sections, we have illustrated some of the dynamics of the movement of ice sheets.

We know by now that ice will move and deform due to its own weight, causing it to flow. In this section we need to modify the ice thickness equation (2.29) under SIA. The volume flux Q can be replaced with horizontal velocities. In turn, these velocities can be deduced from the strain rate which is introduced through Glen's flow law [START_REF] Glen | The creep of polycrystalline ice[END_REF][START_REF] Nye | The distribution of stress and velocity in glaciers and ice-sheets[END_REF]:

= A σ n-1 e e d (2.30)
where is the strain rate tensor, A a temperature dependant rate factor, σ e the effective stress, e d the Cauchy stress deviator, and n a constant. In our case, the normal stress deviators e d x x , e d yy and e d zz are neglected.

Before utilizing (2.30), we start by rewriting the momentum balance equation (2.14) in component form, while neglecting the acceleration term and setting f = -ρge z :

∂e x x ∂ x + ∂e xy ∂ y + ∂e xz ∂z = 0 ∂e xy ∂ x + ∂e yy ∂ y + ∂e yz ∂z = 0 ∂e xz ∂ x + ∂e yz ∂ y + ∂e zz ∂z = ρg (2.31)
We now consider that throughout the ice sheet, the gradients of the shear stresses e xz and e yz are negligible compared to that of the normal stress e zz . Also, because the slopes of the ice surface and bed are considered to be small (due to the small aspect ratio), the gradient of the shear stress e x y is negligible too. Setting all normal stresses equal to the negative of the pressure p greatly simplifies (2.31):

∂e xz ∂z = - ∂e x x ∂ x = ∂p ∂ x ∂e yz ∂z = - ∂e yy ∂ y = ∂p ∂ y ∂e zz ∂z = ∂p ∂ x = ρg (2.32)
The third equation in (2.32) can be integrated from z to h, giving p = ρg(hz). We can now replace p in the first two equations of (2.32):

∂e xz ∂z = ρg ∂h ∂ x ∂e yz ∂z = ρg ∂h ∂ y (2.33)
Due to the small aspect ratio, the surface derivatives are small, which sets the unit normal vector as n = e z . This can cause the stress-free boundary condition for the free surface to become:

p | z=h = 0 , e xz | z=h = 0 , e yz | z=h = 0 (2.34)
Integrating (2.33) and using (2.34):

e xz = -ρg(h -z) ∂h ∂ x e yz = -ρg(h -z) ∂h ∂ y (2.35)
We can now return to Glen's flow law (2.30). By definition, the strain rate tensor can be written in term of velocities:

=        x x xy xz yx yy yz zx zy zz        =         ∂v x ∂ x 1 2 ( ∂v x ∂ y + ∂v y ∂ x ) 1 2 ( ∂v x ∂z + ∂v z ∂ x ) 1 2 ( ∂v y ∂ x + ∂v z ∂ y ) ∂v y ∂ y 1 2 ( ∂v y ∂z + ∂v z ∂ y ) 1 2 ( ∂v z ∂ x + ∂v x ∂z ) 1 2 ( ∂v z ∂ y + ∂v y ∂z ) ∂v z ∂z         (2.36)
As we are only interested in the x-z and y-z components, we select xz and yz from (2.36) and calculate them using (2.30):

xz = 1 2 ∂v x ∂z + ∂v z ∂ x = A σ n-1 e e xz yz = 1 2 ∂v y ∂z + ∂v z ∂ y = A σ n-1 e e yz (2.37)
The two shear stresses e xz and e yz were calculated just above in (2.35). The only missing term is the effective stress σ e which can be found using the trace of the deviatoric stress tensor e d :

σ e = 1
σ e = e 2 xz + e 2 yz = ρg(h -z) | grad(h) | (2.39)
We can now substitute (2.35) and (2.39) in (2.37) to obtain:

1 2 ∂v x ∂z + ∂v z ∂ x = -A ρ n g n (h -z) n | grad(h) | n-1 ∂h ∂ x 1 2 ∂v y ∂z + ∂v z ∂ y = -A ρ n g n (h -z) n | grad(h) | n-1 ∂h ∂ y (2.40)
Due to the property of the aspect ratio, the horizontal derivatives of vertical velocities in (2.40) can be neglected. However, the derivatives of vertical velocities can be integrated from the base z = b to a height z:

v x = v b x -2ρ n g n | grad(h) | n-1 ∂h ∂ x ∫ z b A (h -ζ) n dζ v y = v b y -2ρ n g n | grad(h) | n-1 ∂h ∂ y ∫ z b A (h -ζ) n dζ
where v b x and v b y are the velocities at the base.

(2.41)

SIA does not usually include sliding but due to our interest in this phenomena, we augment SIA with a Weertman sliding law [START_REF] Weertman | On the Sliding of Glaciers[END_REF][START_REF] Pattyn | Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0)[END_REF] of the form:

v b = A s | τ b | m-1 τ b (2.42)
Three new terms are introduced in (2.42): the basal sliding coefficient A s , the basal drag or basal shear stress τ b , and a constant basal sliding exponent m. Basal sliding describes how the ice slides over the underlying bedrock, while τ b represents the friction arising from such movement. Parallel to the basal drag is the driving stress τ d , which can be computed by evaluating (2.35) at the base z = b, giving:

τ d = -ρgHgrad(h) (2.43)
It was noticed in [START_REF] Morlighem | Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model[END_REF] that for a large portion of the Antarctic ice sheet the driving stress is balanced by the basal drag. This is why, for our calculations we set τ b = τ d , and replace it in (2.42):

v b x = -A s ρ m g m H m | grad(h) | m-1 ∂h ∂ x v b y = -A s ρ m g m H m | grad(h) | m-1 ∂h ∂ y (2.44)
Using (2.41) and (2.44), we can now calculate the volume flux Q:

Q = Q x Q y = ∫ h b v x dz ∫ h b v y dz Q x = ∫ h b -A s (ρgH) m | grad(h) | m-1 ∂h ∂ x -2(ρg) n | grad(h) | n-1 ∂h ∂ x ∫ z b A (h -ζ) n dζ dz Q y = ∫ h b -A s (ρgH) m | grad(h) | m-1 ∂h ∂ y -2(ρg) n | grad(h) | n-1 ∂h ∂ y ∫ z b A (h -ζ) n dζ dz (2.45)
During our work, we decided to set constant some of the variables above. The exponent m in (2.42) dictates how plastic sliding can be. In [START_REF] Joughin | Regularized Coulomb Friction Laws for Ice Sheet Sliding: Application to Pine Island Glacier, Antarctica[END_REF], friction laws and the effects of varying m are studied. This exponent can vary from m = 1 for a linear viscous flow, to m → ∞ for a more plastic flow. For our work, we set m = 2 [START_REF] Pattyn | Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0)[END_REF]. As for the exponent in Glen's law (2.30), it is often set to n = 3. Additionally, the temperature dependant factor A usually varies with depth, however we set it constant. We can now rewrite (2.45), keeping in mind that within the main integral the only term dependant on z is the integral term:

Q x = -A s ρ 2 g 2 H 3 | grad(h) | ∂h ∂ x -2A ρ 3 g 3 | grad(h) | 2 ∂h ∂ x ∫ h b ∫ z b (h -ζ) 3 dζ dz Q y = -A s ρ 2 g 2 H 3 | grad(h) | ∂h ∂ y -2A ρ 3 g 3 | grad(h) | 2 ∂h ∂ y ∫ h b ∫ z b (h -ζ) 3 dζ dz (2.46)
Integrating the double integral gives:

∫ h b ∫ z b (h -ζ) 3 dζ dz = 1 5 H 5 (2.47)
Replacing (2.47) in (2.46):

Q x = -A s ρ 2 g 2 H 3 | grad(h) | ∂h ∂ x - 2 5 A ρ 3 g 3 H 5 | grad(h) | 2 ∂h ∂ x Q y = -A s ρ 2 g 2 H 3 | grad(h) | ∂h ∂ y - 2 5 A ρ 3 g 3 H 5 | grad(h) | 2 ∂h ∂ y (2.48)
The expressions for the volume fluxes in the x and y directions are now complete and can be substituted into the ice thickness equation (2.29):

∂H ∂t = ∂ ∂ x (D ∂h ∂ x ) + ∂ ∂ y (D ∂h ∂ y ) + a s -a b with D = A s ρ 2 g 2 H 3 | grad(h) | + 2 5 A ρ 3 g 3 H 5 | grad(h) | 2 (2.49)
The above equation is the ice-thickness equation in its diffusive form under SIA. We will adopt it throughout our work. The variable D is the diffusion coefficient that will be the focus of Chapter 4.

Notice how this variable depends on basal sliding A s .

Before moving on to the next section, we recap below the set of equations and assumptions that have enabled us to obtain (2.49):

• We began with the ice thickness equation as it was defined at the end of Section 2.2:

∂H ∂t = -div(Q) + a s -a b
• The volume flux Q was defined as in (2.45):

Q = Q x Q y
• We adopted the Shallow Ice Approximation which enabled us to neglect some of the stresses in the balance of forces.

• The components of the surface velocity were given as in (2.41):

v x = v b x -2ρ n g n | grad(h) | n-1 ∂h ∂ x ∫ z b A (h -ζ) n dζ v y = v b y -2ρ n g n | grad(h) | n-1 ∂h ∂ y ∫ z b A (h -ζ) n dζ
• We used a Weertman sliding law to augment SIA, thus the components of the velocity at the base were given as in (2.42):

v b = A s | τ b | m-1 τ b
• We assumed that for large portions of the Antarctic ice sheet the basal drag τ b is balanced by the driving stress τ d given by (2.43):

τ d = -ρgHgrad(h)
• The exponent m describing the plasticity of sliding was set equal to 2. Glen's law exponent n was set equal to 3. And the temperature dependant factor A , which can vary with depth, was considered to be constant.

• Combining all the above leads to the final diffusive form of the ice sheet dynamics and gives the formula for the spatially varying diffusion coefficient as shown in (2.49):

∂H ∂t = ∂ ∂ x (D ∂h ∂ x ) + ∂ ∂ y (D ∂h ∂ y ) + a s -a b with D = A s ρ 2 g 2 H 3 | grad(h) | + 2 5 A ρ 3 g 3 H 5 | grad(h) | 2

Lyapunov Theory

In the introductory chapter of this manuscript, we described control as a toolbox mixing mathematics and engineering, thus offering solutions for a large variety of problems. This is why a key element of our work is in essence mathematical. Lyapunov theory [START_REF] Lyapunov | The general problem of the stability of motion[END_REF] was established by the Russian mathematician Aleksandr Mikhailovich Lyapunov in 1892 [START_REF] Mawhin | Alexandr Mikhailovich Liapunov, The general problem of the stability of motion[END_REF]. His original work focused on the stability of motion. Through this method, one can study the stability of a system by looking at the stability of some measure of the energy of that system. If that measure or norm is stable, meaning it dissipates with time, then the system must be stable as well. This resonates with us because the ice thickness equation (2.49) is highly non-linear and one of its parameters, the basal sliding A s , is unknown and non-measurable, causing the diffusion coefficient D to be unknown as well. The aim of our work is to estimate those parameters by studying the stability of (2.49). We elaborate on our methodology in chapters 3 and 4.

For an in depth description of Lyapunov theory we refer the reader to [START_REF] Vidyasagar | Nonlinear Systems Analysis[END_REF][START_REF] Khalil | Nonlinear Systems[END_REF]. We base our review of the theory on one of the chapters in [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation. 1st[END_REF]. The first problem tackled by Lyapunov consisted in investigating the possibility of choosing the initial values of differential equations small enough so that for all subsequent time instants these equations remain bounded by arbitrarily small limits. Consider a dynamical system:

H = dH dt = F(H, t) H(t 0 ) = H 0 H ∈ R n (2.50)
In order for solutions to exist and be unique, the function F(H, t) needs to be Lipschitz continuous [START_REF] Sohrab | Topology of R and Continuity[END_REF] with respect to H, uniformly in t, and piece-wise continuous in t. In other words, and very crudely, F has to be smooth. We also consider that H ∈ R n is an equilibrium point of (2.50) if F( H, t) = 0 ∀ t. The origin of the system can be shifted towards the equilibrium point, so that H = 0. All the definitions regarding stability (locally or globally), equilibrium, continuity and the behaviour of functions (strictly increasing, decreasing ...) can be found in chapter 4 of [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation. 1st[END_REF]. We do however introduce one definition before stating Lyapunov's basic theorem: Definition 2.4.1. Given a system as in (2.50), the time derivative of a function V(H, t) along the trajectories of the system is defined as:

V H=F(H,t) = ∂V ∂t = ∂V ∂H F (2.51)
The basic theorem of Lyapunov, also known as Lyapunov's direct method:

Theorem 3. Let V(H, t) be a non-negative function with derivative V along the trajectories of the system.

1. If V(H, t) is locally positive definite and V(H, t) ≤ 0 locally in H and ∀ t, then the origin of the system is locally stable.

2. If V(H, t) is locally positive definite, decrescent and V(H, t) ≤ 0 locally in H and ∀ t, then the origin of the system is uniformly locally stable.

3. If V(H, t) is locally positive definite, decrescent and V(H, t) is locally negative definite, then the origin of the system is uniformly locally asymptotically stable.

4. if V(H, t) is positive definite, decrescent and V(H, t) is negative definite, then the origin of the system is globally asymptotically stable.

To be able to quantify the rate of convergence of a system, Theorem 3 can be modified. Below is Lyapunov's exponential stability theorem:

Theorem 4. H = 0 is an exponentially stable equilibrium point of (2.50) if and only if there exist an > 0 and a function V(H, t) which satisfy

α 1 H 2 ≤ V(H, t) ≤ α 2 H 2 V H=F(H,t) ≤ -α 3 H 2 ∂V ∂H (H, t) ≤ α 4 H with {α 1 , α 2 , α 3 , α 4 } > 0 and H ≤ (2.52)
The last theorem describes Lyapunov's indirect method where one can show local stability of the original non-linear system by analyzing the stability of its linearization. Let the non-linear system be as in (2.50). We define the Jacobian matrix of F(H, t) with respect to H evaluated at the origin, as:

A(t) = ∂F(H, t) ∂H H=0 (2.53)
The linearization around the origin of (2.50) is subsequently:

H L = A(t)H L (2.54)
The deviation of the linearized system away from the non-linear one is: it is a locally uniformly asymptotically stable equilibrium of (2.50).

∆(H, t) = H -H L = F(H, t) -A(t)

Relevant Works

Readers of this manuscript can see that our work is situated between glaciology and control theory.

The dynamics that represent the main model come from the ice-thickness equation (2.49), while our approach is based on Lyapunov theory. In general, models can be set-up with known parameters and inputs, run forward in time, thus generating an output that can be analyzed. For example, a given ice sheet model can run under various inputs to see how it might evolve with time. This is a forward problem. However, in cases where parameters are unknown and observational data are available, one can seek to estimate those parameters to fit the observations, this is an inverse problem [START_REF] Hasanoğlu | Introduction to inverse problems for differential equations[END_REF].

For a survey of inversion methods and some of their applications we refer the reader to [START_REF] Yaman | A Survey on Inverse Problems for Applied Sciences[END_REF].

Our work falls within the category of inverse problems as we estimate basal sliding and diffusion to match measurements [START_REF] Fretwell | Bedmap2: improved ice bed, surface and thickness datasets for Antarctica[END_REF] of the surface topography of the Antarctic ice sheet. Even though our methodology is new, the question of such parameter inversion is not. A lot of research has been devoted to tackling such issues in glaciology. In [START_REF] Douglas | The basal stress distribution of Ice Stream E, Antarctica, inferred by control methods[END_REF][START_REF] Macayeal | Basal friction of Ice Stream E, West Antarctica[END_REF] the authors sought to calculate basal friction in order to match measurements of the surface velocity of an Antarctic ice stream. They carried out their inversion by defining a cost function dependant on the difference between simulated and measured surface velocities. Their minimization method is detailed in [START_REF] Douglas | A tutorial on the use of control methods in ice-sheet modeling[END_REF]. The same method was applied in [START_REF] Morlighem | Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model[END_REF] but on a continental scale and results were compared to that of [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF]. The Pollard and DeConto method [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] is a simple and elegant method where the misfit between observed and simulated surface elevation of the Antarctic ice sheet are used to calculate a gain that will iteratively update basal sliding. We offer a quick recap of this method in Chapter 3, and compare our method to it. In [START_REF] Le Clec'h | A rapidly converging initialisation method to simulate the present-day Greenland ice sheet using the GRISLI ice sheet model (version 1.3)[END_REF] an inversion was done to find the distribution of the basal drag coefficient that can reproduce measurements of ice thickness. The authors iteratively corrected this coefficient using measurements and simulated values of sliding velocities. In [START_REF] Konovalov | Inversion for basal friction coefficients with a two-dimensional flow line model using Tikhonov regularization[END_REF], similar to the previously mentioned works, the author inverted for basal friction but constructed a different cost function with a smoothing term based on Thikhonov regularization [START_REF] Nikolaevich | On the solution of ill-posed problems and the method of regularization[END_REF][START_REF] Andrey | Solutions of ill-posed problems[END_REF]. In [START_REF] Mélanie | Estimating basal properties of ice streams from surface measurements: a non-linear Bayesian inverse approach applied to synthetic data[END_REF] both surface velocities and surface topography are used to infer basal characteristics. Here the authors used a different approach. They based their inverse method on non-linear probabilistic Bayesian inference [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF][START_REF] Rodgers | Inverse methods for atmospheric sounding: theory and practice[END_REF]. In this approach, a priori information about the basal properties is expressed as a probability density function and combined with the surface measurements to give a posteriori probability distribution describing the final uncertainty of the estimate. The solution of the inverse problem provides a set of solutions from which the most likely one corresponds to the maximum of the a posteriori probability. In [START_REF] Van Pelt | An iterative inverse method to estimate basal topography and initialize ice flow models[END_REF], interest was shown in inverting for bed topography. Here the authors were inspired by [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] as their simple iterative scheme used misfit between observed and simulated surface topography to update bedrock elevation. But instead of doing their update in a multiplicative fashion as in [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF], they introduce change in an additive fashion while using a weight to scale the update.

Lyapunov methods are powerful tools when stability of systems are of concern. However one of the main challenges facing these methods is the choice of the function V [START_REF] Malisoff | Constructions of strict Lyapunov functions[END_REF]. It is common practice to usually start the search by writing down the energy within a system and to then use one of the methods discussed in the previous section to prove stability. But this choice does not guarantee success, and in our case, as we are dealing with ice sheets, such a concept is not clear cut. Another strategy to find such functions is to define some generic form of V then try to resolve the conditions that will guarantee stability [START_REF] Xu | Exponential Stability and Transfer Functions of Processes Governed by Symmetric Hyperbolic Systems[END_REF][START_REF] Coron | A Strict Lyapunov Function for Boundary Control of Hyperbolic Systems of Conservation Laws[END_REF]. Also, we should keep in mind that the notion of stability is not the same in glaciology and control. Whenever the stability of ice sheets is discussed it usually refers to how ice sheets loose mass. While when we talk about instability, we mean that the given system "blows-up". We decided to steer away from trying to find an elaborate Lyapunov function and instead focused on the misfit between observed and simulated ice thickness and built our V around that.

Our choice is not only supported by its simplicity but by the fact that our ultimate goal is not actually stability but parameter estimation. To our knowledge, basal sliding has not been estimated using such an approach.

As for the diffusion coefficient, it is evident that properly defining it in a diffusive system, describing ice sheets or otherwise, is important as this parameter dictates the dynamics of the system. Identifying a varying parameter is a distributed parameter estimation problem [START_REF] Banks | Estimation techniques for distributed parameter systems[END_REF]. In [START_REF] Kravaris | Identification of Parameters in Distributed Parameter Systems by Regularization[END_REF][START_REF] Yu | Identification of parabolic distributed parameter systems by regularization with differential operators[END_REF], and much like the previously cited works, the diffusion coefficient was estimated by minimizing a cost function while imposing regularization to guarantee both stability and uniqueness of the solution. In the case of ice sheets, one can infer the diffusion coefficient using both the surface elevation and ice thickness [START_REF] Jarosch | Restoring mass conservation to shallow ice flow models over complex terrain[END_REF].

However, this comes with its own set of problems as certain parameters like Glen's flow parameter [START_REF] Cuffey | The physics of glaciers[END_REF] heavily depend on ice temperature. We avoid such an issue and use the method introduced in [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF], which consists of using an observer and an update law that enable the estimation of unknown spatially varying parameters, based on either on-line or steady-state measurements. This method was successfully applied to estimating the diffusion as well as the source term in a Tokamak heat transport model [START_REF] Mechhoud | Adaptive Distributed Parameter and Input Estimation in Plasma Tokamak Heat Transport[END_REF]. We expand this method to the two-dimensional (2D) case and search for the diffusion coefficient using measurements of the surface topography of the Antarctic ice sheet.

Estimating the diffusion coefficient for the Antarctic ice sheet has not been addressed in this fashion before.

Conclusion

We finish this chapter with a quick recap of what was mentioned:

• In Section 2.1 a brief review of physics introduces the momentum balance and mass balance equations. These equations are the basis upon which the dynamics of our system are built.

• In Section 2.2 we develop the balance equations and obtain the ice thickness equation in its generic form.

• In Section 2.3 the shallow ice approximation is introduced. This approximation shapes the dynamics of our system and transforms the ice thickness equation. We find that basal sliding and diffusion are in fact dependant.

• In Section 2.4 the various necessary definitions and theorems pertinent to Lyapunov theory are stated.

• In Section 2.5 a review of relevant topics and advancements related to our work are discussed.

In the next chapters we elaborate on our method and explain how to estimate basal sliding for a one-dimensional model in Chapter 3, and diffusion for both one and two-dimensional models in Chapter 4, as to replicate measurements of the surface topography of the Antarctic ice sheet.

C h a p t e r 3 ESTIMATING THE BASAL SLIDING COEFFICIENT IN A 1D MODEL

Antarctica has been the subject of numerous studies as it encompasses a large variety of phenomena that impact reality around us. The understanding and study of past and future behavior of the Antarctic ice sheet requires a good initialization of the models used to simulate the dynamics of ice.

Initializing such models with the observed surface topography and velocity field requires the knowledge of basal characteristics of the ice sheet. Such characteristics cannot be directly measured.

In this chapter, the focus is on the basal sliding coefficient, here denoted as A s , which describes the sliding of ice sheets over the basal bed. Using a constant A s in ice sheet models forward in time leads to a simulated ice sheet that is not in agreement with the observed one. Alternatively, spatial variations of the basal sliding coefficient can be obtained through an iterative method, thus guaranteeing a simulated ice sheet close to the observation. Our method utilizes the misfit between the simulated and observed ice thickness in order to modify A s . It will be compared to that of Pollard and DeConto [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] as their method also uses such a misfit to iteratively update the basal sliding coefficient. This parameter can be constrained by noticing that it acts as a transport coefficient for the ice thickness H, a quantity whose spatial distribution is measured for Antarctica. As we have seen in the previous chapter, the evolution of H can be modeled by a diffusive nonlinear partial differential equation (PDE) (2.49). This fact is exploited here using a Lyapunov-based technique to find the distribution of A s that ensures the stability and exponential convergence of the modeled H toward its measured value.

ice sheet models can vary in complexity, from models solving the full Stokes equations [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] to simpler models where the Shallow-Ice approximation (SIA) [START_REF] Hutter | Theoretical Glaciology[END_REF], the Shallow-Shelf approximation (SSA) [START_REF] Weis | Theory of shallow ice shelves[END_REF],

or a combination of the two [START_REF] Bueler | Shallow shelf approximation as a "sliding law" in a thermomechanically coupled ice sheet model[END_REF] is used. Such models can also vary in the number of spatial dimensions considered, from spatially distributed depth-dependent 3D models [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] to depth-integrated 2D models [START_REF] Partridge | A moving mesh approach to an ice sheet model[END_REF]. In this chapter we deal with the one dimensional (1D) case, and as we are interested in the ice in contact with the bedrock, we decided to adopt the Shallow-Ice approximation. This approximation does not include sliding, as the ice is considered to be frozen to the bedrock. In order to induce sliding, a sliding law can be added by either choosing a nonlinear Weertman sliding law [START_REF] Weertman | On the Sliding of Glaciers[END_REF] or a linear Coulomb friction law [START_REF] Popov | Contact mechanics and friction[END_REF]; we decided to use the former.

We propose a law based on Lyapunov analysis that could pave the way to extend the study to hybrid and two dimensional ice sheet models. It is worthy to mention that our approach does not imply finding a Lyapunov function that will prove convergence relative to a given control law. In fact, the reverse is done as a simple Lyapunov function is chosen and an appropriate update law is calculated in order to fulfill the required conditions for convergence. As a validation, the performance of our proposed method is compared to the method of Pollard and DeConto [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] in terms of cumulative absolute error, difference between the simulated and observed ice thickness, convergence time and the estimated A s . The robustness to the initial guess of basal sliding coefficient A s0 is also investigated.

The content of this chapter is based on our conference paper [START_REF] Mourad | The Initialization Of Basal Sliding Coefficients For Antarctica A Lyapunov Based Approach[END_REF]. However, our findings are more detailed here, including additional tests and simulations that further the understanding of our method.

It also shows the method's capabilities using real data.

In this chapter we:

Introduce the 1D ice dynamics PDE under SIA Linearize the PDE around a given equilibrium Show the Lyapunov analysis and find the appropriate update law Present the results of the linear and nonlinear models Compare our method to that of Pollard and DeConto

Present the results on real data

Reference Model and Problem Formulation

Recalling and dropping one of the dimensions in equation (2.49), the one-dimensional evolution of ice thickness is described by the diffusive PDE:

H t = Dh x x + a (3.1)
We remind the reader that all partial derivatives in x and t are denoted by subscripts.

In Section 2.3, we saw that SIA enables the use of some simplifying relationships, below we recall them expressed in 1D:

D = -vHh -1 x v = v b + 2A 5 Hτ 3 D v b = A s τ b | τ b | τ b ≈ τ d = -ρgHh x (3.2)
Substituting (3.2) into (3.1) gives:

H t = ρ 2 g 2 A s H 3 h x | h x | x + 2 5 A ρ 3 g 3 H 5 h 3 x x + a ∀ x ∈ [x 1 , x 2 ] and t ∈ [0, T] H(x 1 , t) = H(x 1 ) and H(x 2 , t) = H(x 2 ) H(x, 0) = H(x) + ∆H (3.3)
Our estimation problem is formulated as finding the distributed profile A s (x) such that the solution of (3.3) converges toward H (considered to be an equilibrium). The nonlinear nature of (3.3) clearly does not allow us to directly find an analytical solution, we thus propose to use A s as a feedback gain that stabilizes the linearized dynamics through an iterative method.

Linearized Dynamics

Linearizing the dynamics around an equilibrium point is a common practice when dealing with PDEs.

A first order Taylor series expansion around the equilibrium H = H and A s = Ās is carried out. A similar approach was done for another diffusive system in [START_REF] Mavkov | Distributed Control of Coupled Inhomogeneous Diffusion in Tokamak Plasmas[END_REF]. As previously stated, the aim of this chapter is to find Ās . Expanding (3.3) gives:

H t = f (A s , A s x , H, H x , H x x ) = ρ 2 g 2 A s x H 3 h x | h x | +3A s H x H 2 h x | h x | +2A s H 3 h x x | h x | + 2 5 ρ 3 g 3 A 5H x H 4 h 3 x + 3H 5 h x x h 2 x + a (3.4)
The linearization is carried out as follows:

H t = f + (A s -Ās ) ∂ f ∂ Ās + (A s x -Ās x ) ∂ f ∂ Ās x + (H -H) ∂ f ∂ H +(H x -Hx ) ∂ f ∂ Hx + (H x x -Hxx ) ∂ f ∂ Hxx + a (3.5)
Where f = f ( Ās , Ās x , H, Hx , Hxx ).

Denoting Ãs = A s -Ās , Ãs x = A s x -Ās x , H = H -H, Hx = H x -Hx and Hxx = H x x -Hxx , and noticing that Ht = H t (due to the fact that Ht = 0 as it represents the equilibrium of the system) we obtain the linearized dynamics as:

Ht = c 2 Ãs + c 3 Ãs x + c 4 H + c 5 Hx + c 6 Hxx (3.6)
With the following transport coefficients:

c 1 = ρ 2 g 2 Ās x H3 hx | hx | +3 Ās Hx H2 hx | hx | +2 Ās H3 hxx | hx | + 2A 5 ρ 3 g 3 5 H4 Hx h3 x + 3 H5 hxx h2 x + a c 2 = ρ 2 g 2 3 Hx H2 hx | hx | +2 H3 hxx | hx | c 3 = ρ 2 g 2 H3 hx | hx | c 4 = ρ 2 g 2 3 Ās x H2 hx | hx | +6 Ās Hx H hx | hx | +6 Ās H2 hxx | hx | + 2A 5 ρ 3 g 3 20 Hx H3 h3 x + 15 H4 hxx h2 x c 5 = ρ 2 g 2 2 Ās x H3 | hx | +6 Ās Hx H2 | hx | +2 Ās H3 hxx | hx | hx + 3 Ās H2 hx | hx | + 2A 5 ρ 3 g 3 15 Hx H4 h2 x + 6 H5 hxx hx + 5 H4 h3 x c 6 = 2ρ 2 g 2 Ās H3 | hx | + 6A 5 ρ 3 g 3 H5 h2 x (3.7)
Note that c 1 = f + a = Ht = 0.

The boundary and initial conditions become:

H(x 1 , t) = H(x 2 , t) = 0 H(x, 0) = ∆H (3.8)

Lyapunov Analysis and Update Law

The parameter estimation problem is formulated as a stabilization problem by considering Ãs as the controlled variable that should drive H from ( 3.6) to zero. This is done using Lyapunov analysis on our system. This Section is divided into two parts. The first part deals with the choice of an appropriate Lyapunov function, and the second part focuses on finding A s that guarantees the system convergence.

Candidate Lyapunov Function

The Lyapunov function is often chosen to be an energy-like function [START_REF] Khalil | Nonlinear Systems[END_REF] that needs to be dissipated with time. In the case of the system described by (3.6)-(3.8), we would like the error H between the simulated and reference ice thickness to go to zero, thus a choice is the function:

V = 1 2 ∫ x 2 x 1 H2 dx (3.9)
Theorem 6. The time derivative V t of the function V given by (3.9) verifies: Proof. Differentiating (3.9) with respect to time gives:

V t = ∫ x 2 x 1 H(c 3 Ãs ) x dx + 1 2 ∫ x 2 x 1 H2 c 4 dx - ∫ x 2 x 1 c 6 H2 x dx (
V t = ∫ x 2 x 1 H Ht dx = T 1 + T 2 + T 3 + T 4 (3.11)
where:

T 1 = ∫ x 2 x 1 H c 2 Ãs + c 3 Ãs x dx T 2 = ∫ x 2 x 1 c 4 H2 dx T 3 = ∫ x 2 x 1 c 5 H Hx dx T 4 = ∫ x 2 x 1 c 6 H Hxx dx
Integration by parts is performed on T 3 and T 4 to allow H2 to appear inside the integrals. First, T 3 becomes:

T 3 = ∫ x 2 x 1 c 5 H Hx dx = - 1 2 ∫ x 2 x 1 c 5 x H2 dx (3.12)
Notice that due to the boundary conditions, terms like (c 5 H2 )

x 2

x 1 = 0. The same applies to all subsequent integrations.

Second, T 4 can be expressed as:

T 4 = ∫ x 2 x 1 c 6 H Hxx dx = ∫ x 2 x 1 1 2 c 6 xx H2 -c 6 H2
x dx

(3.13)
Thus, (3.11) can be written as:

V t = ∫ x 2 x 1 H c 2 Ãs + c 3 Ãs x dx - ∫ x 2 x 1 c 6 H2 x dx + ∫ x 2 x 1 H2 c 4 - 1 2 c 5 x + 1 2 c 6 xx dx = ∫ x 2 x 1 H(c 3 Ãs ) x dx + 1 2 ∫ x 2 x 1 H2 c 4 dx - ∫ x 2 x 1 c 6 H2 x dx (3.14)
where the last equality is obtained by noticing that c 4 -1 2 c 5 x + 1 2 c 6 xx = 1 2 c 4 . Also note that c 3 x = c 2 due to the linearization process. This can be obtained by recalling c 6 and differentiating with respect to x:

c 6 x = ρ 2 g 2 2 Ās x H3 | hx | + 6 Ās Hx H2 | hx | + 2 Ās H3 ( h) x x | hx | hx + 2A 5 ρ 3 g 3 15 Hx H4 h2 x + 6 H5 ( h) x x hx (3.15)
We introduce a new variable c 7 by subtracting c 5 from (3.15):

c 7 = c 6 x -c 5 = ρ 2 g 2 2 Ās x H3 | hx | + 6 Ās Hx H2 | hx | + 2 Ās H3 ( h) x x | hx | hx -2 Ās x H3 | hx | -6 Ās Hx H2 | hx | -2 Ās H3 ( h) x x | hx | hx -3 Ās H2 hx | hx | + 2A 5 ρ 3 g 3 15 Hx H4 h2 x + 6 H5 ( h) x x hx -15 Hx H4 h2 x -6 H5 ( h) x x hx -5 H4 h3 x = -3ρ 2 g 2 Ās H2 hx | hx | -2A ρ 3 g 3 H4 h3 x (3.16)
Thus, it was found that c 6 xc 5 + c 7 = 0 ∀ x, which consequently implies that c 6 xxc 5 x + c 7 x = 0. The next step is to find the value of c 7 x :

c 7 x = ρ 2 g 2 3 Ās x H2 hx | hx | + 6 Ās Hx H hx | hx | + 6 Ās H2 ( h) x x | hx | + 2A 5 ρ 3 g 3 20 Hx H3 h3 x + 15 H4 ( h) x x h2 x (3.17)
From the above, we found that c 7 x = c 4 which implies:

c 4 - 1 2 c 5 x + 1 2 c 6 xx = 1 2 c 4
(3.18)

Design of the Update Law

Using the Lyapunov function discussed in Section 3.3, the convergence of (3.6)-(3.8) is ensured by the following theorem.

Theorem 7. Choosing:

Ãs = - 1 c 3 H ∫ x x 1 c 4 (l)dl (3.19)
guarantees the exponential convergence of (3.6)-(3.8) and:

∫ x 2 x 1 H(x, t)dx ≤ ∫ x 2 x 1 H(x, 0)dx e -γt (3.20)
where γ is a positive-definite constant.

Proof. Let Ãs be chosen as in (3.19) and replace it in (3.10). This gives:

V t = - ∫ x 2 x 1 H2 c 4 + H Hx ∫ x x 1 c 4 (l)dl dx + 1 2 ∫ x 2 x 1 H2 c 4 dx - ∫ x 2 x 1 c 6 H2 x dx (3.21)
Using integration by parts, the first integral becomes:

-

∫ x 2 x 1 H2 c 4 + H Hx ∫ x x 1 c 4 (l)dl dx = - 1 2 ∫ x 2 x 1 H2 c 4 dx (3.22)
Replacing (3.22) in (3.21) gives:

V t = - ∫ x 2 x 1 c 6 H2 x dx (3.23)
Notice that c 6 (x) ≥ 0 from (3.7). This enables the use of Wirtinger's inequality [START_REF] Dym | Fourier series and integrals[END_REF] on the above integral:

-

∫ x 2 x 1 c 6 H2 x dx ≤ - c 6 min C ∫ x 2 x 1 H2 dx (3.24)
where c 6 min is the minimum of c 6 (x) and C = (x 2 -x 1 ) 2 π 2 . We now have:

V t ≤ -γ ∫ x 2 x 1 H2 dx with: γ = c 6 min C = c 6 min π 2 (x 2 -x 1 ) 2 > 0 (3.25)
This concludes the proof.

Iterative Calculation of the Sliding Coefficient

The results of Theorem 7 are used to iteratively calculate A s . The main idea is that, starting from an arbitrary guess on Ās , we iteratively update the equilibrium as Ās + Ãs until H becomes as small as possible and Ās becomes as close as possible to the ideal true A s . Dividing c 7 = ∫ x

x 1 c 4 (l)dl by c 3

gives:

1 c 3 ∫ x x 1 c 4 (l)dl = 3 h Ās + 2A ρg h|( h + b) x | (3.26)
Equation (3.19) becomes:

Ãs = -H 3 H Ās + 2A ρg H| hx | (3.27)
Since Ās is our variable of interest and not known a priori, we consider an iterative update law that will push the equilibrium of the linearized dynamics (3.6)-(3.8) toward a solution of the nonlinear dynamics (3.3). Also, A s is kept bounded between Āsmin = 10 -10 m a -1 Pa -2 and Āsmax = 10 -5 m a -1 Pa -2 which respectively represent the cold hard bedrock and the slipperiest deformable sediment. The update is done with the following algorithm:

• Start with an initial guess of basal sliding A s 0 and initialize the nonlinear system with H = H.

• Run the system with the last calculated Ās to get close to an equilibrium (e.g. during a time period sufficiently large with respect to γ) and obtain a new H.

• At iteration i calculate:

Ãs i = -Hi 3 H Ās i + 2A ρg H| hx | (3.28)
and update A s using:

                   Ās i+1 = A s i -Ãs i ∀ x ∈ [x 1 , x 2 ] if Ās i+1 (x) < Āsmin then Ās i+1 (x) = Āsmin if Ās i+1 (x) > Āsmax then Ās i+1 (x) = Āsmax A s i+1 = Ās i+1 (3.29)
• Stop the simulation if

∫ x 2 x 1
| H|dx ≤ , where > 0 is an arbitrarily small scalar constant.

Inverse problems can be ill posed as a unique solution might not exist [START_REF] Engel | Regularization of inverse problems[END_REF]. It is a common practice to add a regularizing term in order to impose some degree of smoothness to the solution [START_REF] Petrat | An inexact Gauss-Newton method for inversion of basal sliding and rheology parameters in a nonlinear Stokes ice sheet model[END_REF]. In our case, all spatial derivatives are calculated using a Gaussian derivative operator, and all transport coefficients (c 2 to c 6 ) are smoothed with a Gaussian filter. Additionally, every new A s is a weighted sum between smoothed and unsmoothed versions of it as shown in:

A s = (1 -k) A s + k A s *
where k ≥ 0 is the weight assigned to the smoothed term A s * (3.30) This filtering helps in rejecting high frequency variations and avoids numerical instabilities when solving for H.

Linear Model Results

The general setup of the experiments consists of two phases. In the first phase a profile is built by solving (3.3). This requires choosing a reference basal sliding coefficient Ās , a bedrock profile b, and a profile of the input a. Then, starting from an initial ice thickness H 0 , the non-linear system is allowed to relax until it reaches steady-state, thus giving us a reference ice thickness H that we treat as our observation in the linear system. Then in phase two, H and (3.28)-(3.30) are used to estimate Ās by updating the system after a fixed number of iterations. This process is shown in Fig. 3.1.

The linear system must be initialized with some H. To do that we start with H = 0 but allow the system to run forward without any updates for some arbitrary number of iterations, this will generate an H 0. Then, at every iteration we calculate s H (the sum of the absolute error H) and its slope. If the absolute value of the slope is less than a chosen > 0 the estimation is stopped and the system is allowed to relax for a fixed number of iterations. We decided to use the slope of s H as our criterion to stop the estimation to be able to assess how fast the system converges. The last relaxation phase is to show that the system does in fact converge to a stable state.

The linearized model (3.6)-(3.8) is discretized and solved using a fully implicit scheme. Four different profiles are tested against three initial conditions A s0 . The update step u s (in steps) and the smoothing weight k are kept constant. The setup of these experiments is detailed in Table 3.1 where X can be one of four profiles A, B, C or D. Note that the values of A s0 are chosen to reflect different behaviors at the start of the simulations. A value of 10 -5 means that the bottom of the ice is very slippery thus sliding is at a maximum, the further this value drops sliding decreases.

The main aim of these tests is to show the capability of the method to retrieve a given profile, and to be able to do so regardless of the choice of A s0 . Also, the updates are done at every step; this will change, for reasons that will be detailed later when we test with the non-linear model in Section 3.6.

The results are presented and discussed below.

Profiles A and B: Setup and Errors

We first start with profiles A and B as they are rather similar. Figures 3.2 and 3.3 contain the plots of all the reference variables needed to setup these profiles: the reference ice thickness H (dashed blue) obtained after running forward the non-linear system, the bedrock elevation b (full brown), the reference surface topography h (full pink), the input mass balance a, and the reference basal sliding coefficient Ās . Profile A represents the simplest of cases: a flat b, a constant input a, and a constant Ās . Profile B is slightly different as it has an inclined bed. Note that in order to better visualize the profiles one needs to keep in mind the dimensions of the represented variables. At first glance it might seem that the produced ice sheets have more thickness than width, but in fact they have a maximum elevation of a few kilometers while being several thousand kilometers wide.

The second pair of figures shows the errors for profile A in Fig. 3.4 and profile B in Fig. 3.5. Each figure contains the errors for each of the chosen initial guesses A s 1 0 = 10 -5 (full green), A s 2 0 = 10 -6 (dashed yellow), A s 3 0 = 10 -8 (dotted purple). Focusing on the evolution of the absolute cumulative error s H , we notice that for the case when the initial sliding is highest (A s 1 0 ), the start of s H is highest as well. Remember that each simulation starts with H = 0, and that the deviation of the system depends on sliding. So, if sliding is high then deformation is fast, while if sliding decreases then the ice will need more iterations to deform. This is why with A s 3 0 the initial error is lower than for the other two values. We also notice that for all three cases s H converges to the same value, thus the obtained misfits H are all equal (right side of Fig. 3.4 and Fig. 3.5). This directly hints that the estimated A s are the same and that our estimation method is robust. However, an additional reason for having such a consistent behavior is the fact that we are using a linear system operating around a given equilibrium. The transport coefficients (3.7) dictate how the system is bound to evolve, prompting the recalculation of these coefficients after every estimation. Once our method produces consistent estimates, the system behaves in a physically consistent manner. In Section 3.3, we proved that the error must go to zero, but in practice some error remains due to numerical integration. Most of the error occurred at the boundaries, but H is quite small compared to H; less than 0.3% of the peak value.

Profiles C and D: Setup and Errors

Profiles C and D introduce some changes to A and B. First, the sliding coefficient A s varies throughout the domain, more so for profile D. The aim of choosing Ās as can be seen in figures 3.6 and 3.7 is to test the capability of our method to estimate more complex distributions of basal sliding. Second, while profile C still has an inclined bed b, undulations were added to that of profile D. And last, the input a in profile D varies as well. All these changes ought to make the dynamics of the system more complex and the job of estimation harder. Much like the two previous profiles, all errors converge towards the same values as the misfits H are all the same (Fig. 3.8 and Fig. 3.9). One striking difference however is that the behavior of s H is different. We notice that the transient behavior of s H is not smooth. The reason for this is quite simple to explain. In order to obtain the given reference, the ice had to slide unevenly across the domain. Transitioning from a constant profile to a varying one is bound to cause sudden changes in H. More so, this sudden transition is visible in both figures as peaks appear in s H and the appearance of those peaks is delayed with slower moving ice. For example, for A s 2 0 the peak happens at around i = 10 3 iterations while for A s 3 0 it occurs closer to Figure 3.2: 1D linear system, profile A: reference variables. On the left, the ice thickness H (dashed blue), surface elevation h (full pink), bedrock b (full brown). On the top right, the surface mass balance a. On the bottom right, the basal sliding coefficient Ās . A s i = 10 4 . In a sense, more error implies more information for the update law to exploit, remember that the update formula (3.28) depends on such an error. So, if the ice begins as slow moving, error/information needs time to build up, thus the update needs more time to enforce changes. This particular issue is why when we will deal with the non-linear model updates are not at every step. In fact we will allow longer relaxation times between updates, so that the non-linear system will have time to adapt to changes in A s . For both profiles, H is higher at the boundaries and is orders of magnitude less than H.

All Profiles: Estimation

Having discussed the setup of each profile and the obtained errors, we can move on to evaluating the performance of the update law in the linear 1D case. Figure 3.10 shows the reference sliding coefficient Ās (full blue), its estimate A s (dashed brown), and L Ãs = log 10 (| Ãs |) (full pink), for all profiles. We decided to plot only one of the estimates for each profile, as regardless of the value of A s0 , the estimates are the same. Of course, this shows that our method is robust to initial guesses. The estimates fit quite well their references, though some error is visible at the boundaries.

Profile C has a more pronounced estimation error on one of its boundaries. However, as we saw in Subsections 3.5.1 and 3.5.2, these errors did not cause significant deviations in H.

We know by now that the error H is needed for the update law to be effective. However, other factors Figure 3.6: 1D linear system, profile C: reference variables. On the left, the ice thickness H (dashed blue), surface elevation h (full pink), bedrock b (full brown). On the top right, the surface mass balance a. On the bottom right, the basal sliding coefficient Ās . can influence (3.28), namely hx and Hx , the slopes of the surface topography and ice thickness respectively. It is clear that in (3.28) if hx and/or Hx approach zero less change will be done on A s . These slopes also heavily affect the transport coefficients (3.7). Since we linearized around an equilibrium (3.5), points with low slopes will change slower relative to others, when the linearized system is run forward in time. All these factors play an important role in the quality of the estimates.

For profiles A and B, we notice such an effect near the center of the domain. Referring back to figures 3.2 and 3.3, we see that the top of the ice sheet occurs at these points. For profile C, h is somewhat flat at two points near X = 2 × 10 3 km and X = 4300 km, this is where L Ãs stands out. As for profile D, it is harder to attribute the pattern of L Ãs to such reasons.

Conclusion of 1D Linear Results

We briefly recap the results and observations of this section. We saw that our method performed well and was able to estimate basal sliding with various degrees of complexity. The obtained error on ice thickness was found to be low and concentrated at the boundaries of the domain. We tested against different values of initial guesses of basal sliding and our method proved to be robust in this regard as all estimates were equal. However, the method and profile depend on the produced misfit H and the slopes of surface topography and ice thickness. This is why at points where H and/or these slopes are close to zero the method's effectiveness slightly decreases.

Non-linear Model Results

The tests for the 1D non-linear model (3.3) are similar to those of Section 3.5. In addition to using the same profiles, we compare our method (Method 1 or M1) to that of Pollard and DeConto [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] (Method 2 or M2), and focus on the following aspects:

• robustness to A s 0 ;

• effects of changing the update step u s ;

• effects of changing the smoothing weight k;

• the performances of both methods. Methods 1 and 2 are similar in that they both utilize the misfit in ice thickness to enforce changes in the estimates of basal sliding. However, unlike Method 1 where the change in A s is done in an additive fashion, Method 2 seeks to iteratively apply a varying gain on A s . At each update, every 3.10: 1D linear system, estimations of basal sliding for all profiles: log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference (full pink). point of the new A s is calculated as the product of its current value and a gain dependant on H. Thus at iteration i one can write:

A i+1 s = A i s 10 ∆z with ∆z = max[ -1.5 , min[ 1.5 , H/H inv ] ]
where H inv is a scaling factor (3.31) We see that in (3.31) 10 ∆z is limited by two values (0.03 and 30), the authors enforced such a limit to avoid overshoots and ensure convergence of their method. Also, the scaling factor H inv reflects the weight with which the misfits influence the update. In other words, it affects the speed of convergence, though a lower value of H inv does not necessarily imply faster convergence. In fact, increasingly lowering the scaling factor might cause numerical instabilities in the simulations.

No. Profile Method

= 10 -6 , A s 2 0 = 10 -8 , u 1 s = 10 2 , u 2 s = 5 × 10 2 , k 1 = 10 -3 , k 2 = 5 × 10 -3 , k 3 = 10 -2 , H 1 inv = 4 × 10 3 , H 2 inv = 2 × 10 3 , H 3 inv = 8 × 10 3
. The non-linear system was solved using the MATLAB solver ode23s [START_REF] Shampine | The MATLAB ODE Suite[END_REF].

Robustness to the Initial Guess

It was established in Section 3.5 that Method 1 is indeed robust to initial values of A s 0 . We now show that it is still the case for the non-linear system. In this section we compare the results of experiments 1-4. Between each experiment we have changed the value of A s0 , see Table 3.2. We start by looking at Fig. 3.11 where we show s H for all profiles. We quickly notice that for both methods the errors do not converge to the same value, as was the case for the linear system, though it is clear that the errors have a decreasing trend. However, this does not imply that the methods are not robust. In Fig. 3.12 and Fig. 3.13 we show the reference sliding coefficient Ās (full blue), its estimate A s (dashed brown), and log 10 of their absolute difference when A s 1 0 (full pink) and A s 2 0 (dotted blue), for methods 1 and 2 respectively. The estimation errors are low and quite close, this is why we plot only one estimate for each reference.

Both methods performed well and prove to be robust to initial guesses. Contrary to the linear case, the boundaries do not carry errors. In fact, for the non-linear case, the estimation errors generally dropped for all profiles, except for profile C (M1) where some stiffness emerged around X = 5300 km. The reason for its appearance is numerical as for later tests it is no longer present. Moving to Fig. 3.14, all the misfits for experiments 1-4 are plotted. Though all misfits are very low, M1 is able to produce less error for profiles A, B and D. For profile C, M2 performed better. Also, the effects of the slopes hx and Hx (discussed in the Section 3.5) are no longer present as the non-linear system does not depend on transport coefficients (3.7) to run forward in time.

Sensitivity to the Update Step

We now examine the effect of changing the update step u s . In experiments 2 and 4, updates were done after u 1 s = 10 2 iterations, in experiments 5 and 6 we increase it to u 2 s = 5×10 2 iterations. Doing so will give the system more iterations to evolve, thus accumulating more H before the next update.

Though this aides the estimation process and increases the numerical stability of the system, it comes at the expense of prolonging the convergence time. This can be visible in Fig. 3 (full dark blue), and 6 (dashed purple) to 4 (dashed red). However, the more interesting comparison is between methods. In the previous set of tests, M2 showed some advantage over M1 for profile C. This is no longer the case once the update step was increased. Updating A s using M1 appears to be more aggressive than M2. Allowing the system to relax longer when using M1 will benefit the estimation process.

Both methods produced good estimates, see Fig. 3.17 for M1 and 3.18 for M2. The quality of those estimates did not drastically change. Though for M1, the stiffness that appeared in the estimation of profile C decreased. Hinting once again that increasing u s can benefit the estimation.

Sensitivity to the Smoothing Weight

We now focus on the effect of changing the smoothing weight k. To do so we use the results from experiments 2 (M1, k 1 = 10 -3 ), 4 (M2, k 1 = 10 -3 ), 7 (M1, k 2 = 5 × 10 -3 ), 8 (M1, k 1 = 10 -2 ), 9 (M2, k 2 = 5 × 10 -3 ), and 10 (M2, k 3 = 10 -2 ).

We start with the evolution of the cumulative error, seen in Fig. 3.19. We have plotted s H for all six experiments and for both methods. We find that profiles A and B benefited from more smoothing. This is due to the fact that Ās for these two profiles is constant across the domain, so with more smoothing, the estimate A s became more uniform and approached the reference. This of course is no longer the case for profiles C and D as they have spatially varying sliding coefficients. Logically, Figure 3.17: 1D non-linear system, all profiles, experiments 5-2: estimations of basal sliding. log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference for u 1 s (full pink) and u 2 s (dotted blue). L Ãs is shown in dotted orange. Another confirmation that profiles A and B benefited from an increase in k is the the fact that L Ãs decreases for experiments 7-8 and 9-10. Profiles C and D show the opposite except at the region that showed stiffness for M1 (X = 5300 km) smoothing did indeed remove the high variations in A s .

Sensitivity of M2 to the Scaling Factor

The goal of this subsection is show that M2 is sensitive to the choice of scaling factor H inv . This sensitivity is however non-intuitive, as an increase in H inv does not necessarily imply faster convergence, lower errors hence better estimates. In Fig. log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference for k 2 (full pink), k 3 (dotted blue), and k 1 (dotted orange). log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference for k 2 (full pink), k 3 (dotted blue), and k 1 (dotted orange). But as clearly depicted in Fig. 3.23 H 2 inv produced the least errors for two of the four profiles. It was shown in [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] that the quality of the results depends on the combination of both H inv and u s .

Conclusion of 1D Non-linear Results

The non-linear 1D case was tested much like the linear one. In addition to confirming the robustness of our method to initial guesses A s 0 , we investigated the sensitivity to changes in update step u s , and to changes in the smoothing weight k. We also compare our method (M1) to an already established one, that of Pollard and DeConto [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] (M2).

We stated at the start of Section 3.6 that the reason we chose to compare our method with M2 is because they both utilise the misfit H to update A s . All along analyzing the non-linear results, M2 was present to mainly prove one thing, that if both methods behave in a similar fashion, then M1 is a new valid approach.

Our method can introduce stiff changes into the system, and one way to deal with this is to allow larger relaxation times between the updates by increasing the update step u s . The obvious drawback here is the need for longer runs. However, larger update steps have another benefit. As our update formula (3.28) depends on the misfit H, allowing the system to further run forward produces more accurate misfits, thus more accurate estimations A s . Another way of dealing with numerical stiffness is filtering. We thus defined our new estimates as a weighted sum between smooth and non-smooth terms (3.30). A trade-off arises between updating more frequently and smoothing more, and allowing the system to relax longer after each update (pushing for longer simulation times) and smoothing less. We also briefly tested changing the scaling factor H inv in M2. We noticed that M2's results can be improved by choosing a better combination of H inv and u s .

For all tests, M1 produced results close to or better than M2. This clearly implies that our method is a new viable option for estimating basal sliding. In the next section we prove this point further as we test the method on real Antarctic data.

Experimental Results

With the linear and nonlinear tests done in sections 3.5 and 3.6, we now carry tests on real data.

Early in our work, we selected cross-sections of the Antarctic ice sheet and tried to retrieve the selected profiles of surface topography. For example, a 10 km data resolution holds the data in matrices 561 × 561 large. To get a cross-section we could select the 280 th row of the matrix and simulations would be done to estimate the basal sliding distribution that could reproduce this profile.

This of course proved to be inaccurate because the Antarctic ice sheet does not flow along such dimensions. This is why we turned our attention to flowlines as they are one-dimensional and within them the ice flows following one direction. Flowlines can be useful to describe certain aspects of glaciers [START_REF] Oerlemans | A flowline model for Nigardsbreen, Norway: projection of future glacier length based on dynamic calibration with the historic record[END_REF][START_REF] Riffenburgh | Encyclopedia of the Antarctic[END_REF], and can be sometimes coupled to other geophysical models [START_REF] Gladstone | Calibrated prediction of Pine Island Glacier retreat during the 21st and 22nd centuries with a coupled flowline model[END_REF] or used to study phenomena like flowstripes [START_REF] Glasser | Longitudinal surface structures (flowstripes) on Antarctic glaciers[END_REF]. Our aim is not to model them but to utilize the concept of a flowline in order to produce estimates of basal sliding along one dimension. We use the non-linear model (3.3) and the data from the Bedmap2 data set [START_REF] Fretwell | Bedmap2: improved ice bed, surface and thickness datasets for Antarctica[END_REF].

To extract flowlines we select basins within the Antarctic map. These drainage systems are as defined in [START_REF] Zwally | Antarctic and Greenland drainage systems[END_REF] and are shown in Fig. 3.24. We test our method on flowlines from basins 3, 10, 17, 21 and 22, noting that we treated basins 21 and 22 as one in order to cover a larger area of the map.

For each selected basin, and given a starting point on the map, we use measurements of the vertically averaged surface velocities v and its spatial derivatives to extract the coordinates of the flowline using the stream2 function in MATLAB. Then the corresponding surface elevation h, ice thickness H, bed elevation b, and surface mass balance a are extracted from the map. The data had to be interpolated in order to have evenly spaced data points. The data needs to be processed one more time in order to remove repetitive points and those not belonging to the grounded ice sheet. Thus a total of 10 flowlines are selected for each of the chosen basins. The general setup of the simulations is similar to the previous sections. The update step u s is fixed at 5 × 10 2 iterations to give the system enough time to react to changes in A s , and the value of the smoothing weight k is slightly increased to 5 × 10 -2 . The main difference however is that the data is of course not generated from the non-linear model (3.3). We focus on minimizing the generated misfit H and on the estimates of basal sliding. For each basin, we show the selected flowlines and their mean absolute error Hm , as well as a selection of individual lines in order to better understand the outcome of our method: Fig. 3.25-Fig. 3.27 for basin 3, Fig. 3.28-Fig. 3.30 for basin 10, Fig. 3.31-Fig. 3.33 for basin 17, and Fig. 3.34-Fig. 3.36 for basin 21/22. We set a threshold of 20 m for Hm to qualify the success of our method.

For basins 3, 17, and 21/22 the results are good as most of the obtained Hm are below 20 m. This is not the case for basin 10 for reasons that we are to discuss shortly. But first we look at the results of the estimation of basal sliding for the prior three basins. We notice that in the cases of successful estimation, like in Fig. 3.26, the basal sliding profile always tends to increase from the left to the right of the domain, keeping in mind that the leftmost point of the domain is the starting point of a flowline.

This feature is quite important as it indicates that points closer to the interior of the ice sheet are slower than those closer to the boundary of the grounded ice sheet [START_REF] Rignot | Ice Flow of the Antarctic Ice Sheet[END_REF]. This is quite known and the fact that our method reproduces this phenomenon is a good indicator. Another interesting observation is in the case of basin 21/22. Here, not only are all the estimations successful, but they are consistent as well in the sense that similar topographical characteristics will give results that are close to each other. Looking at Fig. 3.35 and Fig. 3.36 we can see that the bedrock profiles are rather similar and that the paths of flowlines 1 and 2 are indeed in close proximity, Fig. 3.34. Because of this, the two estimated basal sliding profiles have similar distributions. Basin 17 illustrates the limitations of our method. Though most of the flowlines in this basin have a low Hm , we see that the parts of the domain with the larger misfits H tend to correspond to those where the estimated basal sliding has saturated at one of the two limits Āsmin = 10 -10 m a -1 Pa -2 or Āsmax = 10 -5 m a -1 Pa -2 . This is visible in Fig. 3.33 as H is high at two distinct regions where A s saturated at Āsmin . Remember that Āsmin implies that there is no or very little sliding occurring, in other words that the base of the ice sheet is very cold. A very important point to make here is that this saturation always produces a negative misfit, one where the estimated ice thickness is less than its reference. This deficit in ice thickness in places where the ice base is supposed to be colder, pushing the ice sheet to move more slowly thus allowing more accumulation at the top, strongly hints that lateral flows feed these regions. However, in the model we used we do not account for such flows, hence the deficit that causes our method to produce such saturated estimates. It is for this exact reason that the estimates in basin 10 were unsuccessful. We do not show the plots for all the flowlines in basin 10, but in every case the misfit H is negative and A s is stuck at Āsmin . We show the results of flowlines 1 and 4 in Fig. 3.29 and Fig. 3.30. For flowline 1, the estimated A s saturated for the left half of the domain, and looking at the corresponding misfit, we see that it is indeed negative. The results of flowline 4 are more extreme as A s = Āsmin for the entirety of the domain, and once again H is correspondingly negative.

Conclusion of 1D Experimental Results

Having studied test cases in previous sections, we moved on to test our method on real topographical measurements of the Antarctic ice sheet ( [START_REF] Fretwell | Bedmap2: improved ice bed, surface and thickness datasets for Antarctica[END_REF]). One-dimensional flowline data were extracted using the measured vertically averaged surface velocities. From a selection of four different basins Fig. 3.24

we extracted ten different flowlines from each one of them and then proceeded to test our estimation method. The update step u s was fixed at 5 × 10 2 iterations, the smoothing weight k was constant at 5 × 10 -2 . The method produced satisfying results as it highlighted both positive and negative behaviours that were not obtained during the study cases (linear or non-linear).

First, we saw that for three of the four selected basins (basin 3, 17, and 21/22) the results were mostly positive. The misfit in ice thickness was low and the reproduced topographies fit quite well the measurements. An interesting feature is the consistency of the obtained basal sliding profiles for flowlines that are physically close to each other, like in flowlines 1 and 2 of basin 21/22 Fig. 3.35-Fig. 3.36. Another feature is the tendency of the estimated basal sliding profiles to increase from left to right of the physical domain. This happens because the staring point, which is the leftmost point in each graph, is the one closer to the interior of the ice sheet where the ice moves slower. The more we approach the edge of the domain the faster ice moves giving rise to higher values of A s .

Second, we saw one limitation of our method. Whenever the obtained misfit in ice thickness was high and negative (H < H) and the corresponding values of A s were always equal to its minimum Āsmin . In fact, the method tried to compensate the missing inflows of ice into the domain by forcing A s to be very low, which means that there is no sliding. Such flows were not introduced in the test model used, thus the simulations failed in producing viable estimates of basal sliding.

Conclusion

In this chapter we have introduced a new inverse technique to retrieve the distribution of basal sliding coefficients for 1D grounded ice sheets.

• In Section 3.1 the one-dimensional evolution of ice thickness equation (3.1) is set-up under SIA (3.2) thus forming our non-linear system (3.3).

• In Section 3.2 the nonlinear system is linearized around an equilibrium. The linear system is described by the equations (3.6)-(3.8).

• In Section 3.3 we carry a Lyapunov analysis and a new update law (3.27) is formulated such that convergence of the linear system towards a given steady-state is guaranteed as is shown in Theorems 6 and 7.

• In Section 3.4 we expanded on the way A s is to be iteratively calculated.

• In Section 3.5 we tested our method on the linearized system (3.6)- (3.8). This was done by choosing a profile for A s and running forward in time the nonlinear system (3.3) until it reached a steady-state. This generated our fictitious measurements. The procedure is explained through the block diagram in Fig. 3.1. Four different profiles were chosen for A s , and our method was tested against different values of initial guesses of basal sliding A s 0 . For all profiles and for all the different initial guesses, the estimates were exactly the same and fit well the references of basal sliding.

• In Section 3.6 we carried out more extensive tests on the nonlinear system (3.3). Alongside the initial guesses, the update step u s , and the smoothing weight k were changed. We also compared our method to that of Pollard and DeConto [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] as it also updates A s based on the misfit H (3.31). This set of tests is summarized in Table 3.2. Through these tests we saw that both methods give similar results. Both are robust to initial guesses, though our method is more aggressive in applying changes to A s . This caused our method to produce slightly better results but at the expense of some stiffness in the estimates. This can be countered by allowing larger relaxation periods between the updates and/or smoothing further. In other words, our method can benefit from larger u s and k.

We summarize the results of this section in Table 3.3 where the comparisons between the experiments are detailed, and in Table 3.4 to Table 3.7 we show the results for all experiments and all profiles, more specifically we show the steady state cumulative error s H , the average percentage absolute error in ice thickness, and the average percentage absolute error in A s .

• In Section 3.7 tests are done with real measurements of the Antarctic ice sheet [START_REF] Fretwell | Bedmap2: improved ice bed, surface and thickness datasets for Antarctica[END_REF]. In order to reduce the 2D data into 1D we had to turn our attention to flowlines [START_REF] Oerlemans | A flowline model for Nigardsbreen, Norway: projection of future glacier length based on dynamic calibration with the historic record[END_REF][START_REF] Riffenburgh | Encyclopedia of the Antarctic[END_REF]. However, we had to first find such structures in the available data. We selected 10 flowlines for each of four different basins. Representative samples to better evaluate the effectiveness of our method were discussed. We saw that even though a large portion of the results generated low errors, the main disadvantage was the saturation of A s at one of its two predefined limits Āsmin or Āsmax . The main reason this happened is due to the inability of the used model to replicate some of the physics necessary to simulate lateral flows into flowlines. This pushed our method to produce A s at its limits. We propose that further simulations should be done with more dedicated models.

In the next chapter we turn our attention to estimating the diffusion coefficient D in one and two dimensions.

Experiments the steady-state error s H , the average percentage absolute errors in ice thickness and estimated sliding coefficient.

As previously mentioned, inverse problems emerge in a multitude of research topics as they consist of estimating variables of interest using available data. This is particularly true for fields involving natural phenomena as some measurements might not be available or impossible to gather. For example, the friction or sliding coefficients of ice sheets over the underlying bed [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF], the gravity of planetary bodies [START_REF] Tricarico | Global Gravity Inversion of Bodies with Arbitrary Shape[END_REF], or aquifer parameters [START_REF] Yeh | Review of Parameter Identification Procedures in Groundwater Hydrology: The Inverse Problem[END_REF]. Identifying a diffusion coefficient is a classical example of parameter inversion [START_REF] Sakthivel | Inverse problem for the reaction diffusion system by optimization method[END_REF][START_REF] Richter | An Inverse Problem for the Steady State Diffusion Equation[END_REF][START_REF] Isakov | Some inverse problems for the diffusion equation[END_REF]. In chapter 2, we mentioned some methods developed to solve such a problem. We base our approach on the method introduced in [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF], which consists of using an observer and an update law that enable the estimation of unknown spatially varying parameters, based on either on-line or steady-state measurements.

In the previous chapter 3, we found the sliding coefficient using a Lyapunov-based technique for a 1D model. Scaling that technique for the 2D case proved to be particularly hard, this is why we shift our focus to the diffusion coefficient D of the mass continuity equation 4.1. Calculating D can lead back to finding A s .

∂H ∂t = ∂ ∂ x (D ∂h ∂ x ) + ∂ ∂ y (D ∂h ∂ y ) + a (4.1)
with:

D = A s ρ 2 g 2 H 3 | grad(h) | + 2 5 A ρ 3 g 3 H 5 | grad(h) | 2 (4.2)
When dealing with partial differential equations an important aspect is the stability of the system they describe. One approach is based on Lyapunov stability. We used this method in the previous chapter. Such techniques are widely used, for example in [START_REF] Gahlawat | Bootstrap current optimization in tokamaks using sum-of-squares polynomials[END_REF][START_REF] Meyer | Stability analysis of parabolic linear PDEs with two spatial using Lyapunov method and SOS[END_REF] stability of 2D parabolic PDEs with polynomial varying parameters was shown and estimates of the decay rate were given, and in [START_REF] Shirinabadi | Lyapunov stability analysis of special class of PDE systems[END_REF] a Lyapunov stability analysis was conducted on hyperbolic and parabolic PDEs. In our case, stability is inherently addressed through the use of the method from [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] as a Lyapunov function is shown to have a negative time derivative for all time. Often, Linear Operator Inequalities (LOIs) or Linear Matrix Inequalities (LMIs) can be solved numerically in order to verify Lyapunov stability.

Another important issue when dealing with inverse problems and numerically calculating varying parameters is that such problems are ill-posed [START_REF] Kravaris | Identification of Parameters in Distributed Parameter Systems by Regularization[END_REF], and that major oscillations in the estimated coefficients must be avoided [START_REF] Lions | Some aspects of modelling problems in distributed parameter systems[END_REF]. This is why in [START_REF] Kravaris | Identification of Parameters in Distributed Parameter Systems by Regularization[END_REF] a smoothing functional was added to the cost function. In our case, as we do not explicitly minimize a cost function, a simple way to induce some regularity into the estimates is to have a weighted sum between regularized and unregularized versions of the coefficient. For the 1D case, we used Tikonov regularization [START_REF] Andrey | Solutions of ill-posed problems[END_REF] when calculating the diffusion coefficient; the details of this approach are presented in this chapter. As for the 2D cases, we decided to use a low pass filter on the estimates.

Lyapunov Function and Update Law

We recall the 2D ice-thickness equation (2.49), set Dirichlet boundary conditions, and choose some h 0 as the initial condition :

∂H ∂t = ∂h ∂t = ∂ ∂ x (D ∂h ∂ x ) + ∂ ∂ y (D ∂h ∂ y ) + a h(x, y, t) = h b ∀ x,y ∈ Ω b and t ≥ 0 h(x, y, 0) = h 0 ∀ x,y ∈ Ω (4.3)
Here x and y are the two spacial coordinates, Ω and Ω b are the spacial domain and its boundary, respectively, t is time, H(x, y, t) is the ice thickness, h(x, y, t) is the surface topography, h b is the value of h at the boundary, D(x, y) > 0 is an unknown diffusion coefficient, and a(x, y) is a known source term. Just like previous chapters, in what follows all partial derivatives are denoted by subscripts; e.g. ∂h ∂ x = h x . As mentioned before, D is to be estimated in order to replicate as closely as possible measurements of h. The method employed was introduce by Orlov [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF], and consists of introducing an observer with a correction term with gain v 0 , alongside an update law with gain v 1 that uses the gradients of the difference between the observer and the system. Also, we introduce either Tikhonov regularization or a low pass filter into the estimates of D in order to guarantee a degree of smoothness, thus keeping the numerical experiments stable.

Before discussing the adaptive identification process, it is important to note that all derivations, theorems and proofs can be scaled down to the 1D case without any loss of generality.

The main purpose of an observer is to mimic a system and correct for differences caused by unknown states [START_REF] Ellis | Chapter 1 -Control Systems and the Role of Observers[END_REF][START_REF] Meurer | Control and observer design for nonlinear finite and infinite dimensional systems[END_REF]. In our case, the state, which is the ice thickness, is known but we do not know D. This is why for the system in (4.3) we introduce the following pair:

The observer:

ĥt = D ĥx x + D ĥy y + a -v 0 (t) ĥ -h ĥ(x, y, t) = h b ∀ x,y ∈ Ω b and t ≥ 0 ĥ(x, y, 0) = ĥ0 ∀ x,y ∈ Ω (4.4)
and update law:

Dt = 1 v 1 (t) ĥ -h x ĥx + 1 v 1 (t) ĥ -h y ĥy D(x, y, 0) = D0 (x, y) D0 > 0 ∈ C 1 (Ω) (4.5)
with v 0 (t) ≥ 0 ∀ t and 1 v 1 (t) > 0 ∀ t as adaptation gains. Note that (4.4) and (4.5) can be used in both cases when on-line or steady-state measurements are available. The only difference between both cases is that during steady-state ∂h ∂t in (4.3) is set equal to zero, but this does not impact our procedure as the right-hand side of (4.3) is still used in the proof (as will be seen next). Also, note that our cases differ from those in [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] in that only one spatially varying parameter is being estimated as opposed to two or three, and that v 0 and v 1 are functions of time. The authors of [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] do mention that their method is scalable to more than one dimension. However, it is useful and necessary to show the derivations in the context of our problem.

Before showing that (4.4) and (4.5) enable the convergence of D towards D and ĥ towards h, some deviation variables and their time derivatives must be introduced: h = ĥh

Dt = Dt -D t = Dt = 1 v 1 (t) hx ĥx + 1 v 1 (t)
hy ĥy

V t = ∬ Ω -D h2 x -D h2 y -v 0 h2 + v 1 t D2 dx dy ≤ 0 (4.11)
∀ t ≥ 0 along the solutions of (4.3)- (4.5).

In what follows, we assume that the terms within the double integrals are bounded and measurable.

Proof. Differentiating (4.10) with respect to time gives:

V t = T 1 + T 2 + T 3
where:

T 1 = ∬ Ω ht h dx dy T 2 = ∬ Ω 1 2 v 1 t D2 dx dy T 3 = ∬ Ω v 1 Dt D dx dy (4.12)
Integration by parts is performed on T 1 , yielding:

T 1 = ∬ Ω D ĥx x h + hx D x h + D ĥy y h + hy D y h -v 0 h2 dx dy
Dealing with each term inside the double integral individually and denoting by Ω x , Ω y ∈ Ω b the x and y boundaries on the domain respectively, we get:

               ∫ Ω x hx D x h dx = - ∫ Ω x D h2 x dx ∫ Ω x D ĥx x h dx = - ∫ Ω x
hx ĥx Ddx since ĥb = 0 (from the boundaries of (4.4)).

The same applies to the y direction; while keeping in mind that the double integral can be flipped because the terms inside T 1 's integral are bounded and measurable. T 2 remains as in (4.12), while T 3 needs to be expanded as:

T 3 = ∬ Ω hx ĥx D + hy ĥy D dx dy
Replacing the above in (4.12) gives:

V t = - ∬ Ω D h2 x + D h2 y + v 0 h2 - 1 2 v 1 t D2 dx dy (4.13)
In order to guarantee that (4.13) is always ≤ 0, all we need to do is to choose v 1 (t) to be a decreasing positive function.

The condition on v 1 can be satisfied by having v 1 of the form:

v 1 = α e -βt ; α, β > 0 (4.14)
This concludes the proof.

Regularization

A major issue in inverse problems is high variations in the estimates, which can lead to numerical instability and leave artifacts in the estimated variables. In our case, and as the update law (4.5) contains spatial derivatives of the difference between h and ĥ, oscillations in D might appear due to the amplification of the error. Many methods exist to overcome this issue [START_REF] Knowles | Methods for numerical differentiation of noisy data[END_REF]. A common one would be to use Tikhonov regularization [START_REF] Andrey | Solutions of ill-posed problems[END_REF], or a smoothing filter for the estimates as they get generated. In [START_REF] Pattyn | Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0)[END_REF] a Savitsky-Golay filter [START_REF] Savitzky | Smoothing and Differentiation of Data by Simplified Least Squares Procedures[END_REF] was used to remove high-frequency variations in the basal sliding coefficients estimates. In [START_REF] Eilers | A Perfect Smoother[END_REF] penalized least-squares are used to smooth 1D data, and in [START_REF] Eilers | Fast and compact smoothing on large multidimensional grids[END_REF] both 1D and 2D data are handled with P-splines. Even though [START_REF] Eilers | Fast and compact smoothing on large multidimensional grids[END_REF] can handle data stored in large matrices (500 × 500), it still needs some computational time before producing a smooth version of the original matrix. This is not feasible in our case as the smoothing needs to be done frequently. In [START_REF] Chartrand | Numerical differentiation of noisy, nonsmooth data[END_REF] total-variation regularization is used; even though this method proves to be effective in denoising data and could be used in our case to directly regularize the spatial derivatives in (4.5), the method requires a specified large number of iterations before reaching convergence.

An additional reason for seeking to regularize the estimates is that (4.5) depends on the spatial derivatives of h, and if h and ĥ happen to have similar slopes around certain points then (4.5) will fail to generate proper values at those points as identifiability is lost due to the underconstrained nature of the problem. This particular issue was discussed in chapter 3. This is why, for the 1D case, we consider such points as outliers, and smoothing is done by calculating regularized derivatives of D, as defined in [START_REF] Knowles | Methods for numerical differentiation of noisy data[END_REF], then integrating back to obtain a smoothed version D * . As an example, assume that the regularized derivative y x of a function y needs to be calculated on a domain [x 1 , x 2 ]. This can be achieved by minimizing the functional:

E(y x ) = A 1 y x -ŷ 2 + α A 2 y x (4.15)
Where ŷ(x) = y(x)y(x 1 ), A 1 is a discretized integration matrix, α is a regularization parameter, A 2 is a differential operator, and • represents the Euclidean norm. Calculating the spatial derivative can thus be reduced to solving the system of linear equations below. Keeping in mind that for a given α, the term multiplying y x must remain invertible.

A T 1 A 1 + αA T 2 A 2 y x = A T 1 y (4.16)
Using this method, y x can be integrated over [x 1 , x 2 ] to obtain a smooth version of y. For the 2D case, we opted for an averaging filter with a fixed window size. In both cases (as in Section 3.30) a weighted sum of D and D * is used to obtain the regularized estimate:

D = (1 -k) D + k D * (4.17)
where k ≥ 0 is the weight assigned to the smoothed term D * .

Results of the 1D Study Cases

The numerical tests for the 1D study cases are done to first illustrate the method, and second to understand the the effects of changing some of the simulation variables. We focus on the adaptation gain v 1 , the smoothing weight k, the regularization parameter α, and the correction gain v 0 . We kept v 0 = 0 for the first eight experiments as we would like to understand how the update law works and how sensitive it is to the aforementioned parameters. From those eight experiments, we then select the one with the best result and study the effect of changing v 0 on it. Updates are done every 10 steps, the rest of the variables are set as shown in Table 4.1. A summary of the obtained mean absolute error is shown in Table 4.2.

We use the same test profiles described in Section 3.5. One major change is that the system (4.3)

does not explicitly need A s to run forward, but rather uses the diffusion coefficient D. This is why the four profiles obtain their distributions of D by using (4.2). This equation requires distributions of surface elevation h and basal sliding Ās . We get this pair directly from the simulations done in the previous chapter. We thus obtain a distribution D for diffusion. However, we smooth D because of the influence of the slopes of surface elevation on the values of diffusion. As can be seen from This is due to having a high value of v 1 with a low value of k. In other words, higher reliance on updates with less smoothing. We explain further the sensitivities of the results to the simulation variables in the following subsections.

Sensitivity to the Adaptation Gain

In this subsection we focus on the sensitivity of the results to changes in the adaptation gain v 1 , and test two values as noted in Table 4.1. In Fig. 4.3 we show the evolution of s h through the simulations.

The full lines represent experiments 1 to 4, while the dashed lines represent experiments 5 to 8. We notice that s h is always lower for higher values of v 1 . The experiments can be grouped in the following pairs for comparison: (1,5), [START_REF] Huybrechts | The Antarctic ice sheet and environmental change: a three-dimensional modelling study= Der antarktische Eisschild und globale Umweltveränderungen: eine dreidimensionale Modellstudie[END_REF][START_REF] Weertman | On the Sliding of Glaciers[END_REF], [START_REF] Evans | Partial differential equations[END_REF][START_REF] Ribe | Coiling of viscous jets[END_REF], and (4,8). To better illustrate the effect of changing v 1 we select the pair [START_REF] Hutter | Theoretical Glaciology[END_REF][START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF], but note that the same observations apply to the other three. We compare in Fig. 4.4 the distributions of | h |, and in Fig. 4.5 the estimates of diffusion. We can see that the estimates of experiment 5 fit better the chosen distribution of diffusion, thus causing less misfit. This effect is more noticeable for profiles C and D. We also notice that for the two profiles (A and B) with less variations in D, increasing v 1 had less impact on the estimates. This hints that regularization must play a part in the quality of the estimation. We look into this effect in the next subsection.

Sensitivity to the Smoothing Weight

In this subsection we focus on the sensitivity of the results to changes to the smoothing weight k.

This weight is varied as shown in Table 4.1. The experiments are grouped into the following pairs for comparison: (1,3), [START_REF] Huybrechts | The Antarctic ice sheet and environmental change: a three-dimensional modelling study= Der antarktische Eisschild und globale Umweltveränderungen: eine dreidimensionale Modellstudie[END_REF][START_REF] Greve | Comparison of Numerical Schemes for the Solution of the Ice-Thickness Equation in a Dynamic/Thermodynamic Ice-Sheet Model[END_REF], [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF][START_REF] Ribe | Coiling of viscous jets[END_REF], and [START_REF] Weertman | On the Sliding of Glaciers[END_REF][START_REF] Ribe | Stability of liquid rope coiling[END_REF]. As in the previous subsection 4. [START_REF] Evans | Partial differential equations[END_REF].1, we start with Fig. 4.6 by looking at the evolution of s h for the four experimental pairs. Once again, the results show a clear distinction between the experiments with k = 10 -4 and k = 10 -3 , with the lower value producing less errors. In Fig. 4.6, the dashed lines correspond to k = 10 -3 , while the full line to k = 10 -3 . We select the pair [START_REF] Huybrechts | The Antarctic ice sheet and environmental change: a three-dimensional modelling study= Der antarktische Eisschild und globale Umweltveränderungen: eine dreidimensionale Modellstudie[END_REF][START_REF] Greve | Comparison of Numerical Schemes for the Solution of the Ice-Thickness Equation in a Dynamic/Thermodynamic Ice-Sheet Model[END_REF] to display the misfits and estimations of diffusion. In Fig. 4.7 we plot the misfits and notice that indeed for the smaller value of k the misfit h is consistently lower for the four profiles. This is due to the difference in the quality of the estimates between experiments 2 and 4. Increasing k will produce smoother estimates but at the cost of losing detail. Profile D suffered the most as its distribution of D is the one with most variations. In fact, for any of the profiles, the regions of misfit directly correspond to those where the estimated diffusion deviated from D. Showing the diffusion profile to estimate (full black), the estimated profile D5 for experiment 5 (dotted light blue), D9 for experiment 9 (dashed orange), D10 for experiment 10 (dotted green), and D11 for experiment 1 (dash-dotted pink).

in s h, the error evolves slowly as it needs a large number of iterations before reaching a value close to the other three experiments. We saw in (4.11) that the derivative of the Lyapunov function has a term dependant on h2 and v 0 . So, the rate at which the error decreases will increase depending on the value of v 0 . This explains how the first dip in the error is due to the correction in the observer, while the rest of the behavior of the error is due to the update law. Experiments 9 and 10 go through the same phenomenon, though because of the lower values of v 0 , their errors evolve closer to that of experiment 5, with experiment 9 being the closest.

Towards the end of the simulations, during the relaxation phase, we see that the errors start increasing before they reach a steady-state. This increase diverges further from the steady-state of experiment 5, with higher values of v 0 . This happens because too much correction was done in the observer causing the update law (4.5) to be less effective. Looking at Fig. 4.13, this effect becomes more and more noticeable throughout the profiles and the experiments, until for profile D experiment 11, the estimate D11 is quite far from D. We notice that the effect of v 0 is not the same for all profiles.

Though it became evident that having this correction does not necessarily benefit the estimation of diffusion.

Conclusion of 1D Study Cases

In this section we tested a method [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] to estimate one-dimensional diffusion coefficients. We constructed four different profiles of 1D ice sheets by choosing distributions of basal sliding A s and calculating the corresponding distributions of diffusion D using (4.2). Estimating D implies the possibility of retrieving A s . The method consists of having an observer (4.4) to mimic the dynamics of the system (4.3) and an update law (4.5) to produce estimate of D.

We set-up experiments by varying the adaptation gain v 1 , the smoothing weight k, the regularization parameter α, and the correction gain v 0 . A summary of this set-up is shown in Table 4.1. Tests were done on four profiles, A, B, C, and D. The diffusion coefficient to be retrieved for each profile is shown in Fig. 4.2. Updates are done at fixed intervals of 100 steps. We saw that the used method can produce better results by increasing v 1 and lowering k. While increasing α will smooth parts of the estimates with higher slopes. This happens because we implemented a Tikhonov regularization [START_REF] Knowles | Methods for numerical differentiation of noisy data[END_REF] method to smooth the derivatives of the estimates. Varying k has a more global effect on the estimated diffusion coefficients. Setting v 0 > 0 did not provide any benefit to the estimation process, as it made the update law less effective, and thus produced more errors during steady-state.

Results of the 2D Study Cases

Having studied the one-dimensional case in the previous section, we now carry similar tests with two spatial dimensions. This section is kept brief as it mainly serves to show that the chosen method can be used in 2D, and that the conclusions of the 1D case still hold. We start by choosing a diffusive profile D and run both the system and observer in parallel. We then run a series of simulations by varying the correction gain v 0 and the smoothing weight k. The adaptation gain v 1 is kept constant as its effects are quite clear. And although v 0 had no positive impact on the estimation process in the 1D case, it is nonetheless tested here to confirm this point. Note that unlike the previous section, the smoothing of the estimated diffusion coefficients is done using an averaging filter. After having tried different 2D smoothing techniques, we noticed that the most straightforward and predictable method is a simple low-pass filter. Table 4.3 shows the set-up for the 2D experiments.

No. Profile Notice that the value of k is much larger than in the 1D case. For the 2D case, attempting to heavily lower k can introduce artifacts or push the simulations to be numerically unstable. The four profiles for diffusion are randomly generated, while keeping a degree of smoothness as we do not allow sharp changes in D. The spatial domain for the simulations is square. These profiles are shown in Fig 4 .14. The diffusive profiles are varied in complexity in order to test the capability of the method to retrieve details in D, within the proposed set-up. We also varied the bedrock b and the input a for each profile. We show these variables in Appendix A.

v 1 k v 0 1 X 10 -11 1
In experiment 1, no correction was done in the observer, and smoothing was set to its maximum with k = 1. Then in experiment 2, the smoothing weight was lowered to k = 0.5. And in experiment 3, we set k back to its maximum and introduce correction within the observer. We show in Fig. 4.16 the evolution of s h for all experiments and all profiles, and in Fig. 4.17 we compare the percentages of the absolute error in ice thickness and diffusion between experiments. It is abundantly clear now that the correction term is of no positive use for our purposes. In both 1D and 2D this term lowered the efficiency of the update law. Once again, when correction is stopped, during the end relaxation phase of the simulations, the error s h increased. While, in experiment 2, lowering the smoothing weight did benefit the estimation as both the errors in ice thickness and in diffusion are lower than in the other two experiments.

To better understand the results we look at the percentage errors in ice thickness and diffusion for experiment 2, in Fig. 4.18 and Fig. 4.19 respectively. We observe that the percentages in h are always lower than D. This means that the system (4.3) is not very sensitive to changes in the diffusion coefficient, and that more correction is needed in order to decrease h. This point becomes important when no knowledge of the real diffusive profile is available, as will be the case with real data. We also notice that most of the error is at the boundaries of the domain. This is to be expected as we have applied Dirichlet boundary conditions on our system, which makes the estimation harder due to the lack of error. In Fig. 4.15 we show the estimated diffusion coefficients D. They are in close resemblance to the real distributions. The increased complexity of the profiles was not a problem, though the smoothing did soften some of the sharper features of the profiles.

As a conclusion, we saw in the subsection results for the 2D study cases where four diffusive profiles were randomly generated and tested, each with a varied level of complexity. Results were similar to the 1D cases. The use of a correction term in the observer lowered the quality of the estimations as it produced more error. And lowering the weight of the smoothing weight was beneficial.

Though in order to maintain numerical stability, k was set much larger than in the 1D cases. Also, we decided to use a low-pass filter with a fixed window size to smooth the estimates. The estimated diffusion coefficients are in close agreement with their true values. Having found that the chosen 3), diffusion multiplies second derivatives of h, while derivatives of diffusion multiply first derivatives of h. This means that, for any given area of D, a high value implies a higher diffusion rate. While a higher variability in D implies more transport of ice. This is why regions of Antarctica subject to faster ice-flow, hence more transport of ice, exhibit distributions of D with both higher variations and amplitudes. The interior regions of the continent, where the ice is slower, have D with lower variations and lower amplitudes. The distribution of the error h appears to be higher for those same regions. This indicates that regions of fast ice flow must have diffusion coefficients with more variations and that further estimation is required.

Retrieving Basal Sliding

Our main goal in estimating the diffusion coefficient is to use it to find basal sliding. To do so we use (4.2), the estimated diffusion D calculated using the observer (4.4) and update law (4.5), and the measurements of ice thickness H and surface topography h. For a given simulation, at every update iteration, we propose the following simple strategy to estimate A s :

• Calculate D using (4.5).

• Calculate A s using :

A s = D -2 5 A ρ 3 g 3 H 5 h 2 x + h 2 y ρ 2 g 2 H 3 h 2 x + h 2 y = D v 2 - 2 5 A ρ g H 2 h 2 x + h 2 y (4.18) with v 2 = ρ 2 g 2 H 3 h 2 x + h 2 y
Note that in order to do so v2 needs to be different than zero. This term can be zero if either the ice thickness or the slopes of the surface topography are both zero, but in practice this is not a problem.

The models (system and observer) we run are based on solving (4.3) and (4.4) which rely on a diffusion coefficient, this is why in our case we could have simply calculated A s once at the end of the simulation. However, for models that depend on explicitly finding A s instead of D, the above strategy can be used. In such experimental set-ups, an ice-sheet model is to replace the observer.

We illustrate this strategy in Fig. 4.24.

In the previous sections, we set each new estimate of the diffusion coefficient to be a combination of a smoothed version and an unsmoothed version of it (4.17). For the 2D test cases and for real k is set to be bounded between two constants λ and 1 -λ, and is calculated as follows:

k (x, y, t) = | h | max | h | k (k < λ) = λ and k (k > 1 -λ) = 1 -λ D = (1 -k ) D + k D * (4.19)
where 0 ≤ λ ≤ 1.

The new parameter is used as a weight to unevenly smooth the new estimates of diffusion. Areas where error is high, thus k is high, are smoothed more, while areas where error is low, thus k is low, are smoothed less. This allows to fine tune areas where the error is low while still enforcing a satisfactory degree of smoothness in other areas. This ultimately produces less misfit and better estimates of basal sliding. The numerical set-up for the new experiments is shown in Table 4.5. We keep unchanged the adaptation gains v 0 = 0, v 1 = 10 -9 , the smoothing window w = 3 × 3, and the update step u s = 250. The time step is ∆t = 100, and the total number of updates is 500. S h is shown for each experiment.

In experiment 6, the smoothing weight k is kept constant and equal to 1, it is similar to experiment 1 in the previous section. We begin by showing the end part of the evolution of the mean absolute error S h in Fig. 4.25. We see that the best result was for experiment 7 with an average misfit of 45m. In general, using a variable k did bring some improvement to the results compared to experiment 6. Next we compare in Fig. 4.26 the distribution of misfits between experiments 6 and 7.

We see a noticeable improvement throughout the map, even though there is some increase in the error towards the center. We also notice, for either experiments, that the error branches out towards the boundary. This is due to the fact that the adopted shallow-ice approximation cannot reproduce Even though the used resolution is not high enough to show more details, the overall pattern fits with the one obtained in [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF][START_REF] Morlighem | Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model[END_REF]. Though in [START_REF] Morlighem | Inversion of basal friction in Antarctica using exact and incomplete adjoints of a higher-order model[END_REF], the inversion was for basal friction and not sliding, and of course more accurate physics were taken into account when doing the inversion. A very important point to keep in mind is that given any inverse problem, adopted simplifications and uncertainties in the measurements, modelling parameters, and model used, will be thrown on the variables to be estimated. We obtained relatively low misfits using a very simple diffusive model. These results can be improved by testing our method with dedicated ice-sheet models.

Conclusion

In this chapter we have used adaptive parameter identification to find estimates of the diffusion coefficient for 1D and 2D ice sheets.

• In Section 4.1 we reworked the proof for adaptive parameter identification [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] in two dimensions and with time-varying adaptation gains. The method introduces an observer (4.4) and an update law (4.5) to produce an estimate D of the diffusion coefficient D. We also introduced the possibility of having time-varying adaptation gains.

• Due to the stiffness of the problem at hand, regularization is needed to guarantee numerical stability. Another reason to seek regularization is that D needs to be C 1 over its spatial domain. This is why in Section 4.2 we proposed to smooth the estimates using Tikhonov regularization [START_REF] Nikolaevich | On the solution of ill-posed problems and the method of regularization[END_REF][START_REF] Knowles | Methods for numerical differentiation of noisy data[END_REF] for the 1D case. In this approach, we smooth the derivatives of the estimated diffusion. We then integrate back over the domain to obtain the regularized coefficient.

• In Section 4.3 we tested the estimation method on 1D ice sheets. Four different profiles were tested. These profiles were obtained using the coefficients of basal sliding from the previous chapter and (4.2). The obtained diffusion coefficients D are smoothed to avoid sharp changes in the coefficient. We thus obtain the coefficient D which is to be estimated. Then the system (4.3) and observer (4.4) are run together forward in time. We studied the sensitivity of the results to the update law adaptation gain v 1 in Section 4.3.1, to the smoothing weight in Section 4.3.2, to the regularization parameter α in Section 4.3.3, and to the observer correction gain v 0 in Section 4.3.4. We found that increasing v 1 improves the results. This occurs because the estimated diffusion is sensitive to changes in surface topography, while the reverse is not true. Meaning, large changes in D are needed to have relevant changes in h. Thus, during the updates, and with the decrease in error, a given value of v 1 becomes less effective with time.

We also found that lowering the smoothing weight k better captures detail in the coefficient, which produces less error. Adjusting the regularization gain α affects the smoothness of the estimates, but in a different way than k. Increasing α smooths areas of the diffusion coefficient with higher slopes. While k has a more global effect on the estimate. The correction gain v 0 has no positive effects on estimation. Introducing this parameter always produces less accurate coefficients. This happens because, in our problem, the dynamics of the observer and the system are identical. No hidden states exist, so there is no need for correction.

• In Section 4.4 we test the method in two dimensions. Four 2D profiles were tested. Unlike the 1D case, smoothing was done using a low-pass filter with a fixed window size. Each profile had a diffusion coefficient with varying degrees of complexity. Observations similar to the 1D case were obtained. Less smoothing and no correction produce better estimates.

elevation. We found that adopting a variable smoothing weight lowers the misfit. We also proposed a strategy to retrieve the basal sliding coefficient from the estimated diffusion coefficient. This strategy can be adapted to ice-sheet models by simply replacing the observer by the chosen model. The obtained A s showed an overall pattern that fits with other works.

In the next chapter we conclude this manuscript and end with propositions for future works.

the stability of the linear system (3.6)-(3.8) guaranteed, we applied the method to both the linear and nonlinear systems. This new formula can be used iteratively (3.29) to continuously generate better estimates of basal sliding during a simulation. We also introduced a smoothing effect through a weighted sum combining smoothed and unsmoothed versions of the estimate.

To test our method we generated four different profiles. The profile generation was done by choosing profiles of the bedrock, mass balance, and basal sliding, then running the nonlinear model forward in time, thus offering the distribution of ice thickness. This set of variables was considered to be the available measurements, at the exception of basal sliding. The estimation was done for the four profiles in both the linear and nonlinear cases. An ensemble of tests was created by varying the initial guess of basal sliding A s 0 , the update step u s , and the smoothing weight k. We also compared our method to another known method [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF] which uses misfits between simulated and measured ice thickness to update A s . We chose the other method to validate ours. If both methods produce relatively similar results, then our method is valid.

We found that our method is robust to initial guesses of basal sliding. The estimates were very close no matter the value of A s 0 , in both the linear and nonlinear cases. We then investigated the effect of changing the update step u s for the nonlinear case. Increasing u s gives the system more time to relax between updates, thus generating more error, which influence the estimation of A s .

The more error/information is available, the better are the estimates. Our method benefited from allowing longer relaxation times. As for changing the smoothing weight k, our method produced slightly better estimates with slightly higher k. This is because our method introduces stiff changes into the estimates, due to the form of the update law. Here we found that a trade-off exists between u s and k. When compared to [START_REF] Pollard | A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica[END_REF], our method generated similar results. All of these reasons clearly show that our new method is a viable option for estimating basal sliding for 1D ice sheets.

To end Chapter 3, we tested our method on real data. One dimensional flowlines had to be extracted from the map of surface topography of Antarctica. The continent can be divided into 27 different drainage basins [START_REF] Zwally | Antarctic and Greenland drainage systems[END_REF]. We selected a four different basins, set in three groups, and randomly chose ten flowlines for each group. We found that for a basin with a lot of lateral flows, the estimation technique did not work. However this is a limitation of the model used. The nonlinear model under consideration does not account for flows other than in one direction. So, if ice was to move sideways away from or into the 1D profile then a clear misfit will be seen between simulated and observed ice thickness. Saturation of A s occurred in regions where this phenomena is present. However, for successful simulations, we saw that A s will increase as we move away from the interior of the ice sheet. We know that ice in the inner parts of Antarctica is slower than near the edges of the continent. This was reflected in the estimated basal sliding coefficients.

was clear that the results improved when the window size decreased. The estimated distribution of diffusion revealed that for regions where the ice was fastest, not only did D have a higher amplitude but more spatial variation as well. This indicates, that for such regions, both diffusion and advection are dominant. While for the interior of the ice sheet the opposite happened. The estimate showed lower amplitudes and less variation. Most of the error in surface elevation was concentrated in regions of fast moving ice.

We propose the following perspectives as a continuation of this work.

Regarding the 1D inversion for basal sliding, saturation occurred in regions where fast ice flow could not be reproduced. This is why even when A s was at its lowest, the misfit was often negative.

Testing this method on more accurate flowline models can potentially solve the saturation issue in the estimates. This could push for a change in the Lyapunov function as more variables and parameters become involved in the dynamics, this is why an updated version of the update law should be investigated. Another interesting area of research is to test larger batches of flowlines in the various basins of Antarctica. Aggregating the obtained basal sliding coefficients could reveal interesting features for this parameter, especially that we saw that physically close flowlines generate somewhat similar profiles of basal sliding.

Upgrading our basal siding method from 1D to 2D was not achieved with the same Lyapunov function.

This is why we propose to continue the search for an appropriate candidate Lyapunov function that can lead the path to finding a new 2D update formula for A s . However, great care should be given to cross-terms that appear in the time derivative of the Lyapunov function. Getting rid of such terms is not a trivial task. A reformulation of the problem could be necessary.

In Chapter 4, and for the 1D case, we inverted for the diffusion coefficient. The set of diffusion profiles to be retrieved were obtained from the profiles of the test case in Chapter 3. This means that finding the estimated diffusion coefficients can directly lead to those of basal sliding. One interesting test is to estimate diffusion for Antarctic flowlines and then calculate the related profiles of basal sliding.

Less constraints are imposed when looking for diffusion, this is might avoid some of the saturation issues.

Moving to the 2D cases, we found that major factors in estimating diffusion were the used smoothing method and the smoothing weight. An important condition for the diffusion update law to work is that D ∈ C 1 (Ω), meaning that the distribution of diffusion must be smooth all over the spatial domain. This is why the estimates had to be filtered. It is very easy, though inadvertently, to introduce numerical stiffness into the simulations. Great care should be taken when estimating diffusion or basal sliding. The filtering techniques were effective but rather simple, this is why we propose that more elaborate smoothing methods should be investigated. When retrieving A s from D, we tested a variable smoothing weight dependant of the misfit h. This improved our results which suggests that adopting such a strategy can indeed lead to further improvements. We also saw that our method can be easily plugged into a dedicated ice-sheet model by simply replacing the observer with said model. It is only natural that the next step here is to do exactly that. Testing on a higher resolution is also a must. Combining these suggestions will surely produce distributions of basal sliding that reveal more accurate and interesting details. Models describing natural phenomena can depend on parameters that cannot be directly measured, hence the necessity to develop inverse techniques to determine them. Our goal is to utilize such techniques to enable better initialization of ice sheet models for Antarctica. This will help such models to produce better forecasts as part of climate studies. The parameters of interest are the basal sliding coefficient, which characterizes the contact of the ice sheet with the bed underneath, and the diffusion coefficient which dictates the dynamics within the mass-continuity partial differential equation describing the movement of ice sheets. A Lyapunov based approach is proposed to control the convergence of the 1D and 2D inhomogeneous transport models toward a feasible equilibrium matching the measurements of surface topography of the Antarctic ice sheet. Our work offers a new 1D update law for the basal sliding coefficient inversion. We also use adaptive distributed parameter inversion to retrieve basal sliding from diffusion in 1D and 2D models. These two methods are tested on study cases and real data. Our results show that the methods proposed are successful in inverting for sliding and diffusion while replicating the available data. Keywords: partial differential equations, Antarctic, basal sliding, diffusion, inhomogeneous transport phenomena, inverse problem.

viii 4 . 3 . 1 41 3. 2 F 2 GPSe 2 fforces acting on volume elements of ω kg m s - 2 Fof all forces acting on a body kg m s - 2 nP Total momentum of a body kg m s - 1 Qflux m 3 s -1 m - 2 U- 1 vV- 1 x

 431412222213211 23 2D, real data, h for all experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.24 2D, A s estimation block diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.25 2D, real data, S h for experiments 6-11. . . . . . . . . . . . . . . . . . . . . . . . . . 4.26 2D, real data, h experiments 6-7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 2D, real data, estimated basal sliding. . . . . . . . . . . . . . . . . . . . . . . . . . A.1 CH4, 2D, profile A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 CH4, 2D, profile B. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 CH4, 2D, profile C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 CH4, 2D, profile D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF TABLES Number Page Setup of linear model experiments. . . . . . . . . . . . . . . . . . . . . . . . . . . . Setup of non-linear model experiments. . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Non-linear case comparison guide. . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.4 Non-linear case, profile A, numerical results' summary. . . . . . . . . . . . . . . . . 70 3.5 Non-linear case, profile B, numerical results' summary. . . . . . . . . . . . . . . . . 71 3.6 Non-linear case, profile C, numerical results' summary. . . . . . . . . . . . . . . . . 71 3.7 Non-linear case, profile D, numerical results' summary. . . . . . . . . . . . . . . . . 72 4.1 1D, experimental set-up. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.2 1D, summary of results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.3 2D, experimental set-up. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.4 2D, real-data, experimental set-up. . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.5 2D, real-data, basal sliding experiments. . . . . . . . . . . . . . . . . . . . . . . . Initial value of basal sliding coefficientm a -1 Pa -1 Ās Estimated basal sliding coefficient m a -1 Pa -1 Ãs Error in estimated basal sliding coefficient m a -1 Pa -1 Flux of G through the boundary ∂ω -Some additive physical quantity on the entirety of ω Ice-sheet surface elevation in observer system m ĥ0 Initial value of ice-sheet surface elevation in observer system m Internal production of G within ω s h Sum of the absolute error in surface elevation m S h Mean absolute error in surface elevation m External supply of G outside of ω -Internal forces (stresses) acting on ∂ω kg m s -External Sum Unit normal vector to ∂ω -Volume Velocity of ice on base or surface m s Velocity field vector on ∂ω m s -1 xii Velocity of ice base or ice surface m s

Figure 1 . 2 :

 12 Figure 1.2: Diagram of an ice-sheet showing how the ice flows towards the ocean. The ice-shelf is the mass of ice floating on water. Seasonal changes can either cause the thickness of the ice-sheet to increase or decrease. High pressures at the bottom of the ice-sheet lower the melting point of ice causing it to melt.

Figure 1 . 3 :

 13 Figure 1.3: Sketch of a vertical windmill, used for grinding grain, found in ancient Persia (modern day Iran).

Figure 2 . 1 :

 21 Figure 2.1: A material volume ω with ∂ω as its boundary, n is a unit normal vector to the surface, v is the velocity field, x, y and z are the orthonormal basis vectors.

Theorem 1 .

 1 The temporal change of an arbitrary field quantity Γ over a material volume ω is composed of two parts: 1. the local change ∂Γ/∂t within ω 2. the advective flux Γv in the normal direction n across ∂ω

Figure 2 . 2 :e x x e x y e xz e yx e yy e zy e zx e zy e zz

 22 Figure 2.2: A cross section of a material with volume ω, surface ∂ω, internal forces e, and external forces f.

Figure 2 . 3 :

 23 Figure 2.3: The Cauchy stress tensor.

Figure 2 Figure 2 . 4 :

 224 Figure 2.4: Cross section of an ice sheet with ice-shelf extending towards the ocean. The top surface

25 )Figure 2 . 5 :

 2525 Figure 2.5: Illustrations of the top (left) and bottom (right) surfaces of an ice sheet, showing the velocity vector U of the ice on the top or base surfaces, the velocity V of the top or base surfaces, and the unit normal vector n.

3 . 10 )

 310 ∀ t ∈ [0, T] along the solutions of (3.6)-(3.8) with the transport coefficients (3.7).

Figure 3 . 1 :

 31 Figure 3.1: Block diagram of the numerical experiment used to evaluate the estimation method on fictitious data. The blocks in gray represent phase 1 (synthesis of fictitious data), while those in white represent phase 2 (estimation).

Figure 3 . 3 :

 33 Figure 3.3: 1D linear system, profile B: reference variables. On the left, the ice thickness H (dashed blue), surface elevation h (full pink), bedrock b (full brown). On the top right, the surface mass balance a. On the bottom right, the basal sliding coefficient Ās .

Figure 3 . 4 : 3 0

 343 Figure 3.4: 1D linear system, profile A: errors. On the left, the evolution of the cumulative absolute error s H , and on the right the misfit H, for A s 1 0 = 10 -5 (full green), A s 2 0 = 10 -6 (dashed yellow), A s 3 0 = 10 -8 (dotted purple).

Figure 3 . 5 : 1 0 2 0

 3512 Figure 3.5: 1D linear system, profile B: errors. On the left, the evolution of the cumulative absolute error s H , and on the right the misfit H, for A s 1 0 = 10 -5 (full green), A s 2 0 = 10 -6 (dashed yellow), A s

Figure 3 . 7 :

 37 Figure 3.7: 1D linear system, profile D: reference variables. On the left, the ice thickness H (dashed blue), surface elevation h (full pink), bedrock b (full brown). On the top right, the surface mass balance a. On the bottom right, the basal sliding coefficient Ās .

Figure 3 . 8 : 1 0 2 0 3 0

 38123 Figure 3.8: 1D linear system, profile C: errors. On the left, the evolution of the cumulative absolute error s H , and on the right the misfit H, for A s 1 0 = 10 -5 (full green), A s 2 0 = 10 -6 (dashed yellow), A s 3 0 = 10 -8 (dotted purple).

Figure 3 . 9 : 3 0

 393 Figure 3.9: 1D linear system, profile D: errors. On the left, the evolution of the cumulative absolute error s H , and on the right the misfit H, for A s 1 0 = 10 -5 (full green), A s 2 0 = 10 -6 (dashed yellow), A s 3 0 = 10 -8 (dotted purple).

Figure

  Figure3.10: 1D linear system, estimations of basal sliding for all profiles: log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference (full pink).

Figure 3 . 11 : 1 0 2 0

 31112 Figure 3.11: 1D non-linear system, all profiles, experiments 1,2,3,4: evolution of the cumulative absolute error s H . Method 1: A s 1 0 (full dark blue) and A s 2 0 (full blue), Method 2: A s 1 0 (dashed red)

Figure 3 . 12 : 1 0 2 0

 31212 Figure 3.12: 1D non-linear system, all profiles, experiments 1-2: estimations of basal sliding. log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference for A s 1 0 (full pink) and A s 2 0 (dotted blue).

Figure 3 .s 1 0 2 0

 312 Figure 3.13: 1D non-linear system, all profiles, experiments 3-4: estimations of basal sliding. log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference for A s 1 0 (full pink) and A s 2 0 (dotted blue).

Figure 3 . 14 : 2 0 1 0 2 0 53 Figure 3 . 15 :

 31421253315 Figure 3.14: 1D non-linear system, all profiles, experiments 1,2,3,4: misfit H. Method 1: A s 1 0 (full

Figure 3 .

 3 Figure 3.16: 1D non-linear system, all profiles, experiments 2,4,5,6: misfit H. Method 1: u 1 s (full

Figure 3 .Figures 3 .

 33 Figure 3.18: 1D non-linear system, all profiles, experiments 6-4: estimations of basal sliding. log 10 of the reference coefficient Ās (full blue), the estimate A s (dashed brown), and their absolute difference for u 1 s (full pink) and u 2 s (dotted blue).

  Figure 3.19: 1D non-linear system, all profiles, experiments 2,4,7,8,9,10: evolution of the cumulative absolute error s H . Method 1: k 1 (full dark blue), k 2 (full blue), and k 3 (full grey), Method 2: k 1 (dashed red) and k 2 (dashed purple), and k 3 (dashed brown).

Figure 3 . 20 : 57 Figure 3 . 21 :

 32057321 Figure 3.20: 1D non-linear system, all profiles, experiments 7,8,9,10: misfit H. Method 1: k 2 (full dark blue) and k 3 (full blue), Method 2: k 2 (dashed red) and k 3 (dashed purple).

Figure 3 . 22 :

 322 Figure 3.22: 1D non-linear system, all profiles, experiments 4-9-10: estimations of basal sliding.

Figure 3 .

 3 Figure 3.23: 1D non-linear system, all profiles, experiments 3-11-12: evolution of the cumulative absolute error s H . Method 2: H 1 inv (full red) and H 2 inv (dotted purple), and H 3 inv (dashed brown).

Figure 3 . 24 :

 324 Figure 3.24: Antarctic drainage systems numbered from 1 to 27 and each represented by a different color. The grounding line is highlighted by a white line.

Figure 3 . 25 :

 325 Figure 3.25: Antarctic basin 3, the chosen flowlines and their relative mean absolute errors. The threshold of 20 m is shown in dotted green.

Figure 3 . 26 :

 326 Figure 3.26: Antarctic basin 3, flowline 5.In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

Figure 3 . 27 :

 327 Figure 3.27: Antarctic basin 3, flowline 10. In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

Figure 3 .

 3 Figure 3.28: Antarctic basin 10, the chosen flowlines and their relative mean absolute errors. The threshold of 20 m is shown in dotted green.

Figure 3 . 29 :

 329 Figure 3.29: Antarctic basin 10, flowline 1.In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

Figure 3 . 30 :

 330 Figure 3.30: Antarctic basin 10, flowline 4. In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

Figure 3 . 33 :

 333 Figure 3.33: Antarctic basin 17, flowline 8. In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

Figure 3 . 34 :

 334 Figure 3.34: Antarctic basin 21/22, the chosen flowlines and their relative mean absolute errors. The threshold of 20 m is shown in dotted green.

Figure 3 . 35 :

 335 Figure 3.35: Antarctic basin 21/22, flowline 1.In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

Figure 3 . 36 :

 336 Figure 3.36: Antarctic basin 21/22, flowline 2.In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

  ht = ĥth t = D ĥx + hx D x + D ĥy + hy D y v 0 h (4.8)

(4. 10 ) 8 .

 108 TheoremThe time derivative V t of the function V given by (4.10) verifies:

( 4 . 2 )

 42 , the presence of the absolute value of the gradient of h can cause sharp transitions in D, which can lead to numerical complications. In Fig.4.2 we show D and D, the distributions of diffusion before and after smoothing respectively, for the four profiles A, B, C and D. We tackle the case of on-line measurements. This means that both the system (4.3) and the observer (4.4) are running in parallel forward in time. We show in Fig.4.1 bar graphs for the mean of the absolute error in surface topography hm , and the sum of the absolute error in surface topography hss at steady-state, for all experiments and all profiles. It is clear that experiments 5 and 6 have the lowest errors.

Figure 4 . 1 :

 41 Figure 4.1: 1D case, the mean of the absolute error in surface topography hm (blue) and the sum of the absolute error in surface topography hss (red) at steady-state, for experiments 1-8 and all profiles.

Figure 4 . 2 :

 42 Figure 4.2: 1D case, normalized distribution of diffusion for all profiles, before smoothing ( D, blue) and after smoothing (D, dashed orange).

Figure 4 . 3 :

 43 Figure 4.3: 1D case, comparison of experiments based on value of v 1 . The evolution of s h, the sum of the absolute misfit in h, for all profiles.

Figure 4 . 4 :

 44 Figure 4.4: 1D case, comparison of misfits based on value of v 1 for all profiles.The surface topography h (light blue), bedrock b (dark blue) correspond to the left axis. The absolute misfits in surface topography h1 for experiment 1 (dotted orange) and h5 for experiment 5 (dashed dark orange) cor- respond to the right axis.

Figure 4 . 5 :

 45 Figure 4.5: 1D case, comparison of estimated diffusion coefficients based on the value of v 1 for all profiles.Showing the diffusion profile to estimate (full orange), the estimated profile D1 for experi- ment 1 (dotted red), and the estimated profile D5 for experiment 5 (dashed blue).

Figure 4 . 6 :

 46 Figure 4.6: 1D case, comparison of experiments based on value of k. The evolution of s h, the sum of the absolute misfit in h, for all profiles.

Figure 4 . 7 :

 47 Figure 4.7: 1D case, comparison of misfits based on value of k for all profiles.The surface topography h (light blue), bedrock b (dark blue) correspond to the left axis. The absolute misfits in surface topography h2 for experiment 2 (dotted orange) and h4 for experiment 4 (dashed dark orange) cor- respond to the right axis.

Figure 4 . 8 :

 48 Figure 4.8: 1D case, comparison of estimated diffusion coefficients based on the value of k for all profiles.Showing the diffusion profile to estimate (full orange), the estimated profile D2 for experi- ment 2 (dashed blue), and the estimated profile D4 for experiment 4 (dotted red).

Figure 4 . 9 :

 49 Figure 4.9: 1D case, comparison of experiments based on value of α. The evolution of s h, the sum of the absolute misfit in h, for all profiles.

Figure 4 . 10 :

 410 Figure 4.10: 1D case, comparison of misfits based on value of α for all profiles.The surface topography h (light blue), bedrock b (dark blue) correspond to the left axis. The absolute misfits in surface topography h7 for experiment 7 (dotted orange) and h8 for experiment 8 (dashed dark orange) cor- respond to the right axis.

Figure 4 .

 4 Figure 4.11: 1D case, comparison of estimated diffusion coefficients based on the value of α for all profiles.Showing the diffusion profile to estimate (full orange), the estimated profile D7 for experi- ment 7 (dashed blue), and the estimated profile D8 for experiment 8 (dotted red).

Figure 4 .

 4 Figure 4.12: 1D case, comparison of experiments 5 (full light blue), 9 (dashed orange), 10 (dotted green) and 11 (dash-dotted pink). The evolution of s h, the sum of the absolute misfit in h, for all profiles. Zoom on figure for iteration i > 1e4.

Figure 4 .

 4 Figure 4.13: 1D case, comparison of estimated diffusion coefficients based on the value of v 0 for all profiles.Showing the diffusion profile to estimate (full black), the estimated profile D5 for experiment 5 (dotted light blue), D9 for experiment 9 (dashed orange), D10 for experiment 10 (dotted green), and D11 for experiment 1 (dash-dotted pink).

Figure 4 .

 4 Figure 4.14: 2D case, distribution of diffusion D for all profiles.

Figure 4 .

 4 Figure 4.15: 2D case, distribution of estimated diffusion D for all profiles.

Figure 4 .

 4 Figure 4.16: 2D case, the evolution of s h, the sum of the absolute misfit in h, for all profiles. Experiment 1 in full light blue, experiment 2 in dashed orange, and experiment 3 in dotted green.

Figure 4 . 17 :

 417 Figure 4.17: 2D case, the percentage of mean absolute error in surface topography (blue) and the percentage of mean absolute error in diffusion coefficient (red) at steady-state, for experiments 1-3 and all profiles.

Figure 4 .

 4 Figure 4.18: 2D case, experiment 2, percentage misfit in ice thickness h for all profiles.

Figure 4 .

 4 Figure 4.19: 2D case, experiment 2, percentage misfit in diffusion coefficient D for all profiles.

Figure 4 .

 4 Figure 4.20: 2D case, real data. Mask of grounded ice sheet is on the top left, surface topography h on the top right, bedrock elevation b on the bottom left, and the mass balance a on the bottom right.

Figure 4 . 21 :

 421 Figure 4.21: 2D case, real data. The evolution of s h, the sum of the absolute misfit in h, for experiments 1 (full light blue), 2 (dashed orange), 3 (dotted green), 4 (dash-dotted pink), and 5 (full violet).

Figure 4 .

 4 Figure 4.22: 2D case, real data. Distribution of estimated diffusion D for all experiments.

Figure 4 .Fig. 4 . 23 .

 4423 Figure 4.23: 2D case, real data. Misfit in surface elevation h for all experiments.

Figure 4 . 24 :

 424 Figure 4.24: Block diagram to estimate basal sliding.

Figure 4 . 25 :

 425 Figure 4.25: 2D case, real data. The evolution of S h, the mean absolute misfit in h, for experiments 6 (blue), 7 (orange), 8 (green), 9 (gray), 10 (purple), and 11 (red).

Figure 4 . 26 :

 426 Figure 4.26: 2D case, real data. Misfit in surface elevation h for experiments 6 and 7.

Figure 4 . 27 :

 427 Figure 4.27: 2D case, real data. Estimated basal sliding from experiments 7.

Figure A. 2 :

 2 Figure A.2: 2D case, profile B. Diffusion coefficient D is on the top left, surface topography h on the top right, bedrock elevation B on the bottom left, and the mass balance a on the bottom right.

Figure A. 3 :

 3 Figure A.3: 2D case, profile C. Diffusion coefficient D is on the top left, surface topography h on the top right, bedrock elevation B on the bottom left, and the mass balance a on the bottom right.

Figure A. 4 :

 4 Figure A.4: 2D case, profile D. Diffusion coefficient D is on the top left, surface topography h on the top right, bedrock elevation B on the bottom left, and the mass balance a on the bottom right.

  to approach zero uniformly, it is not enough for H to go to zero, but it requires the norm of ∆ to approach zero faster than the norm of H for all time t:

	lim H →0	sup t≥0	∆(H, t) H	= 0	(2.56)
	Lyapunov's indirect method can now be stated as:
	Theorem 5. Consider a non-linear system as in (2.50) and assume (2.56) to be true. Also, let A(t)
	as in (2.53) be bounded. If 0 is a uniformly asymptotically stable equilibrium point of H L (2.54), then

H L (2.55)

For ∆(H, t)

Table 3 .

 3

1: Setup of linear model experiments.

Table 3 .

 3 2: Setup of non-linear model experiments.A total of 48 simulations are done as seen in Table3.2. Experiments 1 and 2 are done on Method 1 to test the robustness to two different initial guesses, the same is done in experiments 3 and 4 for Method 2. In 5 and 6 the effect of the update step is tested on both methods. In 7-8 and 9-10 we study the effect of changing the smoothing weight k on Method 1 and 2 respectively. Experiments 11 and 12 focus only on Method 2 by varying the scaling factor H

inv . As these tests have many changing variables, we name the variables as follows: A s 1 0

Table 3 .

 3 Comparison Method 1, smoothing weight : k 1 , k 2 and k 3 3, 9 and 10 Method 2, smoothing weight : k 1 , k 2 and k 3 3: Comparison guide for the non-linear experiments.

	1 and 2		Method 1, initial guesses : A s	1 0 and A s	2 0
	3 and 4		Method 2, initial guesses : A s	1 0 and A s	2 0
	1 and 3		Method 1 and Method 2
	2 and 4		Method 1 and Method 2
	5 and 6		Method 1 and Method 2
	2 and 5		Method 1, update step : u 1 s and u 2 s
	4 and 6		Method 2, update step : u 1 s and u 2 s
	1, 7 and 8
	Profile s A No. Method 1 1 8.84 × 10 -2 2 1 4.03 × 10 -2 3 2 3.69 × 10 -1 4 2 1.34 5 1 1.83 × 10 -1 6 2 6.21 7 1 8.75 × 10 -2	3.67 × 10 -5 1.76 × 10 -5 1.52 × 10 -4 5.73 × 10 -4 8.03 × 10 -5 2.8 × 10 -3 3.54 × 10 -5	1.22 × 10 -2 2.4 × 10 -3 2.52 × 10 -2 6.66 × 10 -2 9.8 × 10 -3 2.18 × 10 -1 5.6 × 10 -3
	8	1	1.05 × 10 -2	4.21 × 10 -6	6.88 × 10 -4
	9	2	2.88 × 10 -2	1.16 × 10 -5	1.8 × 10 -3
	10	2	5.05 × 10 -2	2.06 × 10 -5	2.5 × 10 -3
	11	2	5.32 × 10 -1	2.21 × 10 -4	5.64 × 10 -2
	12	2	3.22 × 10 -1	1.35 × 10 -4	2.16 × 10 -2

H Average % | H | Average % | Ãs |

Table 3 .

 3 4: Non-linear system, numerical results' summary for profile A, all experiments. Showing the steady-state error s H , the average percentage absolute errors in ice thickness and estimated

	sliding coefficient.

Table 3 .

 3 5: Non-linear system, numerical results' summary for profile B, all experiments. Showing the steady-state error s H , the average percentage absolute errors in ice thickness and estimated sliding coefficient.

	Profile	No.	Method	s H	Average % | H |	Average % | Ãs |
		1	1	4.73	2.1 × 10 -3	3.31
		2	1	4.38	2 × 10 -3	3.20
		3	2	4.35	2 × 10 -3	4.7 × 10 -1
		4	2	4.18	1.9 × 10 -3	4.8 × 10 -1
		5	1	11.54	5 × 10 -3	1.5
	C	6 7	2 1	26.6 8.36	1.23 × 10 -2 3.4 × 10 -3	1.18 5.23 × 10 -1
		8	1	16.45	6.5 × 10 -3	5.24 × 10 -1
		9	2	20.21	8.8 × 10 -3	4.1 × 10 -1
		10	2	40.35	1.76 × 10 -2	7.5 × 10 -1
		11	2	2.95	1.5 × 10 -3	1.46
		12	2	8.08	3.6 × 10 -3	2.6 × 10 -1

Table 3 .

 3 6: Non-linear system, numerical results' summary for profile C, all experiments. Showing the steady-state error s H , the average percentage absolute errors in ice thickness and estimated sliding coefficient.

	Profile	No.	Method	s H	Average % | H |	Average % | Ãs |
		1	1	3.33	1.6 × 10 -3	7.71 × 10 -1
		2	1	2.96	1.4 × 10 -3	7.71 × 10 -1
		3	2	10.34	5.4 × 10 -3	9.2 × 10 -1
		4	2	9.87	5.2 × 10 -3	9.4 × 10 -1
		5	1	6.03	2.8 × 10 -3	1.14 × 10 -1
	D	6 7	2 1	33.16 15.03	1.83 × 10 -2 7.3 × 10 -3	2.48 1.3
		8	1	29.93	1.44 × 10 -2	2.21
		9	2	44.92	2.37 × 10 -2	2.75
		10	2	86.56	4.54 × 10 -2	5.49
		11	2	4.72	2.5 × 10 -3	6.8 × 10 -1
		12	2	21.40	1.11 × 10 -2	1.47

Table 3 .

 3 7: Non-linear system, numerical results' summary for profile D, all experiments. Showing

Table 4 .

 4 1: Setup of 1D diffusion experiments. Showing the adaptation gains v 0 and v 1 , the smoothing weight k, and the regularization parameter α.

	No.	A	B	C	D
	1	3.44	3.16	3.68	5.58
	2	3.18	2.60	2.80	8.15
	3	33.00	39.28	46.62	29.83
	4	31.39	27.04	31.35	39.92
	5	0.36	0.36	0.46	0.76
	6	0.37	0.31	0.41	2.05
	7	3.47	3.58	4.36	4.69
	8	3.22	2.55	2.85	5.52
	9	0.36	0.37	0.47	0.76
	10	1.02	1.07	1.30	0.92
	11	8.16	23.52	129.32	502.87

Table 4 .

 4 2: Mean absolute error in m for 1D diffusion experiments.

Table 4 .

 4

3: Setup of 2D diffusion experiments. X can be one of three profiles: A, B, C or D.

Table 4 .

 4 5: Setup of real data experiments with variable smoothing weight. The mean absolute error

=

-8 (dotted purple).

ACKNOWLEDGEMENTS

I tend to talk too much, so I will keep this brief:

To friends, family, and all those involved in this PhD

In the upper left, the reference ice thickness H, the estimated ice thickness H, reference surface elevation h, estimated surface elevation h, and bedrock b. In the upper right, the misfit H. In the lower left, the cumulative absolute error s H . In the lower right, the estimated basal sliding A s .

C h a p t e r 4 ESTIMATING THE DIFFUSION COEFFICIENT IN 1D AND 2D MODELS

In the previous chapter we described how we found and tested a new method for estimating basal sliding in a 1D ice sheet. The natural follow-up was to expand the method into two dimensions (2D).

However, with our choice of Lyapunov function, and following the same methodology, no new feasible update formula was found. Two difficulties surfaced. First, using the same Lyapunov function, and duplicating the update formula (3.27) by adding a new dimension, introduced new terms into V t .

We were unable to get rid of those terms thus the conditions to prove convergence and stability were not met. Second, seeking a new Lyapunov function that fits our needs then deduce a suitable update law for A s proved to be an arduous task in 2D. This is why our perspective had to change. A more practical approach is to look for the diffusion coefficient. We saw in chapter 2 that, under SIA, the diffusion coefficient D depends on the sliding coefficient A s (2.49). In this chapter we seek to estimate D. The diffusion coefficient is estimated using an observer and an update law, based on adaptive distributed parameter identification. Study tests for the one and two dimensional cases are done on fictitious data to illustrate the performance of the method. And like the previous chapter, more tests are done on real measurements of surface topography of the Antarctic ice sheet.

In this chapter we:

Introduce the method of adaptive parameter identification Justify why data regularization is necessary Test adaptive parameter identification numerically Present the results for the 1D and 2D test cases Present the results for real data Retrieve basal sliding from diffusion in 2D

Sensitivity to the Regularization Parameter

In this subsection, we turn our focus to the regularization parameter α. Note that this parameter first appeared in the functional (4.15) which is to be minimized in the course of applying regularization.

In our case, derivatives of D are calculated, by the means of Tikhonov regularization [START_REF] Knowles | Methods for numerical differentiation of noisy data[END_REF], then integrated to obtain smoothed estimates. Giving more weight to α means adding more penalty on variations in the derivative. This implies that higher values of α will introduce more smoothing into the estimates. However, as we will explain next, the effect is not exactly similar to increasing k. We group the experiments into the following pairs: (1,2), [START_REF] Evans | Partial differential equations[END_REF][START_REF] Greve | Comparison of Numerical Schemes for the Solution of the Ice-Thickness Equation in a Dynamic/Thermodynamic Ice-Sheet Model[END_REF], [START_REF] Hutter | Thermo-mechanically coupled ice-sheet response -cold, polythermal, temperate[END_REF][START_REF] Weertman | On the Sliding of Glaciers[END_REF], and [START_REF] Ribe | Coiling of viscous jets[END_REF][START_REF] Ribe | Stability of liquid rope coiling[END_REF].

In Fig. 4.9 we plot the sum of absolute misfits s h and focus on the iterations closer to the end of the simulations (i ≥ 4000). We notice that, unlike the two previous subsections, the experiments cannot be separated into two groups based on the value of α alone. Similarly looking at | h |, in Fig. 4.10, does not shed light on how α affects the experiments. The error is distributed similarly across the profiles for the two values of α. However, when looking at the estimated diffusion, in Fig. 4.11, we see that profile D was affected the most. We know that α penalises variations in the derivatives of D. And from the four tested profiles, D is the one with the most variations, thus its derivatives will be highly varying as well, this is why the increase in the regularization parameter affected it the most.

More so, it is the region of D with the higher slopes that was smoothed the most. While for the other three profiles the effect is less noticeable.

Sensitivity to the Correction Gain

In this subsection we study the effect of the observer correction gain v 0 . In experiments 9-11, v 0 was increased from 10 -7 to 10 -5 then 10 -3 , while the rest of the variables were set as in experiment 5.

We chose experiment 5 as it is the one with the best results, thus it offers a good benchmark to evaluate the effects of v 0 .

The main drawback of having v 0 not equal to zero is the fact that part of the error, between the observer and the system, will decrease independently of the estimated diffusion coefficient. Meaning, that if no update was done and v 0 > 0, then s h can drop. Remembering that the update law (4.5) is dependent on h, this implies that the misfit needed to produce good estimates of D will be absent.

Thus causing the update law to be less efficient.

In Fig. 4.12 we show the evolution of s h for the selected four experiments. We also enlarge part of the plots to focus on the end of the simulations. We do so because we have allowed some relaxation period in order to see if changes might occur once correction within the observer and update of the diffusion coefficient are stopped. We quickly notice that for experiment 11, after a quick initial drop method if feasible to estimate diffusion in 2D, and having studied the sensitivity of the method to some of the simulation variables, we can move to estimating the distribution of the diffusion coefficient for Antarctica.

Diffusion for Real Data

Estimating the diffusion coefficient using real data [START_REF] Fretwell | Bedmap2: improved ice bed, surface and thickness datasets for Antarctica[END_REF] is done in this section. We carry out the estimation using adaptive distributed parameter identification [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF]. This method was successfully tested on 1D and 2D ice sheets in Sections 4.3 and 4.4. However, unlike the two previous sections, the system 4.3 is considered to be in steady-state. The real data describe this steady-state. Hence, we have distributions of the surface elevation h, the bedrock elevation b, the ice-thickness H, and the input a. The estimated diffusion coefficient D is calculated using these variables, the observer (4.4), and the update law (4.5). Of course, the real diffusion coefficient D is not available.

Because we have already studied the sensitivity of the method to some of the simulation variables, we do not repeat all the tests. Instead, we fix the update gain v 1 and set the correction gain v 0 = 0.

In 2D, a very influential factor is the smoothing of the estimates. In Section 4.4, a 2D low-pass filter of fixed size w was used to smooth D. In this section, we investigate the effects of varying the size of the filter. We saw that lowering the smoothing weight k is beneficial for the estimation, this is why we illustrate this effect again in this section. The numerical set-up of the simulations is shown in Table 4.4. The updates are fixed at every 250 iterations with a time step of 100 iterations. The resolution of the data is 20 km. In Fig. 4.20 we show the system variables, we focus on the grounded ice sheet.

No. 4.4. We can clearly see that the error decreases when using a smaller window size for the filter, and with a lower smoothing weight k. During these experiments, we did try lowering k further, but this introduced artifacts and caused numerical problems. In addition, all the numerical derivatives within the observer and the update law were taken as Gaussian derivatives in order to avoid stiff changes.

C h a p t e r 5

CONCLUSION AND PERSPECTIVES

This thesis represents an interdisciplinary effort merging control engineering and glaciology. The aim was to estimate two spatially distributed parameters pertinent to ice sheet models. The first parameter is the basal sliding coefficient A s . It describes how ice sheets slide over the underlying bedrock. The second is the diffusion coefficient D. It is of a more abstract nature though its distribution heavily influences the dynamics of ice sheets. This parameter describes how ice tries to spread evenly over its spatial domain, and how ice will be transported from one area to another.

These two parameters are dependant, as the distribution of one can imply the other. This kind on problem is generally known as an inversion process, because the estimation is done using available measurements, as opposed to forward problems where a sets of variables are used to project a model forward. The ultimate goal was to reproduce measurements of the surface topography of the Antarctic ice sheet.

In Chapter 2, we offered the necessary theoretical background. We begun with a review of the physics influencing the movement of grounded ice sheets. The conservation laws of mass and momentum, along with the boundary conditions at the top and bottom of ice sheets, and the shallow ice approximation, allowed us to obtain the main dynamical equation, known as the ice thickness equation in diffusive form (2.49). We also defined the equation relating diffusion to basal sliding.

We then followed by stating the formal definitions and theorems from Lyapunov theory. This branch of control is aims at studying the stability of systems by investigating the stability of some measure of the internal energy of that system. In our case we set the integral (over the spatial domain) of the squared difference between the simulated and measured ice thickness to be the measure of stability. Through the proposed inversions, if this measure tends to zero, then the simulations match the measurements, and the estimations of basal sliding and diffusion are considered valid. We ended the chapter with a review of relevant works in parameter estimation in both glaciology and control.

We offer two contributions laid out in chapters 3 and 4.

In Chapter 3, we proposed a new method to estimate basal sliding for one dimensional ice sheets.

We started by linearizing the nonlinear dynamics (3.3) around an equilibrium, represented by the available measurements of surface topography. We showed that the chosen Lapunov function (3.9) is proven to be always decreasing with time if the update of basal sliding is chosen as in (3.27). With

In Chapter 4, we use adaptive distributed parameter system identification [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF] to find estimates of the diffusion coefficient D in both 1D and 2D cases, and for study cases and real data. We started by defining our problem and defining the system (4.3), the observer (4.4), and the update law (4.5).

The goal here is to find estimates of diffusion D, such that the misfit between the observer and the system is zero. Two new variables are introduced. The first is the correction gain v 0 . It serves as a gain to reduce the difference between the observer surface topography ĥ and the system's h. The second is the adaptation gain v 1 . This gain influences the weight of the corrections done to D. In other words, the influence of each update. We lay out the proof for the Lyapunov theory as done in [START_REF] Orlov | Adaptive distributed parameter systems identification with enforceable identifiability conditions and reduced-order spatial differentiation[END_REF], but in two spatial dimensions and time varying gains v 0 and v 1 . Regularization was introduced for both the 1D and 2D cases. Though for 1D we used Tikhonov regularization [START_REF] Knowles | Methods for numerical differentiation of noisy data[END_REF] to smooth variations in D. This introduced a the regularization parameter α. And like the previous chapter, smoothing is a weighted sum between smoothed and unsmoothed versions of D. The smoothing weight is denoted by the variable k.

For the 1D study cases, we selected the same four profiles as in Chapter 3. The basal sliding coefficients, bedrock and surface topography profiles were used to calculate diffusion. However, we used an averaging filter to remove sharp changes in D and avoid numerical problems during the simulations. For the tests we varied v 1 , k, α, and v 0 . We found that lowering v 1 produced less misfit for all profiles, and that the diffusion profiles were reproduced with good fidelity. Increasing k caused the estimates to lose detail, as expected. This of course caused more error. Increasing α had similar effects, though regions where the slopes of D are highest were smoothed the most, as opposed the more global effect k has. This is due to the employed Tikhonov regularization method. The correction gain v 0 had no positive impact. Because its presence inherently reduces the misfit h, the more v 0 increases the less effective is the update law. At the end of each simulation a relaxation period, where v 0 and v 1 are set to zero, shows if the steady-state error remains constant. This was not the case when v 0 was used.

For the 2D study cases, a similar logic was followed. Four diffusive profiles were randomly generated and were to be estimated. We did not repeat the whole ensemble of tests. Instead we focused on the effects of changing k and v 0 . The same observations as the 1D case were noted. Increasing k removed some of the finer variations in the diffusion coefficient. And including correction in the observer generated more errors. Looking at the distribution of error in both D and ĥ, we noticed that though low, it was mostly concentrated at the boundaries of the domain, and at points were the original coefficient D varied the most.

We then moved to the real data. Having found that smoothing is an important factor in the estimation process, we focused solely on varying the window size of the low pass filter used to smooth D. It

A p p e n d i x A

2D STUDY PROFILES IN CHAPTER 4

Here we show the four profiles used for the 2D study cases in chapter 4. The surface topography h is the steady-state obtained from running the system (4.3) forward in time.