Dr Fabien Gillet-Chaulet

Prof Adrian Luckman

Amber Nathan Laura

Jai

Julian

Julian, Jai, Pete Jordi Nathan

Juan Pe

Robin

Wei Marc

Moham- Mad Karol

Lucas

David

Martina Chiara

Climate Change has already caused large losses on natural occurring ice in the last decades and continues to do so. Antarctica, the largest source of land ice, contains an equivalent of several tens of meters of sea water level rise, which is induced not only by melting but also by out flow into ice shelves, the floating tongues of glaciers. In the last decades with recorded longer melt seasons and warming water in Antarctica, ice shelves are thinning and breaking up. With that, speed up of shelf inflowing glaciers are recorded, highlighting the buttressing force ice shelves have on glaciers and therefore their importance to the Antarctic Ice Sheet.

Many factors play a role in the integrity of ice shelves, besides melting on the surface and from below, for example waves, tides and pressure imbalances at the calving front (buoyancy force) induce bending moments on short time scales.

While ice is mainly studied as purely viscous, elastic or even plastic material. Certain ice shelf behaviour, like the tidal deformation near the grounding line or drainage of melt ponds, can only be sufficiently simulated with an viscoelastic model -representing the viscous thinning of ice as well as short term elastic bending deformations. Viscoelastic studies of ice so far have used commercial proprietary software, which require expensive licences and often lack inside in the used algorithms or flexibility to adapt the code as needed. In this thesis a viscoelastic model was developed in the open source FEM C++ environment Rheolef. First and foremost I would like to express my thank and gratitude to my supervisors for their support. Especially to Steph and Fabien for their patience, gentle pushes and continues inspiration throughout this thesis, without whom I would have given up. To Manu for his guidance

R ésum é

Le changement climatique a engendré des pertes considérables pour les glaces continentales au cours de ces dernières décennies et perdure toujours aujourd'hui.

L'Antarctique, la plus grande réserve de glace naturelle recèle à lui seul l'équivalent d'une hausse potentielle du niveau marin mondial de plusieurs dizaines de m, laquelle résulte non seulement de la fonte mais aussi de l'écoulement de la glace au travers les plateformes flottantes des glaciers émissaires, les 'ice shelves'. Depuis les dernières décennies caractérisées par des saisons de fonte plus longues et un océan de plus en plus chaud les 'ice shelves' s'amincissent et se désintégrent plus facilement. Parallèlement, une accélération des glaciers nourrissant ces shelves s'observe, mettant en évidence leur caractère d'arc-boutant sur l'écoulement de la glace vers la mer et leur importance sur la calotte glaciaire Antarctique.

De nombreux facteurs participent à la perte de cohésion des ice shelves. En plus de la fonte en surface et à la face inférieure, les vagues, les marrées ainsi que les différences de pression au niveau du front de vêlage induisent des moments de torsion à haute fréquence.

La glace est généralement modélisée comme étant soit purement visqueuse, soit élastique voire même plastique, cependant, certains des comportements du shelf comme la déformation liée aux marrées près de la ligne d'échouage ou la purge des lacs de surface ne peuvent être correctement simulés qu'avec un modèle viscoélastique reproduisant à la fois l'amincissement visqueux de la glace et le flambage élastique à court terme de la plaque.

Jusqu'à présent, les modélisations visco-élastiques de la glace ont fait appel

Motivation

The ice shelves are part of the cryosphere, which is one out of five Earth's climate system components -the others being atmosphere, lithosphere, biosphere and hydrosphere. These five spheres interact and influence each other in various ways, and it is a certain balance between them that ensures human existence on earth.

The cryosphere (Fig. 1.1) in general describes frozen elements on Earth, namely ice, permafrost and seasonal snow, and can be found around the poles as well as in high elevation areas like the Himalayan mountains. While permafrost stands for permanently frozen ground, and snow describes the frozen water vapour precipitation as well as the accumulation of fallen snow on the ground; ice encompasses all types of frozen water and can be divided into two groups -'sea ice' and 'land ice' -based on its formation. Lake, river and sea ice develop by phase change due to low temperatures; whereas ice sheets, ice caps and glaciers develop by densified accumulations of snow, i.e. meteoric ice. The latter, so-called 'land ice', spreads out or creeps down mountain slopes due to gravity, a process also described as ice flow. Once the ice reaches water, it starts to float because of its lower density compared to that of water (densities: glacier ice: ∼ 910 kg m -3 , freshwater (0°C):

1 000 kg m -3 , seawater: ∼ 1 028 kg m -3 [START_REF] Scambos | Ice shelf disintegration by plate bending and hydro-fracture: Satellite observations and model results of the 2008 Wilkins ice shelf break-ups[END_REF][START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF][START_REF] Cuffey | The Physics of Glaciers[END_REF][START_REF] Vaughan | Observations: Cryosphere[END_REF]). These floating tongues of marine terminating glaciers are called ice shelves. Land ice covers approximately 10 % of Earth's land surface and consists of frozen freshwater equal to ∼ 66 m of sea water level rise [Cuffey and Paterson 2010, p. 1], [Vaughan et al. 2013, p. 321]. Since ice shelves are floating, their contribution to sea level is already accounted for due to the displacement of water they induce [START_REF] Church | Sea Level Change[END_REF][Church et al. , p. 1172]]. But observations have shown that ice shelves hold back the glaciers they are fed by -a phenomenon often described as buttressing. One of the most thoroughly documented ice shelf break ups is the 2002 Larsen B ice shelf disintegration on the Antarctic Peninsula. Studies by, e.g. , [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF] and [START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF] showed an up to 8-fold increase in ice flow velocities with a concomitant thinning of shelf-inflowing glaciers compared to pre-break up years; while nearby glaciers uneffected by the break-up event did not demonstrate any comparable behaviour. However it does not take an extreme event like the complete disintegration of an ice shelf; the calving of a large iceberg can already remove enough buttressing [START_REF] Fürst | The safety band of Antarctic ice shelves[END_REF] to increase the velocity of an ice shelf inflowing glacier and hence cause thinning of it.

Glacier or ice shelf thinning induced by either buttressing removal or melt on their surface or below by warm waters can set the grounded glacier front afloat and thereby cause the grounding line to migrate inland. A grounded glacier describes ice that sits on bedrock which is below sea water level. This occurs if the ice is thicker than the bedrock's submergences below sea water level, under consideration of the density ratio between ice and water. The grounding line is generally referred to the location where grounded ice starts to float (bold vertical dashed line in Fig. 1.2). In the specific case of the glacier sitting on an inland sloped ('retrograde') bedrock, as sketched in Fig. 1.2, the glacier conditions are unstable, which can lead to continuous retreat of the grounding line until the bedrock slopes upward again. This processes of grounding line retreat is call Marine Ice Sheet Instability (MISI) and can cause a rapid reduction of the grounded ice sheet if the retrograde sloped area is very large, unless the glacier is stabilized by a sufficiently large ice shelf that exerts enough buttressing to reduce the ice flux [START_REF] Goldberg | Grounding line movement and ice shelf buttressing in marine ice sheets[END_REF][START_REF] Gudmundsson | Ice-shelf buttressing and the stability of marine ice sheets[END_REF]. Already in the 1970's the question of the stability of West Antarctic Ice Sheet was raised [START_REF] Hughes | Is the west Antarctic Ice Sheet disintegrating[END_REF][START_REF] Weertman | Stability of the Junction of an Ice Sheet and an Ice Shelf[END_REF]] and in the last few years it has emerged that specifically the West Antarctic Ice Sheet is prone to MISI [DeConto and Pollard 2016;[START_REF] Pattyn | The paradigm shift in Antarctic ice sheet modelling[END_REF]. Therefore it can be concluded that the disintegration or reduction of ice shelves affects inflowing glaciers and hence indirectly contributes to sea level rise [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF]].

Ice does not only store a lot of water in solid form on land, it also reflects large amounts of incoming short wave radiation back into space. While ice and snow have a reflectivity of 60 -90 %, depending on their surface cleanness and internal structure; open waters, i.e. the oceans, absorb ∼ 93 % of the incoming energy [Lubin and Massom 2006, p. 324, 351]. Therefore, the reduction of the ice surface on earth does not only reduce the area that can reflect energy but also, especially in case of floating ice, increases the area that can absorb and store heat.

The raising ocean temperatures cause a gain in average ocean volume due to thermal expansion, which has been stated as one of the main contributors to global sea level rise in the last decades by IPCCs 5 th Assessment Report in 2013 [START_REF] Church | Sea Level Change[END_REF][Church et al. , p. 1139[Church et al. , 1143]]. A warming ocean also increases the melting of the floating ice from below. The released freshwater, i.e. meltwater, to the oceans decreases its salinity and density; and together with the warming temperatures can slow down the thermohaline circulation. The thermohaline circulation -also known as the global oceanic conveyor belt -is an important part of the overall oceanic current system, that exchanges warm surface water with cold deep water due to horizontal differences in temperature and salinity. A slow down of the conveyor belt or even complete shutdown, as some scientists controversially say, will have a large impact on the overall ocean currents and further to our climate [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF][START_REF] Rafferty | Thermohaline circulation[END_REF][START_REF] Gordon | Ocean current[END_REF].

Approximately 90 % of the ice on earth is located in Antarctica, containing of about 70 % of the world's freshwater [Massom and Lubin 2006, p. 4]. The IPCCs 5 th Assessment Report [Vaughan et al. 2013, p. 320] highlighted that the Antarctic ice sheet has been losing mass during the last two decades, mainly due to the acceleration of outlet glaciers. They also found that two out of three identified processes linked to dynamical changes of the Antarctic ice sheet due to Climate Change are related to the stability of ice shelves (destabilisation through melt ponds and sub-shelf melting) [START_REF] Church | Sea Level Change[END_REF][Church et al. , p. 1172]]. While the amount of ice in ice shelves is relatively small compared to the grounded ice, but considering that most of Antarctica is surrounded by ice shelves [Cuffey and Paterson 2010, p. 373] together with the impact they have on the stability of the grounded ice (buttressing) makes them very important.

With the increase of computational power over the last decades, modelling of ice flow evolved from very simplified analytical models [START_REF] Nye | The flow of glaciers and ice-sheets as a problem in plasticity[END_REF][START_REF] Bödvarsson | On the Flow of Ice-Sheets and Glaciers[END_REF],

to complex 3-dimensional non-linear non-Newtonian approaches in finite elements, e.g. [START_REF] Rosier | The Interaction Between Tides, Ice Shelves and Ice Streams[END_REF]. Modelling of a process like ice flow is used to simulate the behaviour of a material to predict its unknown future or past state, but also to verify our understanding of the processes themselves, by simulating a period in the past of which the final state is known. With most naturally occurring ice being in remote locations and difficult to access, modelling gives us a good and non destructive way to verify our knowledge of its behaviour as well as help us to identify any lack thereof.

Most ice flow modelling approaches in the past have focused on purely viscous flows as it describes ice sheet dynamics over very long time scales, i.e. years, well enough. Unfortunately, these viscous models can only simulate the deformation causing thinning of the ice but can not replicate bending, e.g. deformations due to tidal movement or melt pond drainage, as do elastic models. On some occasions elastic models are used to replicate bending deformation [START_REF] Vaughan | Tidal flexure at ice shelf margins[END_REF][START_REF] Rignot | Tidal motion, ice velocity and melt rate of Petermann Gletscher, Greenland, measured from radar interferometry[END_REF]], yet these neglect the viscous thinning and are therefore only accurate on very short time scale.

In the last two decades, since the Larsen B ice shelf break up, studying ice shelves with a viscoelastic model has become increasingly popular. Starting with the more complex 4-element viscoelastic Burgers model by [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF] and a paper by [START_REF] Gudmundsson | Ice-stream response to ocean tides and the form of the basal sliding law[END_REF] indicating that the simpler 2-element Maxwell model might be sufficient enough. These were followed by studies on melt water ponds on ice shelves [START_REF] Macayeal | A model of viscoelastic ice-shelf flexure[END_REF], tidal forcing impacts, hence grounding line migration [START_REF] Rosier | The Interaction Between Tides, Ice Shelves and Ice Streams[END_REF]], calving processes [START_REF] Christmann | Viscoelastic Modeling of Calving Processes at Antarctic Ice Shelves[END_REF]] and local grounding of ice shelves on so called ice rises [START_REF] Macayeal | Treatment of ice-shelf evolution combining flow and flexure[END_REF]]. Unfortunately, most complex fluid models use proprietary software, e.g. COMSOL [MacAyeal et al. 2015[START_REF] Luckman | McDonald ice rumples on Brunt ice shelf Antarctica[END_REF],

MATLAB [START_REF] Christmann | Viscoelastic Modeling of Calving Processes at Antarctic Ice Shelves[END_REF]] and MSC.Marc [START_REF] Gudmundsson | Ice-stream response to ocean tides and the form of the basal sliding law[END_REF][START_REF] Rosier | The Interaction Between Tides, Ice Shelves and Ice Streams[END_REF],

where the usage of these models is restricted to licence availability and often lack transparency of what is inside the algorithms used or lack flexibility to adapt the code as needed.

This thesis aim was to develop a viscoelastic model for the ice shelf problem using the open source software Rheolef. This program had successfully been used

by [START_REF] Dansereau | A Maxwell-Elasto-Brittle model for the drift and deformation of sea ice[END_REF] for a viscoplastic model for sea ice.

Chapter 2 Background

Antarctica and its ice

Antarctica, the 5 th largest continent on Earth covers an area of ∼ 12.3 × 10 6 km 2

excluding ice shelves and islands [BAS 2005], this is about half the size of Africa. It is located far south at and surrounding the south pole and it has no land connection above sea water level to any other continent. Antarctica is generally split into 3 parts: the Antarctic Peninsula, West Antarctica and East Antarctica (see Fig. terface is largely below sea water level [START_REF] Morlighem | Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet[END_REF]], which makes it prone to sea water undercutting.

Since the AIS is spreading out towards the sea, the ice eventually starts to float. This way Antarctica has accumulated a large number of ice shelves (shown as grey areas in Fig. 2.1). Approximately 74 % of Antarctica is fringed with ice shelves [START_REF] Hogg | Impacts of the Larsen-C Ice Shelf calving event[END_REF], the largest being the Ross Ice Shelf and the Ronne-Filchner Ice Shelf with an area of 472 900 km 2 and 422 420 km 2 respectively [START_REF] Scambos | MODIS-based Mosaic of Antarctica (MOA) data sets: Continent-wide surface morphology and snow grain size[END_REF]]. In the past decades ice shelves have significantly reduced in size, especially around the Antarctic Peninsula due to their 'low' latitude extent, i.e. distance to the pole [START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF]]. In early 1995, Larsen A Ice Shelf (∼ 1600 km 2 (Fig. 2.2 -left); formerly known as the northern Larsen Ice Shelf)

and the ice shelf in Prince Gustav Channel collapsed [START_REF] Rott | Rapid Collapse of Northern Larsen Ice Shelf, Antarctica[END_REF][START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF]]. In 1998 Wilkins Ice Shelf suffered a large retreat (∼ 1100 km 2 ) [START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF]] and in 2002 Larsen B Ice Shelf collapsed within days (∼ 3200 km 2 ) (Fig. 2.2 -right) [START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF]]. Wilkins Ice Shelf retreated again significantly (∼ 1800 km 2 ) in 2008 [START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF]] and in 2017 Larsen C Ice Shelf had a large iceberg broken off (∼ 6000 km 2 ) [START_REF] Hogg | Impacts of the Larsen-C Ice Shelf calving event[END_REF]]. This retreat is mainly attributed to warming air and ocean temperatures and longer melt periods [START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF]].

These recent ice shelf collapses highlighted the importance they have on the grounded ice sheet. Ice shelves represent a large mass that is floating on water and is generally enclosed in bays which causes lateral drag on their sides. Furthermore, ice shelves are sometimes grounded on small islands, so-called pining points, that cause additional drag. This drag or back force restricts the ice flow from the inland and consequently slows down the ice shelf inflowing glaciers (buttressing).

In the event of an ice shelf collapse, it was documented that formerly ice shelf inflowing glaciers sped up [START_REF] Rignot | Accelerated ice discharge from the Antarctic Peninsula following the collapse of Larsen B ice shelf[END_REF][START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF]], which caused a larger flux from grounded to floating ice and hence a larger increase in sea level. The acceleration of ice flow velocity of those glaciers also caused a thinning of these glaciers reaching far inland, hence contributing to mass loss of the AIS.

Following the realisation of the importance of ice shelf buttressing to the AIS, a study performed in 2016 by Fürst et al. investigated how much each ice shelf contributes to buttressing its inflowing glaciers and how much it does not, and hence is passive and whose loss will not have an effect on those glaciers. They found (Fig. 2.3) that only around 13 % of ice in all ice shelves is passive and that the ice shelves in the Bellingshausen Sea has the least passive ice.

Factors influencing Ice Shelves

Ice shelves depend highly on their surrounding environment, e.g. form of the embayment they are located in and their underlying bathymetry, inflow from grounded ice, water and air temperatures, atmospheric pressure, and many more [START_REF] Cuffey | The Physics of Glaciers[END_REF]. To understand and to model an ice shelf's behaviour, these environmental factors need to be considered. While some of those factors can be included as a variable in an equation, e.g. temperature, others have their own underlying process, e.g. tides.

In this section, I present an overview of physical processes and sources that have been identified in literature to particularly influence the flow of an ice-shelf as well as the shelf's stress state. Processes indicated with an asterisk (*) cause deformations on short time-scales and require an extended modelling approach, which will be discussed in the following chapter. This theory, first introduced in the 1960's and 70's [START_REF] Weertman | Can a water-filled crevasse reach the bottom surface of a glacier[END_REF][START_REF] Robin | Depth of water-fitted crevasses that are closely spaced[END_REF], explains purely the crack advancing through water-induced stress but not the creation of a crack, nor does it include the water source and its way into the crack, or any relation to seasons [START_REF] Benn | Calving processes and the dynamics of calving glaciers[END_REF][START_REF] Scambos | Ice shelf disintegration by plate bending and hydro-fracture: Satellite observations and model results of the 2008 Wilkins ice shelf break-ups[END_REF][START_REF] Nick | A physically based calving model applied to marine outlet glaciers and implications for the glacier dynamics[END_REF]]. water on ice has a lower albedo than ice, hence it absorbs more energy. This causes a stronger melt at the water-ice interface [START_REF] Banwell | Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes[END_REF][START_REF] Macayeal | The flexural dynamics of melting ice shelves[END_REF], which increases the pond covered area and depth, as well as the melt water content. Further, the accumulated melt water acts as an extra load on top of the ice shelf, causing the underlying shelf to bend down locally due to gravity. This increases the depression and allows for more water to be collected [START_REF] Banwell | Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes[END_REF][START_REF] Macayeal | The flexural dynamics of melting ice shelves[END_REF]. The local bending of the ice shelf results in additional stresses in the ice. Particularly, tensile stresses on the shelf surface within a limited distance from the lake's edge and on the shelf base at the center of the lake. These tensile stresses can cause fracturing and fracture propagation [START_REF] Banwell | Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes[END_REF][START_REF] Macayeal | The flexural dynamics of melting ice shelves[END_REF][START_REF] Macayeal | A model of viscoelastic ice-shelf flexure[END_REF]]. Hydro-fracturing at the base of the lake or fracture propagation from the shelf base can cause a sudden drainage of the pond [START_REF] Banwell | Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes[END_REF][START_REF] Macayeal | The flexural dynamics of melting ice shelves[END_REF]. The removal of the water (i.e. removal of the additional load)

Melt water ponds

on the shelf's surface causes hydrostatic rebound, inducing strong tensile stresses in the horizontal plane of the shelf, which can initiate crevassing [START_REF] Banwell | Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes[END_REF][START_REF] Macayeal | The flexural dynamics of melting ice shelves[END_REF][START_REF] Macayeal | A model of viscoelastic ice-shelf flexure[END_REF].

MacAyeal and Sergienko [2013] noted that surface lakes often occur in complex arrays or patterns and are usually less than a kilometre wide, but occasionally can reach a diameter of several kilometres. Further observations showed that the filling of the lakes takes ∼ 60 days or more (based on observations at Larsen B and Wilkins) [START_REF] Macayeal | The flexural dynamics of melting ice shelves[END_REF][START_REF] Scambos | The link between climate warming and break-up of ice shelves in the Antarctic Peninsula[END_REF]] and the draining due to completely penetrated ice fractures is assumed to be ∼ 6 hours (based on observations on Greenland) [START_REF] Macayeal | A model of viscoelastic ice-shelf flexure[END_REF]. A study by [START_REF] Banwell | Direct measurements of ice-shelf flexure caused by surface meltwater ponding and drainage[END_REF] on the McMurdo ice shelf showed that not all melt ponds drain rapidly through fractures but also slowly (within a few days) by overflow. Satellite observations of major break-up events on the Larsen B ( 2002) and Wilkins (2008/09) ice shelf have suggested that the sudden rapid drainage of an array of melt ponds can trigger an ice shelf to disintegrate. [START_REF] Scambos | Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment, Antarctica[END_REF] stated that on the Larsen B ice shelf "areas with melt ponds disintegrated; adjacent areas with few or none remained".

Although [START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF] concluded that melt ponds did not play a role in the Wilkins ice shelf break-up, as some areas with melt ponds were uneffected.

Drained melt ponds result in so-called dolines (dry lake bed depressions), due to incomplete rebound of the ice [START_REF] Macayeal | A model of viscoelastic ice-shelf flexure[END_REF]], as well as potential new crevasses, especially in case of melt pond arrays.

Buoyancy forces (*)

The buoyancy force is generally defined as the upward force acting opposed to gravity on an object partially or fully submerged in a fluid. The buoyancy force is in equilibrium if the hydrostatic pressure (of the fluid -here water) is equal to the lithostatic pressure (of the object -here ice) on the bottom of the object. The pressure is generally defined by density Ö gravitational acceleration Ö height or depth.

Since glacier ice has a lower density than seawater, the ice shelf is (partly immersed) floating on water. The difference in density also causes a pressure difference depending on depth on the ice shelf front -the lithostatic pressure of ice increases slower than the hydrostatic pressure of water (see Fig. and, with time, produces a damped oscillating profile with a wavelength of about double the ice thickness [START_REF] Scambos | Ice shelf disintegration by plate bending and hydro-fracture: Satellite observations and model results of the 2008 Wilkins ice shelf break-ups[END_REF][START_REF] Reeh | On The Calving of Ice From Floating Glaciers and Ice Shelves[END_REF]]. The tensile stresses can become large enough to produce small cracks in ice shelf surface, parallel to the cliff. In the absence of surface water, the fractures and ice flow will theoretically reduce the bending stress. In case of the existence of melt water, hydro-fracturing in the new cracks can cause rifting parallel to the ice shelf front. [START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF] concluded that buoyancy forces (in combination with basal melting and coalescing flow units) were responsible for tensile mode rift formations before Wilkins IS break-up in February 2008. or far away storms in the deep ocean [START_REF] Bromirski | Tsunami and infragravity waves impacting Antarctic ice shelves[END_REF]]. Bromirski et al.

[2010] linked all break-up events on Wilkin ice shelf 2008 with IG-wave energy arrival from the coast of Patagonia. Very long period gravity waves are waves with a period larger than 300 s (> 0.003 Hz) [START_REF] Bromirski | Tsunami and infragravity waves impacting Antarctic ice shelves[END_REF].

Waves in general are assumed to induce vibrations or episodic displacements that can fatigue and fracture ice at weak locations, e.g. expand existing fractures [START_REF] Macayeal | Transoceanic wave propagation links iceberg calving margins of Antarctica with storms in tropics and Northern Hemisphere[END_REF][START_REF] Bromirski | Tsunami and infragravity waves impacting Antarctic ice shelves[END_REF]. IG and VLP waves are considered a potential explanation for winter break-up [START_REF] Bromirski | Transoceanic infragravity waves impacting Antarctic ice shelves[END_REF], since the absence of melt water excludes hydro-fracturing as an explanation. Swells on the other hand have been shown to be damped by sea ice and hence were considered to be too weak. While swells only affect the ice shelf near its front, due to their short period and therefore short wavelength [START_REF] Bromirski | Transoceanic infragravity waves impacting Antarctic ice shelves[END_REF], VLP waves can reportedly propagate to the grounding zone [START_REF] Bromirski | Tsunami and infragravity waves impacting Antarctic ice shelves[END_REF]]. Floating ice shelves are naturally affected by the tides, where the periodic sea level change (Fig. 2.5 -top panel) forces an up-down movement on the shelves near the calving front as well as around the grounding line (G in Fig. 2.6). This vertical movement is causing the shelf's grounding line to migrate horizontally within an area dependent on the bedrock slope [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF][START_REF] Rosier | Insights into ice stream dynamics through modelling their response to tidal forcing[END_REF][START_REF] Le Meur | Two independent methods for mapping the grounding line of an outlet glacier -An example from the Astrolabe Glacier, Terre Adélie, Antarctica[END_REF]], hence the grounding line is actually a grounding zone and the vertical movement is used to determine its position. But the displacement induced by the tides is not limited to the floating shelf alone, stress propagation causes a bending of the ice in a damped oscillation, with small vertical displacements several kilometres upstream of the grounding line (F in Fig. 2.6) [START_REF] Fricker | Ice shelf grounding zone structure from ICESat laser altimetry[END_REF][START_REF] Rignot | Antarctic grounding line mapping from differential satellite radar interferometry[END_REF][START_REF] Walker | Ice-shelf tidal flexure and subglacial pressure variations[END_REF]]. This makes it very difficult to locate the actual grounding line (/zone). Studies have shown that its neighbouring points can be detemined by either laser altimetry (ICESat) [START_REF] Fricker | Ice shelf grounding zone structure from ICESat laser altimetry[END_REF]], SAR interferometry (DInSAR) [START_REF] Rignot | Antarctic grounding line mapping from differential satellite radar interferometry[END_REF] Tides, specifically fortnightly ocean tides, appear to also affect ice flow velocities several kilometres upstream of the grounding zone (Fig. 2.5 -bottom panel) [START_REF] Gudmundsson | Fortnightly variations in the flow velocity of Rutford Ice Stream, West Antarctica[END_REF][START_REF] Rosier | Insights into ice stream dynamics through modelling their response to tidal forcing[END_REF][START_REF] Robel | Tidal modulation of ice shelf buttressing stresses[END_REF]] and can potentially influence subglacial hydrology [START_REF] Rosier | Insights into ice stream dynamics through modelling their response to tidal forcing[END_REF]].

Tides

Subglacial melt

Subglacial melt describes the process of ice shelf base melting and hence thinning of the shelf from below. It depends on the ocean circulation below the shelf also known as the "ice pump" and on the bathymetry. The "ice pump" is a basic overturning circulation based on melt water having a lower density than seawater, and buoyancy forces transporting the melt water along the upward facing slope of the shelf to the surface. The outward movement of water near the shelf causes an inflow of water from below (see Fig.

Ice stream

Ice streams are corridors of fast flowing ice within a glacier or ice self. The ice flow within an ice stream is orders of magnitudes faster compared to its adjacent ice [START_REF] Bennett | Ice streams as the arteries of an ice sheet: their mechanics, stability and significance[END_REF]] and this can create cracks parallel to the flow direction in the shear margin between the fast and slow flowing ice [START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF]]. On ice shelves, fast flowing ice units are mostly extensions of inflowing glaciers or grounded ice streams. [START_REF] Glasser | A structural glaciological analysis of the 2002 Larsen B ice-shelf collapse[END_REF] found that rifts between fast and slow flowing ice preconditioned and seemed to set the extent of the Larsen B ice shelf disintegration in 2002.

2.2.8 Ice rises and ice rumples -also known as: pinning points Ice rises and ice rumples are a local grounding on the floating ice to the sea floor (Fig. 2.8) and can occur anywhere under the ice shelf [START_REF] Doake | Ice-shelf Stability[END_REF][START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF][START_REF] Matsuoka | Antarctic ice rises and rumples: Their properties and significance for ice-sheet dynamics and evolution[END_REF]. While ice rises are considered tiny ice caps with ice flowing down either side, ice rumples are lower with the ice shelf flowing over it.

The ice flows slower when grounded due to the friction between ice and bedrock.

Hence, ice piles up upstream of an ice rumple (see Fig. 2.9), which thickens the shelf in this area and can potentially stabilise the ice shelf temporarily [START_REF] Hughes | On the Disintegration of Ice Shelves: The Role of Fracture[END_REF][START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF][START_REF] Doake | Ice-shelf Stability[END_REF]]. Ice rises have been noted to divide an ice flow, causing the flow to slow down and creating a fracture zone in its wake [START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF]]. Rifts originated in those fracture zones have been identified to be able to propagate upstream of the ice rise, resulting in a destabilisation of the shelf in this area [START_REF] Doake | Ice-shelf Stability[END_REF][START_REF] Braun | Changes of Wilkins Ice Shelf over the past 15 years and inferences on its stability[END_REF]].

Ice flow

Ice and its movement has been of interest for a long time. In the 18th century, it was first suggested that gravity is the source of movement, and sliding as well as viscous deformation were proposed as form of movement. With the increase of means of travel -first expeditions to Greenland and Antarctica in the beginning of the 20th century -as well as observation techniques -from rods in the ice and photogrammetry to remote sensing -our knowledge of ice and its movement has increased [START_REF] Cuffey | The Physics of Glaciers[END_REF].

Common observation techniques currently are remote sensing, especially by satellite or airborne, or GPS and can measure ice flow velocities, ice extent, elevation and mass. These methods include analysing satellite images for crevasses and extent (e.g. calving front as shown in Fig. 2.2) and tracking visible features on the ice over time for speed measurements [START_REF] Luckman | ERS SAR feature-tracking measurement of outlet glacier velocities on a regional scale in East Greenland[END_REF]]. Laser altimetry is used for elevation measurements [START_REF] Fricker | Ice shelf grounding zone structure from ICESat laser altimetry[END_REF], interferometry (DInSAR or InSAR) for observations on surface height or speed changes [START_REF] Rignot | Tidal motion, ice velocity and melt rate of Petermann Gletscher, Greenland, measured from radar interferometry[END_REF][START_REF] Rignot | Antarctic grounding line mapping from differential satellite radar interferometry[END_REF]] and ice mass changes can be retrieved from the GRACE mission [START_REF] Wouters | Global Glacier Mass Loss During the GRACE Satellite Mission (2002-2016)[END_REF]. Ground-based (kinetic) GPS stations can record more precise position and elevation changes in shorter time intervals (compare to satellite measurements) [START_REF] Vaughan | Tidal flexure at ice shelf margins[END_REF][START_REF] Le Meur | Two independent methods for mapping the grounding line of an outlet glacier -An example from the Astrolabe Glacier, Terre Adélie, Antarctica[END_REF]. Especially satellite remote sensing, due to the remoteness of most ice on Earth, has the advantage of providing regular large scale data and hence allowing for a periodic monitoring of ice.

While the above mentioned techniques give a snapshot of the current state of a glacier or ice shelf, numerical models are needed in order to simulate the flow of the ice and hence predict its behaviour and extent. Simulations require knowledge of the nature of the flow as well as material specific properties.

The following subsections give an introduction of fluid classification from a rheology point of view, followed by ice rheology and common model approaches.

Flow Modelling

The study of the flow of a material, also known as rheology, is part of the fluid mechanics field, where the motion, hence deformation of a material is generally expressed through stress and strain related to time.

Overall the motion of a fluid is mathematically defined by a set of field and The dashpot on the other hand symbolises the viscous part, and shows an increasing displacement over time with an applied force, but no recovery on force removal. The dashpot represents the viscosity η -the resistance to deformation of a fluid, and presents a purely viscous flow in a single system, also known as the Newton model (Fig. 2.10.b).

There also exists a plastic component, which describes the non-recoverable elastic deformation. The point on which deformation is not recoverable is defined by the yield stress, but this is not presented due to not being relevant to this work. As the name already suggests, the mathematical definition -still based on Navier-Stokes as general flow equations -is more complex and relies on simplifications of the problem as well as time and/or space approximations [START_REF] Saramito | Complex fluids -Modeling and Algorithms[END_REF]].

Ice rheology

Laboratory and in situ experiments have shown that ice is not purely elastic, viscous or plastic, but shows a distinct creep and elastic behaviour (see Fig. 2.11) [START_REF] Glen | Experiments on the Deformation of Ice[END_REF][START_REF] Nye | The flow law of ice from measurements in glacier tun-nels, laboratory experiments and the Jungfraufirn borehole experiment[END_REF][START_REF] Jellinek | Viscoelastic Properties of Ice[END_REF][START_REF] Gold | Engineering Properties of Fresh-Water Ice[END_REF][START_REF] Smith | Viscous relations for the steady creep of polycrystalline ice[END_REF][START_REF] Budd | A review of ice rheology for ice sheet modelling[END_REF][START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF].

The creep behaviour -a time-dependent permanent deformation under load -of ice is described to be occurring in three main phases: a decelerating primary creep and a constant secondary and tertiary creep [START_REF] Smith | Viscous relations for the steady creep of polycrystalline ice[END_REF][START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF]. The secondary phase represents the lowest (minimum) shear rate, followed by an accelerating phase and a significantly larger tertiary creep. While both secondary and tertiary creep are suggested to be constant, the shear angle In 1955 Glen already formalised a power law viscosity function based on the minimum creep rate of ice, also known as the secondary creep, to better describe the flow of ice in simulations:

η = 1 2 A -1 n ε1-n n e (2.1)
where η is the effective viscosity, εe is the effective strain rate, A the creep pa- Glen's stress-strain-rate relationship, also known as Glen's flow law, is the most commonly used relation for the viscosity of ice [START_REF] Cuffey | The Physics of Glaciers[END_REF]. Other viscosity functions have been suggested over the years, e.g. , by [START_REF] Colbeck | A Flow Law for Temperate Glacier Ice[END_REF], [START_REF] Smith | Viscous relations for the steady creep of polycrystalline ice[END_REF] or [START_REF] Doake | Flow law for ice in polar ice sheets[END_REF], but these have never been widely used.

The instant elastic deformation shown Fig.

λ = η E (2.2)
where E is the elastic modulus, which is equal to the Young modulus in a 1Dsystem.

MacAyeal and Sergienko [2013] have described the relaxation time for ice to be around ∼ 10 days, but this is dependent on the material parameters used as input (see next subsection). [START_REF] Rignot | Tidal motion, ice velocity and melt rate of Petermann Gletscher, Greenland, measured from radar interferometry[END_REF][START_REF] Vaughan | Observations: Cryosphere[END_REF]] are often used. An ice sheet investigation with a thickness of ∼ 10 3 m magnitude would use a higher density than an ice shelf study with a thickness of ∼ 10 2 m magnitude.

Ice parameters

The viscosity (η) -resistance to deformation or flow -of ice depends largely on grain-scale structure, temperature and water content [Cuffey and Paterson 2010, p. 7] and ranges between 10 13 -10 17 Pa s for temperatures between 0 --20 °C, according to Greve and Blatter [2009, Fig. 4.6]. For linear viscoelastic models the following values have been found to fit observations best: 10 13.7 Pa s [START_REF] Wild | Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone[END_REF]], 2 -3 × 10 13 Pa s [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF]] and 10 15 Pa s [START_REF] Walker | Ice-shelf tidal flexure and subglacial pressure variations[END_REF].

The Young modulus (E) -instantaneous stiffness or the resistance to elastic deformation -is defined as stress over strain in a one-dimensional scenario and is only meaningful within the elastic range, which is usually defined by the Maxwell time (see subsection 2.3.2). According to [START_REF] Gold | Engineering Properties of Fresh-Water Ice[END_REF] the elastic modulus depends on ice temperature, density, type, purity, as well as stress frequency and load application direction. Common measured or retrieved values can be found in Tab. 2.1.

The Poisson's ratio (ν) is the ratio between relative contraction to relative expansion of an elastic material, with values mostly between 0 -0.5, where 0.5 stands for a elastic deformation of an incompressible material like rubber and 0 is close to cork. For ice values between 0.

Stokes approximations for ice

Studying ice flow with numerical models is very computationally expensive, because it is often necessary to be studied over a large area, requiring a large mesh as spatial geometry approximation, and over long time scales. Hence, there has always been the urge to reduce and simplify the flow equations to the absolute necessary components. Following, I introduce a few of the most common simplifications of the governing field equations.

Also, ice has been considered to be viscoelastic [START_REF] Jellinek | Viscoelastic Properties of Ice[END_REF][START_REF] Sinha | Short-Term Rheology of Polycrystalline Ice[END_REF][START_REF] Smith | Viscous relations for the steady creep of polycrystalline ice[END_REF] for some time, the elastic part is often neglected, because of its relevancy only at short time-scales (days) and specifically Paleostudies or future predictions often look at decade or century time-scales. Additionally, ice is a very slow flowing fluid, hence it has a very high viscosity (see subsection 2.3.3). This causes a very small (a/g << 1) inertia term in the Navier-Stokes equations and this leads to this term often being neglected. The Navier-Stokes equations without the inertia term is also known as the (Full-)Stokes equations [Greve and Blatter 2009, ch. 5].

The Stokes equations determine the full stress field, consisting of 9 components (or 6, due to the symmetry imposed by the conservation of angular momentum), at all time steps. Many flow models only approximate the stress field by neglecting non-relevant stress components, as demonstrated in Fig. 2.12. This often reduces the usability of a model to a specific scenario or area, but reduces the computational cost significantly. The Shallow Ice Approximation (SIA -Fig. 2.12.b) neglects stretching and compression in flow direction (longitudinal stresses) and lateral drag on the sides against valley walls or slower flowing ice (transverse stresses), as well as the stress gradient in vertical direction, because it assumes a small depth to width ratio.

This means SIA is based on the driving gravitational stress and the opposing basal drag [Morland 1984][Greve andBlatter 2009, ch. 5.4], which results in a lower velocity at the base of the ice compared to the top surface. This is demonstrated as black flow profile in the top image of Fig. The Shallow Shelf Approximation (SSA -Fig. 2.12.c), similarly to the SIA, assumes a small depth to width ratio as indicated by the word shallow in the name.

In contrast to the SIA, the SSA assumes no basal shear stress but a dominance of longitudinal stresses, and represents what is called a Plug-flow with a constant horizontal velocities over depth (green flow profile in the top image of Fig. 2.13) [Greve and Blatter 2009, ch. 6.3]. Altogether SSA is reduced to a 2D plane view problem with 2 degrees of freedom and is also a rather strong simplification of the Stokes equations, which makes it computationally cheap. It was designed specifically for ice shelves -where there is no basal drag, and is therefore not suitable for other areas. Further, SSA is a 2D depth-integrated model, hence it can not represent strong vertical variations, which makes it unsuitable at or near the grounding line.

SIA and SSA work for their designated areas over long time-scales, i.e. ice sheet and ice shelf respectively, but not for the 'transition zone' in between (see top image of Fig. 2.13), namely faster flowing ice streams and areas around the grounding line. Both approximations (SIA and SSA) have been coupled to a so-called Hybrid model for these areas. They are common in paleo-ice sheet modelling, because they are computationally cheap but also less accurate than Full-Stokes [START_REF] Kirchner | Shallow ice approximation, second order shallow ice approximation, and full Stokes models: A discussion of their roles in palaeo-ice sheet modelling and development[END_REF].

There are also higher-order models, like the Blatter-Pattyn model [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][START_REF] Dukowicz | Consistent approximations and boundary conditions for ice-sheet dynamics from a principle of least action[END_REF], which include more stress components and are therefore more accurate than the previously mentioned approximations while still being computationally cheaper than the Stokes equation. This makes them more suitable for 3-dimensional or more complex problems like the grounding line. The three graphs in Fig. 2.13 show the effect the amount of stress components considered have on the basal traction and hence give an idea of the accuracy of the result. The basal traction very much depends on the type of ground the ice sits on, e.g. rock or till, and the amount of water. It is expressed in so-called gliding, sliding or friction laws, which are not detailed in this work.

Ice shelf flow modelling approaches

A viscous power law regime, e.g. SSA or Stokes (subsection 2.3.4) with Glen's flow law (subsection 2.3.2), is deemed to be sufficient for ice flow for long-term studies (years or decades), where the viscous deformation is dominant (Fig. 2.11 on the right). But this regime cannot represent short-term responses like tidal movements (factors marked with * in section 2.2), which occur within hours or days and are therefore in the elastic dominant time range [START_REF] Jellinek | Viscoelastic Properties of Ice[END_REF][START_REF] Sinha | Short-Term Rheology of Polycrystalline Ice[END_REF][START_REF] Morland | Non-linear viscoelastic response of ice[END_REF][START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF][START_REF] Christmann | Viscous and viscoelastic stress states at the calving front of Antarctic ice shelves[END_REF]].

In fact, in the literature it can be quite often found that the bending moment induced to an ice shelf by the buoyancy force (subsection 2.2.3) [START_REF] Wagner | On the role of buoyant flexure in glacier calving[END_REF] or the tidal movement at the grounding line (subsection 2.2.5) [START_REF] Vaughan | Tidal flexure at ice shelf margins[END_REF][START_REF] Sayag | Elastic dynamics and tidal migration of grounding lines modify subglacial lubrication and melting[END_REF][START_REF] Le Meur | Two independent methods for mapping the grounding line of an outlet glacier -An example from the Astrolabe Glacier, Terre Adélie, Antarctica[END_REF]] is modelled with a purely elastic thin-beam or -plate model. But [START_REF] Wild | Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone[END_REF] suggests that the purely elastic model can only capture the impact of spring tides well, for lower tidal oscillations a viscoelastic model performs better compared to observed data than a purely elastic model.

Already [START_REF] Jellinek | Viscoelastic Properties of Ice[END_REF] and [START_REF] Sinha | Short-Term Rheology of Polycrystalline Ice[END_REF] suggested that ice would be best represented with a viscoelastic model, but it took until [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF] to really be investigated and used. [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF] started with a Burgers model (Fig. 2015] and bending of the ice shelf front due to buoyancy forces [START_REF] Christmann | Viscous and viscoelastic stress states at the calving front of Antarctic ice shelves[END_REF]].

It can be quite often found that for the viscoelastic models the thin-beam or -plate theory is adapted [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF][START_REF] Walker | Ice-shelf tidal flexure and subglacial pressure variations[END_REF][START_REF] Macayeal | A model of viscoelastic ice-shelf flexure[END_REF][START_REF] Christmann | Viscous and viscoelastic stress states at the calving front of Antarctic ice shelves[END_REF][START_REF] Wild | Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone[END_REF], which assumes a constant thickness of the shelf and hence simplifies the equation, and only a few use the Full Stokes approach [START_REF] Gudmundsson | Ice-stream response to ocean tides and the form of the basal sliding law[END_REF][START_REF] Rosier | Insights into ice stream dynamics through modelling their response to tidal forcing[END_REF][START_REF] Rosier | Temporal variations in the flow of a large Antarctic ice stream controlled by tidally induced changes in the subglacial water system[END_REF][START_REF] Rosier | On the interpretation of ice-shelf flexure measurements[END_REF].

Studies with a linear viscoelastic model gave the best fitting results to observations with a viscosity between 10 13 -10 15 Pa s [START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF][START_REF] Walker | Ice-shelf tidal flexure and subglacial pressure variations[END_REF][START_REF] Wild | Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone[END_REF]] (see subsection 2.3.3), while it was found that results with viscosities of 10 16 Pa s [START_REF] Wild | Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone[END_REF] or 10 17 Pa s [START_REF] Walker | Ice-shelf tidal flexure and subglacial pressure variations[END_REF]] resembled the purely elastic results. It should be noted that these studies were performed with the thin-plate theory assumption.

The Young modulus shows a much larger range of value in Tab. 2.1. Noticeable here, is the large difference between E values retrieved in lab experiments (∼ 10 GPa) and elastic models fitted to observed data (1 MPa -3 GPa), which has a value range of about 3 orders of magnitude itself. [START_REF] Sayag | Elastic dynamics and tidal migration of grounding lines modify subglacial lubrication and melting[END_REF] argued that this is due to the hard bed assumption in the elastic model, while they retrieve good fitting results with a elastic model and a Young modulus of 9.3 GPs, equal to the laboratory value, and a soft bed assumption. Further, [START_REF] Gudmundsson | Ice-stream response to ocean tides and the form of the basal sliding law[END_REF] noted that in situ data, due to the time it was measured over, contains a combination of Chapter 3

Model development

Mathematical and Physical Background

In this section the general stress tensor and stain rate tensor (subsection 3.1.1) will be presented first, followed by the governing equations for the motion of a fluid consisting of the field equations (subsection 3.1.2) expressing conservation laws and the constitutive equation (subsection 3.1.3) describing the relationship between stress and strain. This section is then concluded with the retrieval of the viscous and viscoelastic Stokes equation in their strong (subsection 3.1.4) and weak formulation (subsection 3.1.5).

Stress & Strain tensor

The state of stress is completely defined at a particular point when expressed as a second-order Cartesian tensor, also known as Cauchy stress tensor σ.

σ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ xx σ xy σ xz σ yx σ yy σ yz σ zx σ zy σ zz ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (3.1)
The Cauchy stress tensor consists of 9 components, the 3 normal stress components:

σ xx ,
σ m = tr(σ) = σ xx + σ yy + σ zz (3.2)
and is directly related to the isotropic pressure p, the negative average of all normal stresses:

p = - σ xx + σ yy + σ zz 3 = - 1 3 σ m (3.3)
The deviatoric stress components act parallel to a plane and relate to the shape change or deformation. The deviatoric stress tensor τ is defined as the total stress (σ) minus the average of all normal stresses (σ m ):

τ = dev(σ) = σ - 1 3 σ m I (3.4)
where I is is the identity matrix. Under consideration of eq. 3.3 the deviatoric stress than writes:

τ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ τ xx τ xy τ xz τ yx τ yy τ yz τ zx τ zy τ zz ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = σ + pI = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ xx + p σ xy σ xz σ yx σ yy + p σ yz σ zx σ zy σ zz + p ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (3.5)
and following eq. 3.5 the Cauchy stress tensor can also be expressed as deviatoric stress (τ ) minus the pressure (p):

σ = τ -pI (3.6)
The stress tensor and the deviatoric stress tensor are often considered symmetric σ = σ T and τ = τ T (σ xy = σ yx , σ yz = σ zy and σ xz = σ zx ), due to the conservation of angular momentum, which essentially is reducing the number of unknowns from 9 to 6.

The strain is approximated with the infinitesimal strain tensor ϵ, which holds for small shape changes and is defined as the symmetrical part of the displacement gradient. The strain tensor is often expressed in its time derivative form, the strain rate tenor ε, which consequently equals to the symmetrical part of the velocity gradient.

The velocity gradient grad(u) can be expressed as the sum of the symmetrical part D(u), representing the deformation rate, and the asymmetrical part W (u),

representing the vorticity or rotation rate:

grad(u) = D(u) + W (u) (3.7)
where u is the velocity vector, and the symmetrical and asymmetrical can be expressed as:

D(u) = grad(u) + grad(u) T 2 (3.8) W (u) = grad(u) -grad(u) T 2 (3.9)
Hence the strain rate writes:

ε = D(u) = grad(u) + grad(u) T 2 (3.10)
Similarly to the stress tensor, the strain and strain rate tensor can be expressed in their decomposed form with its principal (tr(ϵ) and tr( ε)) and deviatoric (dev(ϵ) = ε and dev( ε) = ε) parts:

ϵ = ε + 1 3 tr(ϵ)I (3.11) ε = ε + 1 3 tr( ε)I (3.12)

Field equations

The field equations in fluid mechanics describe the conservation of mass and momentum in a closed system.

The conservation of mass, also know as the continuity equation, states that in a closed system the total mass inflow into that system must be equal to the total mass outflow plus the accumulated mass in the system. This expressed in the differential form says the sum of gain or loss of mass in the system over time and the difference between inflow and outflow of the system is 0:

∂ρ ∂t + div(ρu) = 0 (3.13)
where ρ is the density of the fluid, t is the time and u is the velocity vector field.

In case of an incompressible material, an assumption often made for ice, the density of the fluid is constant (ρ = const.) and eq. 3.13 simplifies to:

div(u) = 0 (3.14)
The conservation of linear momentum is based on Newton's second law of motion, stating that the sum of all forces is equal to the change of linear momentum over time, expressed as mass times velocity.

ρ ( ∂u ∂t + (u.∇)u) Lagrangian derivative + ∂ρ ∂t + div(ρu) mass = div(σ) surface forces + ρg volume forces (3.15)
where g is the gravitation and the Lagrangian derivative, the so-called material derivative or inertia, expresses the acceleration over time and space of a material.

Under consideration of the conservation of mass (eq. 3.13) the conservation of linear momentum (eq. 3.15) simplifies to:

ρ ( ∂u ∂t + (u.∇)u) = div(σ) + ρg (3.16)
In case of very slow moving fluids, like ice, the inertia term is so small that it can be neglected, hence the linear momentum conservation equation (eq. 3.16) simplifies further to: div(σ) + ρg = 0 (3.17)

Constitutive equation

The constitutive equation is material specific and describes the deformation behaviour based on the applied forces, hence the relationship between deviatoric stress (τ ) and strain (ε).

The simplest model to describe a viscoelastic fluid is the 2-element Maxwell model (Fig. 2.10.c), where both elements, one representing the elastic part and the other the viscous part, are connected in series. Following this configuration the total stress (σ total ) and strain (ϵ total ) can be expressed in dependens of their elastic and viscous part as following:

σ total = σ elastic = σ viscous (3.18) ϵ total = ϵ elastic + ϵ viscous (3.19)
The elastic part, with the spring (Fig. 2.10.a) as its mechanical representative, is described by Hooke's law. Hooke's law for a spring states that the extension is proportional to the force, and can be extended for a continuous elastic material to the strain is proportional to the stress by the material specific elastic modulus.

The linear relationship (E = const.) can be written as:

σ = 2Gϵ (3.20)
with the shear modulus G being equal to the Young modulus E in an one-dimensional system, but dependent on the Poisson's ratio ν for higher dimensions:

G = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ E for 1D case (3.21a) E 2(1 + ν) for 2-3D case (3.21b)
where eq. 3.21b is the simplified expression for isotropic materials [START_REF] Vincent | Basic Elasticity and Viscoelasticity[END_REF]].

Hooke's law in eq. 3.20 can be expressed for deviatoric strain (ε) using the decomposition of the strain tensor (eq. 3.11) and of the stress tensor (eq. 3.4 and

3.2): 2ε = τ G (3.22)
The viscous part, depicted by the dashpot (Fig. 2.10.b), is also called the Netwonian part. It describes the deformation over time based on the applied stress and the material specific viscosity η. The linear relationship (η = const.) can be expressed in relation to the deviatoric strain rate ( ε) as

2 ε = τ η (3.23)
Hence, the Maxwell elastic-viscous relationship (eq. 3.19) expressed as its derivative over time (strain rate) writes:

2 ε = τ η + ▿ τ G (3.24)
where ▿ τ is the upper-convected derivative of the deviatoric stress tensor, which is commonly used for slow flows like ice. This time derivative is used because the constitutive equation has to be frame invariant and the Langrangian derivative of a tensor itself is not frame invariant as shown by Saramito [2016, ch. 4.2]. A vector or tensor is frame invariant, or objective, if it is the same for an observer at rest or one that is animated by rigid motion. Here a term describing the fluids deformation and rotation is added to the Langrangian derivative to reach objectivity:

▿ τ = ∂τ ∂t + (u.∇)τ Lagrangian derivative -grad(u) τ -τ grad(u) T
objective term

(3.25)

Flow equations

By combining the field equations (subsection 3.1.2) and constitutive equation (subsection 3.1.3) we obtain a set of equations describing the flow of a material.

Because ice is a very slow moving fluid and in our case considered incompressible, the flow equations are derived for this simplified case.

For the viscous case, the Newtonian stress-strain relationship (eq. 3.23) needs to be expressed for the deviatoric stress:

τ = 2ηD(u) (3.26)
where the strain rate ε is written in dependence of the velocity u (eq. 3.10). The Newtonian deviatoric stress can then be inserted in the decomposed stress tensor (eq. 3.6):

σ = 2ηD(u) -pI (3.27)
which, when inserted for the stress tensor in the conservation of linear momentum (eq. 3.17), results in the linear viscous Full-Stokes equation:

-div(2ηD(u)) + div(pI)

grad(p) = ρg (3.28)
For the viscoelastic case, the Maxwell stress-strain relationship (eq. 3.24) is expressed for the shear stress:

τ = 2ηD(u) - η G ▿ τ (3.29)
where the strain rate ε is written in dependence of the velocity u (eq. 3.10) and inserted into the decomposed stress tensor (eq. 3.6):

σ = 2ηD(u) - η G ▿ τ -pI (3.30)
When the above equation is inserted into the conservation of linear momentum (eq. 3.17) it results in the linear viscoelastic Full-Stokes equation:

-div(2ηD(u)) + div ( η G ▿ τ ) + grad(p) = ρg (3.31)
The viscous or viscoelastic Full-Stokes equation (eq. 3.28 or eq. 3.31 respectively) together with the conservation equation for mass (eq. 3.14) are the strong formulation of the flow problem of a very slow moving and incompressible fluid.

Weak formulation of flow equations

In cases when solving a complex problem with e.g. the Finite Element Method (see next section 3.2), the so-called strong form of equations, as introduced in the previous section (subsection 3.1.4) need to be converted into their weak form. The weak expression of a problem relaxes the continuity requirement by introducing test functions, initial and boundary conditions, which makes it easier to solve.

This means the solution might not be exact but a good approximation that satisfies the equations. The conversion from the strong form to its weak or variational formulation is done by multiplying the equation by a test function of the same order (scalar, vector, tensor) as the equation and then is integrated over the model domain Ω.

The mass conservation (eq. 3.14) is a scalar equation, hence it is multiplied by a scalar test function q and integrated over the domain Ω:

∫ Ω q div(u) dx = 0 (3.32)
The linear momentum conservation (eq. 3.17) is a vector equation and is therefore re-written with a vector test function v:

-∫ Ω v.div(σ) dx = ∫ Ω v.(ρg) dx (3.33)
by using integration by part, the left hand side can be re-written as following:

-

∫ ∂Ω v.(σ.n) ds + ∫ Ω D(v) ∶ σ dx = ∫ Ω v.(ρg) dx (3.34)
where the first term is the integration of a stress that is acting normal on the boundary. Hence, this term defines a boundary condition and can be moved to the right hand side of the equation.

Now the constitutive equation (eq. 3.27 and eq. 3.30) can be inserted into the weak form of the momentum equation (eq. 3.34) to gain the weak formulation of the Full-Stokes:

viscous case:

∫ Ω D(v) ∶ (2ηD(u)) dx -∫ Ω D(v) ∶ (pI) dx = ∫ Ω v.(ρg) dx + ∫ ∂Ω v.(σ.n) ds (3.35)
viscoelastic case:

∫ Ω D(v) ∶ (2ηD(u)) dx -∫ Ω D(v) ∶ ( η E ▿ τ ) dx -∫ Ω D(v) ∶ (pI) dx = ∫ Ω v.(ρg) dx + ∫ ∂Ω v.(σ.n) ds (3.36)

Finite Element Method

For complex problems, e.g. time and space dependent problems, and/or complex geometries, an analytical solution is often not possible. A common method is the Finite Element Method (FEM) to find an approximated solution for such kind of problems.

FEM is a numerical method for solving partial differential equations (PDE) on smaller discrete sub-domains, so-called elements. This requires the continuous geometry (Ω in Fig. 3. Nodes are points with a distinct position in space, defined by coordinates, that have the field variables and sometimes their degrees of freedom defined on them.

The element can be in different shapes and sizes. The most common shape for a 2-dimensional element are a 3-node triangle (as shown in Fig. 3.2) or a 4-node quadrilateral.

The field variables or unknowns of a set of equations are calculated on the nodes of the element. The calculation of the field variables on the nodes is done through the weak or variational formulation of the PDE with a so-called shape or basis function often described as polynomials. The most common shape functions in fluid dynamics are the 1 st -order polynomial (P1) for a linear interpolation and a 2 nd -order polynomial (P2) for a quadratic interpolation within an element between nodes. For a problem with the velocity and pressure as their unknowns, a P2-P1 combination (P2 for the velocity, P1 for the pressure) also known as Taylor-Hoodpair, is the most used.

Further, initial and boundary conditions are needed to be defined. Boundary conditions (BC) form the connection between the modelling domain and its surrounding environment and are defined at the boundary of the domain where they act on. These boundary conditions are either set directly on the unknown field variable as fix value, as for example the inflow velocity of a fluid. This type of BC reduces the degrees of freedom at a node and are referred to as Dirichlet BC.

Forces, like the water pressure, that act normal to the boundary are described as surface forces (σ.n in eq. 3.34 -eq. 3.36)

Time dependent and independent problems

The time dependency of a problem depends primarily on the flow equations and secondarily on the problem definition, i.e. geometry, material properties and boundary conditions.

A problem is time independent if the time derivatives in the equations can be set to zero, e.g. the viscous Stokes equation (eq. 3.28) and if the problem definition, i.e. input parameters, are constant and do not change over time. This often simplifies the equation so that the solution can be directly computed, which leads to the so-called steady state solution.

In complex multi-physics problems, especially with non-linearities, the time derivative terms in the equation are important and cannot be set to zero, e.g. the upper-convected derivative term (eq. Transient problems, where one or more input parameters change with time, e.g. the sea water level due to tides, usually require an initial condition, hence the problem is started from a (near) steady state solution.

In the problems presented in this project, the viscous Stokes problem (eq. 3.28) has no time derivative and a steady state can be computed directly, while the viscoelastic Stokes problem (eq. 3.31) contains a time derivative in the upperconvected derivative term (eq. 3.25) and needs to be computed with the transient analysis approach.

For the transient analysis approach, the divergence rather than the convergence of the solution was tested. The model would run for a set amount of iterations or until a set total time was reached and only aborts the model run if the solution diverged. The divergence is evaluated by adding together the absolute maximum difference of the shear stress and velocity between the current and previous solution, and comparing it to the previous absolute maximum difference. If the current difference is 10 times larger than the previous difference, the solution is said to be diverging.

Rheolef

For this work I used the latest version of a software called Rheolef (v7.1). Rheolef example [Saramito 2020, ch. 2.1.4] as starting point for the viscous problem and his 'Oldroyd theta-scheme' example [Saramito 2020, ch. 5.3.3] for the viscoelastic problem.

The advantage of this open source program, compared to commercial software, is the flexibility and inside in the used methodology. The θ-scheme algorithm, used for the viscoelastic problem, is introduced in the next section 3.2.3. The disadvantage is the required knowledge in C++ as well as a deeper understanding of FEM.

θ-scheme algorithm

The θ-scheme algorithm is an operator splitting time discretisation algorithm For the ease of presentation, let U be the field variables, U 0 their initial value, M the coefficients of time derivatives, A the operator containing the flow equations -the constitutive equation eq. 3.29 (with eq. 3.25 and eq. 2.

2), the momentum conservation law eq. 3.17 (with eq. 3.6) and the mass conservation law eq. 3.14, and F the external force, containing the right hand side of eq. 3.17, be defines as:

U = ⎛ ⎜ ⎜ ⎜ ⎝ τ u p ⎞ ⎟ ⎟ ⎟ ⎠ , U 0 = ⎛ ⎜ ⎜ ⎜ ⎝ τ 0 u 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ , M = ⎛ ⎜ ⎜ ⎜ ⎝ λ 2η 0 0 0 -ρ 0 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ , A(U) = ⎛ ⎜ ⎜ ⎜ ⎝ τ -2ηD(u) + λ ((u.∇)τ -grad(u) τ -τ grad(u) T ) div(τ ) + grad(p) div(u) ⎞ ⎟ ⎟ ⎟ ⎠ and F = ⎛ ⎜ ⎜ ⎜ ⎝ 0 -ρg 0 ⎞ ⎟ ⎟ ⎟ ⎠
The problem, based on Newtons 2 nd law, then writes without boundary condition:

find U such that M ∂U ∂t + A(U) = F MU(t = 0) = MU 0
In the θ-scheme algorithm, the operator A is split into 2 convenient parts, namely the Full-Stokes (A 1 ) and the transport problem (A 2 ):

A(U) = ⎛ ⎜ ⎜ ⎜ ⎝ τ 2η -D(u) div(τ ) + grad(p) div(u) ⎞ ⎟ ⎟ ⎟ ⎠ A 1 (U) + ⎛ ⎜ ⎜ ⎜ ⎝ λ 2η ((u.∇)τ -grad(u) τ -τ grad(u) T ) 0 0 ⎞ ⎟ ⎟ ⎟ ⎠ A 2 (U)
These 2 sub-problems (A 1 and A 2 ) are used to compute the overall problem in 3 steps, for which the time step (∆t) is split into 3 sections by a factor θ, as shown in Fig. 3.3. The factor θ can take any numerical value between ]0, 1 2 [, but Saramito recommends:

θ = 1 - 1 √ 2 (3.37)
where the Full-Stokes problem (A 1 ) is solved in step 1 and 3 and the transport problem (A 2 ) in step 2. The 3 steps with the 2 sub-problems unpacked then write for m = 0, U 0 is given and for m ≥ 0, U m is known:

Step 1 (t .. t + θ∆t) find τ m+θ , u m+θ and p m+θ (1a) with τ m and u m known, solve the Full-Stokes problem for u m+θ and p m+θ :

ρ θ∆t u m+θ -div (2 ( θ∆t λ + θ∆t ) ηD(u m+θ )) + grad(p m+θ ) = f m (3.38) -div(u m+θ ) = 0 (3.39)
with

f m = ρg + ρ θ∆t u m + 1 λ + θ∆t div (λτ m -λθ∆tγ m ) (3.40) γ m = (u m .∇)τ m -grad(u m ) τ m -τ m grad(u m ) T (3.41) (1b
) with τ m and u m+θ known, τ m+θ can be compute explicitly from:

τ m+θ = 1 λ + θ∆t (λτ m + 2ηθ∆tD(u m+θ ) -λθ∆tγ m ) (3.42)
Step 2 (t + θ∆t .. t + (1 -θ)∆t) find τ m+1-θ and u m+1-θ (2a) with u m and u m+θ known, u m+1-θ can be compute explicitly from:

u m+1-θ = 1 -θ θ u m+θ - 1 -2θ θ u m (3.43)
(2b) with τ m+θ , u m+θ and u m+1-θ known, solve the tensorial transport problem for τ m+1-θ :

τ m+1-θ -τ m+θ (1 -2θ)∆t + γ m+1-θ = 1 λ (2ηD(u m+θ ) -τ m+θ ) (3.44) with γ m+1-θ = (u m+1-θ .∇)τ m+1-θ -grad(u m+1-θ ) τ m+1-θ -τ m+1-θ grad(u m+1-θ ) T (3.45)
The tensorial transport problem requires a special treatment in FEM which is commonly solved by the discontinuous Galerkin method. This is already implemented in Rheolef and is explained in detail in Saramito [2020, Ch. 4].

Step 3 (t + (1 -θ)∆t .. t + ∆t) find τ m+1 , u m+1 and p m+1 :

This step solves the same problem as step 1, hence by replacing the indexes m and m + θ by m + 1 -θ and m + 1 respectively in eq. 3.38 -eq. 3.41.

Linear viscoelastic ice shelf problem

In the this section the specific problem of this work is presented. First the simplified modelling domain of an ice shelf and its mesh are introduced (subsection 3.3.1), followed by the problems specific boundary conditions (subsection 3.3.2).

Model domain

The Here, the complete length of an ice shelf is explored, starting upstream at the grounding line all the way downstream to the calving front. The grounding line (left hand side in Fig. 3.4) is the location where the ice shelf starts to be detached from the bedrock (brown area in Fig. 3.4) and starts to float. The calving front (right hand side in Fig. 3.4) is the outer end of an ice shelf, facing the open ocean and where icebergs break off recurrently. The geometry of the simplified ice shelf is designed to be at floatation equilibrium. This means that the vertical placement of the ice shelf ramp is determined by the shelf thickness ('H') and the buoyancy force exerted by the water. Hence the ratio of ice and water density prescribe the amount of submergence below sea water level ('h swl ' in Fig. 3.4).

h swl = ρ i ρ w H (3.46)
where h swl is the depth below sea water level, ρ i is the ice density, ρ w is the water density and H is the thickness of the ice shelf.

For FEM the geometry needs to be further discretised into a mesh. Here I have used a common 3-node triangular element mesh, as shown in Fig. 3.4, with a varying element sizes over the length of the shelf. The general mesh set up consisting of an equal amount of elements distributed over the thickness at grounding line and calving front (see Fig. Considering the large difference between the length of the geometry and its thickness, as well as the difference between the thickness at grounding line and thickness at calving front, a uniform element size seemed insufficient, because it would cause a very large mesh, in terms of amount of node in the mesh, which would cost a lot of computational power to compute.

Based on literature (e.g. [START_REF] Reeh | On The Calving of Ice From Floating Glaciers and Ice Shelves[END_REF]) and planed model set up it is to be expected to have a bending of the downstream part of the ice shelf, near the calving front (see Fig. 2.4), due to the difference in lithostatic and hydrostatic vertically distributed pressure forces, whereas no large effects are to be expected near the grounding line, because the grounding line migration is not allowed in this scenario.

This calls for a finer mesh size near the calving front to account for the expected deformation, while a less fine mesh size would be suitable near the grounding line.

Hence the approach of equal mesh density (late referred to as 'MD'), equal amount of elements per shelf thickness at the grounding line and calving front, to balance a fine enough mesh in places where needed and small enough mesh to keep the computational cost low.

Boundary Conditions

In this investigation the spread of the ice without any additional inflow is examined. Hence, the inflow at the grounding line is expressed as a zero horizontal velocity (u x = 0). To remove an additional degree of freedom the ice shelf was pinned to the bedrock at the bottom left point, representing the grounding line where the ice can not move into bedrock. The vertical velocity u y is therefore set to zero in this point. The tensorial transport problem in step 2b of the θ-scheme algorithm (see subsection 3.2.3) is solved using the discontinuous Galerkin method.

The implementation in Rheolef required the shear stress (τ ) to be defined at the upstream (inflow) boundary. In this investigation, I defined the shear stress at the grounding line to be zero (τ upstream = 0). The atmospheric pressure is neglected in this scenario, leaving the top boundary stress free. Further, I do not solve for the evolution of the geometry, hence the geometry is kept constant. The interplay of these two conditions at the water boundary needs to be represented. The pressure of ice is indirectly included in the momentum conservation equation (eq. 3.17), leaving the water pressure to be defined as Neumann boundary condition. Since pressure act normal on the surface the equation writes:

f b = σ.n = p w .n (3.47)
The water pressure acts only on the face below sea water level on the right boundary, and on the whole bottom boundary, which is below sea water level. Hence, the boundary condition for the right and bottom boundary writes This plunging effect shows that the assumption of hydrostatic equilibrium does not hold due to the models dependency to the initial geometry and results in nonrealistic velocities. The geometry is not allowed to deform over time, hence the vertical displacement is not considered and therefore, the water pressure at the bottom of the shelf is not able to balance the weight of the ice. This is also a known problem in viscous Stokes models [START_REF] Durand | Full Stokes modeling of marine ice sheets: influence of the grid size[END_REF]].

p w = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0 if y ≥ swl ρ w gh if y < swl
The problem was resolved by introducing a new vector field to account for and track the displacement (ξ). The displacement is computed after each θ-scheme iteration (subsection 3.2.3) over the same time step (∆t) and its vertical component predictions [START_REF] Reeh | On The Calving of Ice From Floating Glaciers and Ice Shelves[END_REF]].

This approach is neglecting the distortions of the elements, which could only be solved by an adaptive mesh. all models have the top boundary defined as freeboard with an implicit definition of σ.n = 0 and an inflow velocity of 0 m a -1 over the whole thickness of the shelf at the upstream boundary (grounding line). The definition differs at bottom and calving front boundary. While all models have a free-slip condition on the bottom boundary, the shallow shelf approximation implicitly assumes the hydrostatic pressure in normal condition on both boundaries and the viscous and viscoelastic model have the Robin boundary conditions ((eq. 3.47, 3.48 and 3.50)), describing the buoyancy force with displacement adaptation, implemented. Additionally, the SSA equations used from [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF] are not fixed at the grounding line, while the viscous and viscoelastic model are.

The viscoelastic model was run for 6000 iterations (∼ 100 h) to ensure convergence (subsection 3.2.1). All the solutions experiences an over all thinning, which is most dominant at the grounding line (left), and an outward extension at the calving front (right). While the SSA solution appears to be floating in the bounding box, the viscous ('v') and viscoelastic ('ve') solution are at the bottom of the bounding box. This is due to them being fixed at the grounding line on the bottom corner. The viscous ('v') and viscoelastic ('ve') model are only at floatation at the calving front, while the SSA is fully at floatation. Further, the viscoelastic model, different to the other solutions, shows a damped oscillation profile starting from the calving front, which was expected due to the bending moment exerted at the calving front by the buoyancy force (subsection 2.2.3) and under consideration, that the viscoelastic model under constant forcing requires a transient analysis approach, hence it evolves over time. The vertical velocity component of SSA has been excluded, because the fact that it was not pinned at the grounding line has a large effect on the vertical velocity and is therefore not comparable to the other models. The horizontal velocity component appear very similar in the profile representation (Fig. 4.2 -left), with a continuous increase over the length of the shelf. Starting with 0 m a -1 at the grounding line, as prescribed, and peaking at the calving front The vertical velocity components show an overall lower velocity compared to the horizontal velocities by about 2 orders of magnitude and differ largely between the models. The viscous solution shows a zero vertical velocity at the bottom boundary, which decreases vertically with distance to the bottom surface. It has its lowest value with -6.4 m a -1 at the top corner at the grounding line, which is where the shelf is the thickest. This effect is related to the viscous flow, hence thinning of the shelf and due to the shelf being fixated at the bottom corner at the grounding line, which does not allow uplift as adaptation to floatation (see Fig.

Results

Summary

In this section I have presented a comparison between the analytical Shallow Shelf Approximation, the linear viscous Stokes model and the linear viscoelastic Maxwell model under constant forcing. All solutions show a similar overall velocity, which had a closer fit at the horizontal velocities than on the vertical velocities.

The viscoelastic model showed an expected damped oscillation of the velocities, due to the bending moment. Since the velocity of the viscoelastic model oscillations around the velocities of the viscous model, we can assume that the viscoelastic model performs well.

Model sensitivity to geometry

This chapter explores the sensitivity of the previously introduced linear viscoelastic flow model (section 3.3) implemented in Rheolef with regards to mesh resolution and geometrical setting, to specifically analyse the geometry's impact on the model solution as well as model accuracy.

In this analysis I vary the input geometries, which are presented in each sub- For each geometry configuration, the model was run with the same material properties and a constant forcing for 6000 time steps to retrieve a near steady state result (subsection 3.2.1). Each model run started with a time step of ∆t = 2×10 -5 s, which was then increased to 60 s based on good convergence of the residual term.

On average it took 30 iterations for the model to ramp up to and stabilise at ∆t = 60 s.

Sensitivity to mesh resolution

Here we study the sensitivity of the model to the mesh resolution. For this analysis an ice shelf ramp geometry with a length (L) of 5 The horizontal velocity component (u x -red) increases with distance to the grounding line starting from the prescribed value of 0 km a -1 at the inflow. While the rectangle geometry with 400 m thickness at the calving front shows a close to linear increase in velocity in x-direction with a steeper increase in the last ∼ 1500 m, the other geometries (ramp profiles) show a damping behaviour towards the calving front, which is stronger with decreasing shelf thickness at the calving front.

The vertical velocity component (u y -blue) is close to constant with a minimalistic increase towards the calving front and an oscillation in the frontal ∼ 6500 m.

The oscillation is the strongest for the rectangle shaped geometry (400 m) and can be related to the bending moment caused by the pressure difference between water and ice at the calving front (see subsection 2.2.3). A higher calving front causes a larger moment, due to the larger depth of the ice and therefore larger pressure.

In Fig. 4.9 we plot the maximum velocity as a function of shelf thickness. The

Sensitivity to ice shelf length

For the analysis of the length (L) impact on an ice shelf only grounded at the grounding line, the thickness at grounding line (H gl ) is set to 400 m and we keep a constant surface slope (α) of 0.002 rad between the simulations, to compare the same ice-shelf with different calving front position. According to the analytical SSA solution of a viscous model, the velocity should be independent of the shelf front position. As shown in Tab. 4.4 here I look at different length profiles varying from 1 km to 15 km. The vertical velocity component (u y -blue) indicates a very small linear increase towards the calving front and an oscillation in the frontal part, which shows a larger amplitude with shorter length. The oscillation and its amplitude dependency on the shelf thickness at the front have already been discussed in the previous subsection 4. 2.2). Due to the shelf being thicker with shorter length in this scenario, because of the chosen fix surface slope, the oscillation amplitude is larger with shorter shelves. For short ice shelves, like the 1 km geometry used, the length is too short to compensate for the induced bending moment at the calving front and is likely to show its impact further upstream on the grounded glacier.

The horizontal velocity component of the viscoelastic model is close to the SSA assumption of being independent from then shelf length. Yet the vertical velocity component differs largely near the calving front.

With the prescribed model configuration it was not possible to compute a solution for unpinned ice shelf geometries longer than 15 km

Summary

In this chapter, I have shown that the geometry as well as the mesh configuration do have an impact on the results of the model.

For a very simplified case of an unpinned ice shelf only considering the effects of the water pressure at the bottom and front of the shelf, as presented here, a simple mesh of MD2 can be used. But when considering other, specifically local effect like a loading on top of the shelf, a finer mesh, especially near area of local force impacts will be necessary.

The resulting flow velocities show a strong dependency on the shape of the chosen geometry. The horizontal velocity increases nearly linearly with distance to the grounding line but shows a small decrease followed by a stronger increase in velocity near the calving front. While the overall horizontal velocity shows a steeper increase in velocity with larger shelf thickness at the calving front compared to thinner shelves, the vertical velocities show an overall slow increase with a dominant damped oscillation near the calving front. The oscillation amplitude and wavelength is larger with thicker calving front. The very short shelf example of L = 1 km shows a strong decrease in vertical velocity instead of an oscillation profile. A strong relationship between horizontal and vertical velocity component and the viscosity can be seen in Fig. 4.12.b+d, respectively. In fact, every viscositygeometry pair (with varying shear modulus) appears to have the same (with potentially very small variations) velocity component value (u x and u y ). This is supported by the neighbouring shear modulus graphs (Fig. 4.12.a+c), where all velocity component values are distributed in a horizontal line, aligning with a viscosity value.

Parameter study results

Overall, there seems to be no visible impact by the shear modulus on the velocity component values (Fig. 4.12.a+c). A strong inverse exponential correlation between horizontal velocity component and viscosity can be seen in Fig. 4.12.b, where the 'ramp' (orange) geometry shows larger velocities at smaller viscosity compared to the 'rectangle' geometry. The vertical velocity component (Fig. 4.12.b) shows a exponential relationship as well, except for the smallest viscosity. But this cor-relation is not very meaningful, due to the shift in vertical velocity field, as shown in Fig. 4.13. The point was chosen on the calving front at sea water level, but as apparent in Fig. 4.13 and Fig. 4.14, picking a representative point is difficult. no apparent differences between solutions of varying shear modulus, as shown in Fig. 4.12. Tables with all velocity component profiles can be found in Appendix A.

The vertical velocity component profiles (Fig. 4.13) show the largest and most complex changes with viscosity. For very low viscosities, the profile looks similar to a viscous profile. Where the largest downward (negative) velocities are near the surface at the grounding line, and the largest upwards (positive) velocities are along the bottom and calving front boundaries. For the largest investigated viscosities, a velocity profile constant with depth can be seen, with vertical velocities near 0 at the grounding line and the largest downward (negative) velocities near the calving front. Vertical velocity profiles in between the smallest and the largest viscosities show a transition from one state to the other, where the velocity field shows an oscillation pattern, who's maximum velocity moves with increasing viscosity from near the calving front to the grounding line. Two maximum velocity peaks can be seen for viscosities of η = 10 13 Pa s and η = 10 14 Pa s, and it appears that at η = 10 15 Pa s, the second peak is just beyond the grounding line. Additionally, a generally negative velocity exactly at the calving front is observed. This suggests a damped oscillation of the shelf due to the bending moment that is induced by the buoyancy force at the calving front (see subsection 2.2.3).

The horizontal velocity profile (Fig. 4.14) shows a near constant horizontal velocity with depth and an increase with distance to the grounding line for lower viscosities, which is typical for a viscous plug-flow. For viscosities larger η = 10 14 Pa s a continuous change with increasing viscosities of the horizontal velocity field can be observed, so that the largest velocity is experienced only at the top right corner and velocities at the bottom boundary are slower. This is very likely related to the strong impact of the bending moment at high viscosities, as seen in Fig. 4.13 for the vertical velocities. This strong downward bending over the whole shelf causes slower velocities in horizontal direction at the bottom and higher velocities at the top boundary. The remaining model and geometry configurations are kept as in subsection 4.3.1

and each geometry ('ramp' and 'rectangle'), material (Fig. 4.11 and Fig. 5.5) and tidal (Tab. 5.1) parameter combination was run for a time period of 5 days. As common for time-depended problems, the model was started from a (near) steady state solution (subsection 3.2.1), in this case the solutions of the parameter study under constant forcing (subsection 4.3.1) were used.

Model sensitivity to tidal oscillation

In this section, I compare the impact of different tidal forgings (see Tab. the whole domain -follows the sinusoidal velocity responds at the calving front closely, but is slightly higher/lower, and hence shows that the maximum vertical velocity is not at the calving front in any of the 4 cases. In Fig. 5.3.c,d,g,h it can be seen that the maximum vertical velocity is high at the calving front but its maximum is at around 4 km from the grounding line. All 4 tidal configurations show a very similar maximum velocity distribution pattern in Fig. 5.3, which only By comparing the velocity profiles for different shear modulus, no difference in velocity pattern or velocity values can be found. This highlights again that the shear modulus does not affect the model solution, only the successfulness of finding a solution.

The viscosity, on the other hand, strongly effects the velocity pattern and the location of the largest velocities. The largest horizontal velocities can be found at the shelf front for small viscosities, while viscosities of 10 14 and 10 15 Pa s cause the largest velocities to be near the grounding line. The velocity maximum shifts to occurring only at the top surface near the calving front for larger viscosities.

The vertical velocity component profiles show always high velocities at the calving front but the peak velocity moves slowly from the grounding line towards the calving front with increasing viscosities, leaving a rapidly decreasing velocity trace behind. stress is equal to the viscous stress, hence there is no elastic effect, the circle forms a perfect 45 °line, which is indicated as a thin dotted line in the graphs.

The graphs show that both the total and viscous stress are very similar. The biggest difference can be seen at the lowest investigated viscosity (η = 10 10 Pa stop row), where the stresses at the calving front (10 km -darkest colour) do not form a perfect line, but appear slightly elliptic and tilted off the 45 °line. The stresses at the grounding line (0 km -lightest colour) form a slim ellipse that is also not centred around the 45 °line, but the tilt is less than for the stresses at the calving front. While all the stresses in-between these locations form a line on the 45 °line. Also, noticeable is that the stresses at the calving front cycle through 0 and are sometimes positive and sometimes negative. But the closer the location is to the grounding line, the larger the stresses are -positive for horizontal and negative for vertical stresses. The stress variation (circle diameter or line length in the graphs) is the largest at the calving front, second largest at the grounding line and for the in-between locations the stress variation is comparably small.

The stresses at the viscosity of η = 10 14 Pa s (mid row) cycle through 0 but are not exactly centred. All stresses except the ones at the grounding line appear to form a line in a 45 °angle, while the stresses at the grounding line are largest and form a slim off centric ellipse.

The stresses at the highest investigated viscosity of η = 10 18 Pa s (bottom row) appear to cycle through and are centred around 0. While all stresses form lines they appear to be slightly tilted to the 45 °line. The stress variation decreases with distance to the grounding line.

Overall, there is a strong increase in stresses with increase in viscosity.

Long time study

Following the slow decrease of the horizontal velocity over time for lower viscosities, especially in the rectangular geometry, in the previous section 5. The horizontal velocity experiences a fast decline from ∼ 10000 m a -1 to ∼ 2000 m a -1 within the first ∼ 60 days, and continues to decrease at a lower rate for 2.5 -3 years, where it reaches a mean velocity of ∼ 0 m a -1 . After the 3 year mark, it starts to slowly increase in velocity again. During the whole study period, the horizontal velocity appears to have a constant velocity range of ∼ 1500 m a -1 .

Overall, these results suggest that there might not be a solution where stresses are at equilibrium, or it takes many more simulated years to reach that point. But essentially, it highlights that future viscoelastic studies with low viscosities, need to consider long model ramp up times.

Summary

In this chapter, I have studied the linear viscoelastic model with different simplified tidal configuration and analysed these with different material parameters.

The tidal amplitude and tidal cycle ('diurnal' and 'semi-diurnal') influences the model's vertical and horizontal velocity components, both in velocity amplitude and oscillation wavelength.

A material parameter study with tidal forcings suggests the existence of three different flow regimes when looked at the velocity components within the investigated parameters. Three patterns (I, II, III) have been found for the horizontal velocity changing with increasing viscosity. Fitting to this is the encountered phase change (pattern IV) between the vertical velocity and water level velocity.

Small viscosities have been shown to have large horizontal velocity, that decrease over time. A nearly 5 year time period investigation has not produced a constant horizontal velocity oscillation. It is not clear if there is an equilibrium state.

Again, as in the previous chapter, we could only find a shear modulus influence on the model success rate, but no effect on the actual solution. The model success rate appears to be driven largely by viscosities, because model runs with low viscosities, were generally quite successful. But larger viscosities did only produce results for very high shear modulus.

The shear stress comparison suggests an overall viscous dominant flow, with some elastic relevancy at lower viscosities.

Chapter 6

Discussion & Conclusion

The aim of this thesis was to develop an open souce viscoelastic model for ice shelf dynamics. This was done in a finite element method environment called Rheolef (subsection 3.2.2), which already had the major functionalities, i.e. solvers for viscous and viscoelastic flow, implemented. The main task here was to understand the general viscoelastic problem (mathematically and numerically), to understand the existing code and to modify it according to our needs to represent the dynamics of ice shelves.

The results of the developed linear viscoelastic Maxwell model were compared in terms of geometry and mesh size sensitivities as well as material parameter sensitivities. The later was studied with constant forcing and simplified tidal forcing.

The key findings are:

• The study demonstrates that the viscoelastic model is able, compared to a viscous model with the same boundary conditions, to capture the vertical oscillation induced by the bending moment at the calving front, as well as the oscillation induced by tidal forcing.

• The data suggest that the amplitude and wavelength of vertical velocity oscillation are largely governed by the viscosity and the geometry, i.e. thickness at calving front.

• The results also indicate that the same parameters, viscosity and geometry or shape, drives the horizontal velocity (outflow).

• The study indicates that the oscillation amplitude and wavelength of the velocities induced by the tidal forcing, varies with distance to the grounding line.

• The study suggests that there are different flow types, which are primarily driven by the material's viscosity.

• Further, the material parameter analysis demonstrates that convergence of the model is determined largely by the shear modulus, governed by the Young modulus and Poisson's ratio, for a steady state scenario; and by the shear modulus and the viscosity for the tidal forcing scenario.

• The study indicates that the shear modulus does not affect the solution, i.e. velocities.

The comparison between the analytical solution (SSA), a viscous model and the viscoelastic model (section 4.1) showed an overall similar horizontal velocity increase over the length of the shelf, with each model peaking at the calving front with ∼ 120 m a -1 . SSA over-estimates the horizontal velocity by up to 2 m a -1 compared to the viscous solution, and the viscous solution over-estimates the velocity by up to 4 m a -1 at the calving front compared to the viscoelastic solution. Considering the overall horizontal velocity, these differences appear very small.

The vertical velocity difference between the viscous and viscoelastic solution is largely due to the bending moment. [START_REF] Reeh | On The Calving of Ice From Floating Glaciers and Ice Shelves[END_REF] suggested that the pressure imbalance at the calving front induces a bending moment, which in turn would be visible in the velocity field. Since both models are compared under constant forcing, the viscous problem reduces to a time independent problem, where the evolution of the deformation over time is not considered. While the viscoelastic problem stays time dependent and solving for a near "steady state" requires an iterative computation over time. [START_REF] Mosbeux | Viscous and elastic buoyancy stresses as drivers of ice-shelf calving[END_REF] investigated the buoyancy stresses at the calving front of ice shelves with a viscous and an elastic model, and demonstrated a instantaneous bending for the elastic model and a delayed bending for the viscous model. Hence, since our viscous model does not evolve over time, we do not see the delayed impact of the pressure imbalance at the calving front and therefore do not see a bending of the shelf front. While the time dependent viscoelastic model produces a damped oscillating profile, similar as suggested by [START_REF] Reeh | On The Calving of Ice From Floating Glaciers and Ice Shelves[END_REF], that appears to oscillate around the viscous solution, which suggest, based on the mentioned previous study and considering the short time scale, that we see instantaneous elastic bending.

This gives rise to the assumption that the introduced linear viscoelastic model does give reasonable results for ice.

The investigation of different shelf geometries under constant forcing (section 4.2)

showed that the shape and size of the ice shelf profile has a large influence on the modelled ice flow velocity. The observed over all increasing horizontal velocity with distance to the grounding line as well as the horizontal velocity increase rate with shelf thickness at the calving front, can be related to the viscous thinning as has been shown in purely viscous models previously, e.g. [START_REF] Greve | Dynamics of Ice Sheets and Glaciers[END_REF].

The damped increase in horizontal velocity in the mid section of the ice shelf and its followed stronger increase at the calving front as well as the oscillation of the vertical velocity near the calving front, can be attributed to the imbalance of the hydrostatic and lithostatic pressure. This imbalance at the calving front causes a moment, which in turn causes a bending in a damped oscillated form, as proposed by [START_REF] Reeh | On The Calving of Ice From Floating Glaciers and Ice Shelves[END_REF] and shown in Fig. 2.4.

The increase in oscillation amplitude and wavelength of the vertical velocity with increasing shelf thickness at the calving front, can be connected to the before mentioned bending moment as well. A larger depth below sea water level increases the water pressure and hence the moment. Therefore, the bending moment increases with increasing shelf thickness.

The shortest investigated ice shelf (1 km) has shown a larger increase in horizontal velocity and a low vertical velocity without oscillation compared to longer shelf profiles. It is likely that the impact of the bending moment at calving front could not be compensated for and will be noticed upstream beyond the grounding line. [START_REF] Gudmundsson | Fortnightly variations in the flow velocity of Rutford Ice Stream, West Antarctica[END_REF]] observed such an effect on the velocities caused by the tidal movement upstream of the grounding line (see Fig. 2.5).

Unfortunately, the model with current parametrisation was not able to produce results for shelves longer than 15 km. I think that this is due to the 2D-profile and the free floating nature of the shelf geometry. The usage of a more realistic geometry, e.g. 3D-profile, which would introduce buttressing on the sides, or the introduction of pinning points, could stabilise the problem and therefore allow for longer shelf geometries.

This investigation used highly simplified geometries, yet it highlights the importance of using a representative geometry length and thickness profile. A 'rectangle' shaped geometry, hence a constant thickness, as assumed in the thin-beam or thinplate theory might be a good simplification for investigating the ice shelf front, but will potentially over-estimate velocities if considered from the grounding line.

In section 5. The vertical velocity components of all 4 configurations showed a velocity change with time that followed the induced water level velocity neatly at the calving front.

While the vertical velocity at the calving front is large, the maximum velocities were found to be near the middle of the shelf. The horizontal velocity component showed that the velocity at the calving front oscillates proportionally with the tide.

But it also showed that the faster tidal changes of the semi-diurnal cycle can not be compensated by the material elastic modulus, causing large horizontal surface velocities at the top and bottom of the shelf near the fixed grounding line. This is most dominant for the extreme cases of 5 m amplitudes, which are much larger to the tidal amplitudes found in Antarctica (upto 1.5 m). Additionally, the fixed shelf to the grounding line, causes strong bending near this location, where in reality it would move up and down with the tide. Hence those velocities would not be found on a real ice shelf and are artificially due to the model configuration. On the other hand, it suggests, that in a more realistic scenario, where the grounding line is not fixed, velocities of the ice shelf align with the tidal oscillations and might be strongest were movement is restricted or these velocity oscillations might be transported upstream of the grounding line, which was seen in observations at the Rutford ice stream [START_REF] Gudmundsson | Fortnightly variations in the flow velocity of Rutford Ice Stream, West Antarctica[END_REF]].

The study of material parameters was presented for the viscosity and the shear modulus under constant (section 4.3) and tidal forcing (section 5. The viscosity appears to not only drive the velocity, larger viscosity results in smaller velocities, but in the case of the viscoelastic model also appears to drive the overall flow regime, which can be identified by the phase shift between input and output velocities. In accordance to other studies we would be looking at a viscous dominant flow at lower viscosities (10 10 -10 12 Pa s), a viscoelastic flow at medium viscosities (10 13 -10 15 Pa s) and an elastic dominant flow at higher viscosities (10 16 -10 18 Pa s). [START_REF] Wild | Viscosity and elasticity: a model intercomparison of ice-shelf bending in an Antarctic grounding zone[END_REF]; [START_REF] Walker | Ice-shelf tidal flexure and subglacial pressure variations[END_REF] found that viscosities of 10 16 Pa s and 10 17 Pa s respectively correspond to a purely elastic response and viscosities of 10 13 -10 15 Pa s to a viscoelastic regime. Yet, the total shear stress to viscous shear stress comparison suggests an overall viscous dominant flow and slightly stronger elastic influence at lower viscosities. This would explain the lack of influence of the elastic modulus on the solution that was mentioned earlier, but raises the question on where the phase shift between the velocities comes from.

Published observations of the daily tidal cycle by [START_REF] Robel | Tidal modulation of ice shelf buttressing stresses[END_REF], using the same GPS measurements as [START_REF] Gudmundsson | Fortnightly variations in the flow velocity of Rutford Ice Stream, West Antarctica[END_REF] (see Fig. The mesh density approach was originally chosen to ensure enough nodes near the calving front, where an impact of the bending moment was suspected. A mesh density of 2 elements per thickness appeared to give a reasonable small error of ∼ 1 m a -1 in horizontal velocity compared to the finest investigated mesh of MD 10 and under consideration of the computational cost. The comparison was done at sea water level and in hindsight I realized that this is the point of the smallest impact of the bending moment induced by the pressure imbalance. Hence the error is potentially larger. Also, only 3 nodes at the calving front appears now rather few to capture the impact.

In the Robin boundary condition only the vertical displacement is considered, and not the horizontal. This can potentially cause an uncertainty in the solution.

For the elastic part only the shear modulus was considered, which becomes smaller with larger Poisson ratio. The fact that the bulk modulus was neglected might be one reason why I was not able to detect an significant difference in the solutions with different shear moduli.

The assumption of fixing both models to the grounding line was done to remove a degree of freedom and to enable the model to find a solution. This assumption holds for static (unchanging) geometries at short time scale investigation, but is not ideal, especially when looking at tidal displacements, and a different way should be found, e.g. including the upstream part of a glacier, using the 3D capabilities of Rheolef and include buttressing from the sides, or include pinning points.

In conclusion, it was demonstrated in this work that the linear viscoelastic model, developed in the open source program Rheolef, can simulate viscous thinning and elastic bending of an ice shelf induced both by the pressure imbalance at the calving front as well as by tidal oscillation of sea water level. Even so, a clear classification of the flow was not possible due to contradictory results in velocities and shear stresses.

The model with current parametrisation was not able to produce results for shelves longer than 15 km, which may be improved by including e.g. pinning points and/or buttressing.

The numerical scheme fails at some points within the parameter envelope for ice.

The model only produced results under tidal forcing for a realistic shear modulus of 1 -2 × 10 9 Pa and a viscosity of 10 13 Pa s but not for 10 14 Pa s, hence it does not cover all reasonable material parameters for ice. This could be related to other parameters, like the chosen mesh, geometry or boundary conditions.

Overall, and especially with the tidal forcing experiments, I showed that both viscous and elastic strains should be taken into account for ice shelf investigations.

Chapter 7

Outlook

As a next step I would make use of the 3D capability of Rheolef. It would allow for more realistic ice shelf geometries, where buttressing from the sides and pinning points can be accounted for, and would potentially allow to solve for larger shelves. Additionally, I would include Glen's flow law, to consider more realistic ice rheology. With this the impact of several scenarios could be studied, e.g. the impact of basal crevasses, meld ponds and underwater protrusions.

I further would suggest to use a larger domain that includes both ice stream and ice shelf. With this the fixation of the shelf to the grounding line could be removed and would allow the grounding line to migrate. At which point the investigation of the full tidal cycle, the consideration of the 2-weekly tidal variation and not just the daily oscillations, could be interesting.

Additionally, a deeper investigation on the elastic parameters and its impact on the model solution, as well as an investigation of the dependency of the model's convergence to the material parameters should be done. Further, a study of the phase shift between e.g. tide and solution, and therefore the dominant flow regime would be insightful.

Following this, one could think of introducing a moving mesh to fully model shelf deformation and allow for studies at longer time scales.

If, specifically the study of more realistic case does give acceptable results an, an investigation of alternative numerical schemes needs to be considered.

à

  des codes commerciaux coûteux aux algorithmes figés et non adaptables aux spécificités du problème. Cette thèse propose le développement d'un modèle viscoélastique de la dynamique des shelves à l'aide de l'environnement 'open 'source' Finite Element Method C++ Rheolef. v

Figure 1 . 1 :

 11 Figure 1.1: Elements of the Cryosphere -with observed changes summarised in the 5 th IPCC Assessment Report in red. Source: Vaughan et al. [2013, Fig.4.25] (adapted)

Figure 1 . 2 :

 12 Figure 1.2: Marine Ice Sheet Instability (MISI) -schematic of the grounding line retreat due glacier and ice shelf thinning on a retrograde bedrock slope. Source: Pattyn [2018]

Figure 2 . 1 :

 21 Figure 2.1: Overview map of Antarctica. Source: LIMA Project [BAS 2007] (adapted).

Figure 2 . 2 :

 22 Figure 2.2: Recent ice shelf collapse events. (left -source: Rott et al. [1996]) Larsen A ice shelf ERS-1 image taken on 16. February 1993. It shows the stages of collapse with the retreat of the ice shelf front (numbered lines) that progressed quickly between October 1994 (5) and March 1995 (9). The thick white line beyond line 9 indicates the grounding line. (right -source: Scambos et al. [2004]) Larsen B ice shelf MODIS images taken on the 1. November 2003. It shows the ice shelf extent around 2000-2001 in yellow, the extent after the collapse in 2002 in red and the grounding line in thick black.

Figure 2 . 3 :

 23 Figure 2.3: The buttressing effect of ice shelves in Antarctica. The buttressing is shown as percentage, areas in dark orange (close to 1) show a strong contribution to buttressing of the inflowing glaciers and areas in blue, with little to no buttressing effect are considered passive ice. A buttressing threshold was estimated to calculate the amount of passive ice, given in percentage in blue numbers per ice shelf, coast and total Antarctica. Source:[START_REF] Fürst | The safety band of Antarctic ice shelves[END_REF]

2. 2 . 1

 21 Hydro-fracturing -also known as: water-assisted crevasse propagation Hydro-fracturing describes the process of crack deepening through the accumulation of water in the cracks lowest part. With the increase of water in the crack, the pressure at crack bottom increases and eventually causes deepening of the crack. This process can penetrate an ice shelf in a fairly short time, according to[START_REF] Scambos | Ice shelf disintegration by plate bending and hydro-fracture: Satellite observations and model results of the 2008 Wilkins ice shelf break-ups[END_REF] within half-an-hour. They showed that hydro-fracturing depends on the amount of water in the crack as well as the water density. While seawater (1 028 kg m -3 ) and brine (1 084 kg m -3 ) infiltration can only fill the crack up to sea water level, melt water, i.e. freshwater (∼ 1000 kg m -3 ), can fill the crack up to the surface. According to their model, the crack filling up to sea water level does not cause enough pressure to cause crack deepening.

  2.4.a). Hence, a bending moment is exerted on the ice shelf front as shown in Fig. 2.4.b, forcing the upper part of the cliff forward (in flow direction) and downward. This rotation of the ice shelf front leads to tensile stresses on the shelf's surface (near the cliff)

Figure 2 . 4 :

 24 Figure 2.4: Buoyancy forces acting on the ice shelf front. (a) Demonstration of the pressure imbalance (red) at the ice shelf front, that causes a bending moment on the front. (b) Demonstration of the ice shelf oscillation caused by the bending moment. Both sketches are based on Figures by Reeh [1968].

  (*) -also known as: tidal-flexure zone, grounding line migration Tides, a periodically repetitive change of sea level, are highly dependent on the constellation of the moon and sun relative to the earth. Most tides follow a semi-diurnal (M2, S2) cycle, while in some places a diurnal (O1, K1) cycle is dominant [Rosier et al. 2014; Jourdain et al. 2019]. Both semi-diurnal and diurnal cycles show a repetitive fortnight (Mf, Msf) pattern [Gudmundsson 2006].

Figure 2 . 5 :

 25 Figure 2.5: The impact of tides on the ice flow velocity of the Rutford ice stream. The Rutford ice stream is a glacier feeding the Rønne-Filchner Ice Shelf in West Antarctica. The top panel shows the vertical tidal amplitude measured on top of the ice shelf, 20 km downstream of the grounding line (GL). The bottom panel shows the recorded surface velocities downstream of the GL (on the ice shelf -red), at the GL (green) and upstream of the GL (light and dark blue). The numbers indicate the distance to the GL. Source: Gudmundsson [2006].

  or (kinematic) GPS [Le Meur et al. 2014] for point 'F' in Fig. 2.6 or analysing the slope change in DEMs [Scambos et al. 2007; Bindschadler et al. 2011] for point 'Ib'. Le Meur et al. [2014] used measurements of 'F' and an elastic beam model to retrieve a closer approximation of the grounding line position.

Figure 2 . 6 :

 26 Figure 2.6: Grounding zone schematic, indicating important points in the transition zone between grounded and floating ice. Source: Davies [2022] (adapted).

  2.7). Depending on the bathymetry mainly warm, salty Circumpolar Deep Water (CDW), cold, fresh Antarctic Surface Water (AASW) or cold, salty Shelf Water (HSSW) is pulled in. The dominating water type pulled into the ice shelf cavity affects the amount of melt near the grounding line and refreezing mid-shelf to near the ice shelf front [Lewis and Perkin 1986; Grosfeld and Sandhäger 2004]. Thinning of the ice shelf causes a retreat of the grounding line and hence sea level rise.

Figure 2 . 7 :

 27 Figure 2.7: Subglacial shelf water (HSSW) currents and its affect on an ice shelf, an example of the subshelf-water interface process. Source:[START_REF] Grosfeld | The evolution of a coupled ice shelf-ocean system under different climate states[END_REF]

Figure 2 . 8 :

 28 Figure 2.8: Ice rises and ice rumples -schematic showing the difference between an ice rise and ice rumple.

Figure 2 . 9 :

 29 Figure 2.9: McDonald Ice Rumples on Brunt Ice Shelf, with the ice flowing from the top right corner to the bottom left. The images was taken on the 5th of January 2020 with Sentinel 1. Source: Luckman [2021]

  constitutive equations. The field equations are general conservation laws describing the preservation of physical entities like mass or momentum and combined are best known as the Navier-Stokes equations. While the constitutive equations describe the material specific stress-strain relationship, which is often depicted with the mechanical elements: spring and dashpot (examples in Fig. 2.10).

Figure 2 .

 2 Figure 2.10: Mechanical elements representing the stress-strain relationship of a onedimensional linear fluid flow, where E is the elastic modulus and η is the viscosity. (a) Hook model -elastic flow, (b) Newton model -viscous flow, and in combination presenting viscoelastic flow in varying complexity: (c) Maxwell model, (d) Oldroyd model, and (e) Burger model.

The

  spring and dashpot in combination represent the viscoelastic flow. The most simple combination is one element each in series, known as the Maxwell model (Fig. 2.10.c). Sometimes more elements are combined to better represent the behaviour of a flow, most common here the Oldroyd model (Fig. 2.10.d), where η 1 is the viscosity at macroscopic scale and η 2 is the viscosity of the solvent [Saramito 2016], or Burger model (Fig. 2.10.e), representing instantaneous elastic strain, delayed elastic response and viscous deformation [Gudmundsson 2011]. Based on this, fluids are often categorised in two major groups: Newtonian and non-Newtonian fluids. Newtonian fluids (Newton model -Fig. 2.10.b) with a linear stress-strain-rate relationship can very precisely be described by the Navier-Stokes equations. But unfortunately, the majority of fluids (including ice) fall into the category of non-Newtonian fluids -also known as complex fluids -e.g. viscous fluids with a non-linear stress-strain-rate relationship or viscoelastic fluids (Fig. 2.10.c-e).

Figure 2 .

 2 Figure 2.11: The deformation behaviour of ice over time. (left) Schematic of 2D shear stress τ and its resulting shear angle γ. (right) The shear angle or deformation of ice over time. Source:Greve and Blatter [2009, Fig. 4.4]

  rameter and n the creep exponent. The creep parameter depends strongly on ice temperature, but also on grain size, impurity content and water content. The creep exponent is considered between n = 2 -4 but is commonly set to 3 [Cuffey and Paterson 2010].

3 [

 3 Material specific properties needed to simulate the flow of ice are for example ice density, viscosity and/or elastic modulus. These are often retrieved from in situ or laboratory measurements or by tuning a numerical model to an observation. This subsection gives an overview of the most commonly used parameters:The ice density (ρ i ) ranges from 830 -923 kg m -3 . The lower value refers to the density where interconnecting passageways between grains are cut off, also known as pore close-off, and the existing air inside the ice is trapped. The point defines the boundary between firn and ice. The higher value refers to the maximum known possible density for ice and can be found in ice sheets at mid depth with low temperature and medium pressure [Cuffey and Paterson 2010]. Because ice densities vary with depth, hence pressure, and temperature, average values of 900 kg m -Scambos et al. 2009], 910 kg m -3 [Greve and Blatter 2009] and 917 kg m -3 [Rignot

Figure 2 .

 2 Figure 2.12: The shallow approximations compared to Stokes, showing the influential elements considered in the stress tensor for (a) the Full Stokes equation, (b) the Shallow Ice Approximation and (c) the Shallow Shelf Approximation. Based on: Davies [2020]

  2.13. While the Stokes equation is a 3D problem with 4 degrees of freedom, SIA is reduced to a local problem, only requiring vertical integration at a given location to find the velocity, which makes SIA the most simplified version of the Stokes equations with a very low computational cost in comparison. It works well for ice sheets and glaciers, but not for ice shelves -where there is no basal drag, or fast flowing ice streams -which require consideration of longitudinal stresses.

Figure 2 .

 2 Figure 2.13: Comparison of Stokes and its approximations (top) Shows the segmentation of ice sheet, ice shelf and the transition zone between them and their characteristic velocity profiles. The transition zone includes faster flowing ice streams and the area around the grounding line. (bottom) The three graphs show the complexity of the stress field based on the consideration of basal traction in flow direction. The top graph shows the combination of the most simplest models: SIA and SSA, the middle graphs represents the coupled Hybrid model, while the bottom graph shows the stress transfer for Higher-Order and Full Stokes models. Source: [Nowicki and Seroussi 2018]

  2.10.e), a 4-element viscoelastic model, representing instantaneous elastic strain, delayed elastic response and viscous deformation[START_REF] Gudmundsson | Ice-stream response to ocean tides and the form of the basal sliding law[END_REF]]. They were able to demonstrate a good fit between a linear-viscoelastic beam model and tidal bending. Gudmundsson showed in his 2011 paper that this behaviour of ice can be reconstructed with a simpler Maxwell model (Fig.2.10.c) -2-element viscoelastic model -for loading times greater than 100 s. Rosier, Gudmundsson, and Green moved to a 3-D model[START_REF] Rosier | Insights into ice stream dynamics through modelling their response to tidal forcing[END_REF]] and further to a 3-D Full Stokes model[START_REF] Rosier | Temporal variations in the flow of a large Antarctic ice stream controlled by tidally induced changes in the subglacial water system[END_REF] to better represent tidal flexure around the grounding line. Other scenarios investigated with a viscoelastic model include the ice shelf bending due to filling and draining of melt ponds[MacAyeal et al.

  elastic and creep deformation. Hence fitting an elastic model to this kind of data results in a lower Young modulus compared to laboratory experiments, which was measured at very short time-scales and hence can be considered to capture the purely instant elastic modulus. Gudmundsson was able to fit his Maxwell model with Glen's flow law to Reeh et al.'s more complex Burgers model with a Young modulus of 4.8 GPa while[START_REF] Reeh | Tidal bending of glaciers: a linear viscoelastic approach[END_REF] used an instant elastic modulus of 9.3 GPa and a delayed elastic modulus of 10 GPa, noting that the Maxwell model is a good approximation for loading times larger 100 s but still a simplification.

Figure 3 . 1 :

 31 Figure 3.1: Components of the stress tensor.

  σ yy and σ zz (Fig.3.1 -red arrows) and the 6 deviatoric stress components: σ xy , σ xz , σ yx , σ yz , σ zx and σ zy (Fig.3.1 -blue arrows), where the first index refers to the plane normal to the axis the stress is acting on and the second index refers to the direction the stress in acting on that plane.The normal stress components, also known as principal or hydrostatic stress components, act perpendicular to a plane and relate to the volume change. Positive normal stresses are extensive stresses and negative normal stresses are compressive stresses. The sum over all normal stresses is called mean stress σ m :

  2.a) to be turned into a discrete representation, a so-called mesh (Ω d in Fig. 3.2.b). A mesh consists of a finite amount of non-overlapping elements (Ω k in Fig. 3.2.c) connected to each other at points on their boundaries, so-called nodes.

Figure 3 . 2 :

 32 Figure 3.2: Discretisation of 2-D Geometry into a mesh. Source: Tekkaya and Soyarslan [2014]

  in the variational formulation. These are weakly imposed in the variational formulation of the flow equation and are referred to as Neumann BC. A combination of the Dirichlet an Neumann condition is referred to as Robin BC.

  uses FEM to solve PDEs and is freely available under the Gnu General Public license. It is designed as a set of C++ libraries written mainly by Pierre Saramito, who also wrote a book about complex fluids[START_REF] Saramito | Complex fluids -Modeling and Algorithms[END_REF]], which I have used largely as reference and which contains explanations and proof of most equations related in this work, as well as an extensive collection of benchmark and code examples for Rheolef[START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF]]. I specifically used his 'Stokes problem'

  for viscoelastic fluids introduced by Saramito in 1994. The algorithm is already implemented in Rheolef and explained in detail for the Olroyd-B problem in Saramito [2016, ch. 4.6 ff]. It uses a discontinuous P1 shape function for the field variable p and τ and a continuous P2 for the velocity u. In this section I give an overview of the algorithm for the upper-convected Maxwell problem, which is the reduced Oldroyd problem with η 2 = 0 (see Fig. 2.10).

Figure 3 . 3 :

 33 Figure 3.3: θ-scheme algorithm for time discretisation. Source:Saramito [2016, fig. 4.11]

  model domain for a 2-dimensional ice shelf is represented by a vertical cut along the flowline. The flowline is the central line of the fastest flowing part of an ice shelf that is usually fed by an inflowing glacier upstream and extents to the calving front downstream. Because an ice shelf is flowing towards the calving front with little resistance, the shelf is stretched out and hence experiences thinning. The simplest geometry to represent this behaviour is a ramp as shown in Fig. 3.4.

Figure 3 . 4 :

 34 Figure 3.4: Modelling domain of the 2D ice shelf geometry at floatation equilibrium, cut vertically along the flowline. The light blue ramp is representing the simplified geometry of a floating ice shelf, the darker blue area: water, and the brown area below: bedrock. The left side of the ramp is the grounding line, where the shelf is the thickest (H gl ). The right hand side is the calving front with the thinnest part of the ice shelf (H cf ). 'L' stands for the length of the shelf from the grounding line to the calving front and the flow direction is from the grounding line (left) to the calving front (right). The ice shelf is set at floatation equilibrium to sea water level (swl), causing the surface (top) of the shelf to be sloped towards the calving front at an angle 'α'. The triangles inside of the ramp, show the mesh structure used.

  3.4 for an example of 2 elements per thickness). Since the ice shelf is less thick at the calving front compared to the thickness at grounding line, the elements are smaller near the calving front and larger near the grounding line. The element size changes gradually over the length of the ice shelf.

Figure 3 . 5 :

 35 Figure 3.5: Boundary conditions of the 2D ice shelf , with the inflow velocity u x and the grounding line pinning point defined with no velocity in y-direction (u y ) on the left hand side, and the bouyancy force of water f b acting on the bottom and right boundary. Additionally, the overall acting gravity g and the densities of ice (ρ i ) and water (ρ w ) are represented.

  Fig. 3.6 on the left hand side -for the ramp geometry on the top and rectangle geometry on the bottom.

(

  ξ y ) is then added to h, similarly toChristmann [2017, ch. 3.3], replacing eq. 3.49 with h = h swl + ξ y(3.50) Different to[START_REF] Christmann | Viscoelastic Modeling of Calving Processes at Antarctic Ice Shelves[END_REF], who uses the displacement at the end of each time step, I use the trapezium rule, as simplification of the time integral over u, for a better approximation. The current displacement is then added to the previous displacement.ξ new = ξ old + ( u m + u m+1 2 ∆t) (3.51)where u m is the velocity field at the beginning of the time step and u m+1 the velocity field at the end of the time step ∆t. Due to the displacements dependency to the unknown variable of velocity, the displacement function transforms the boundary condition from a Neumann condition into a Robin condition.

Fig. 3 .

 3 Fig. 3.6 on the right shows the viscoelastic model results with the described Robin boundary condition (eq. 3.47, 3.48 and 3.50). The deformation shown in the figure is not the deformation of the mesh, but shows what would have been obtained by multiplying the velocity field with a fictive time step. In the comparison with the Neumann condition results, the Robin condition shows maximum velocities to be 1 magnitude smaller and more realistic (based on the comparison to the analytical viscous model presented in the next chapter 4) than the results with the Neumann condition. The deformation of the Robin condition results shows a floating shelf with a downward bending of the shelf front and a small damped oscillation of the shelf as expected from observations and

Figure 3 . 6 :

 36 Figure 3.6: Boundary Condition (BC) comparison -Neumann versus Robin. The results of the Neumann BC set up are shown on the left, and the results for the Robin BC set up on the right. The top images represent the ramp geometry with L = 10 km, H gl = 400 m and H cf = 200 m, and the bottom images the rectangle geometry with L = 10 km and H gl = H cf = 400 m. The colouring of the four images represent the vector norm of the velocity (u) in m a -1 , with darker colours representing a faster flow. Note: each image has its own scale. The length of the geometry was scaled by 0.02, and the deformation, which is based on the velocity vector, was enhanced for better visibility by 0.02 and 2 for the Neumann case and Robin case respectively.

Fig. 4 .

 4 Fig. 4.1 shows an overall similar norm velocity (light to dark colour), which increases towards the calving front in all three solutions, peaking at 120 m a -1 . The thin black lines in each figure is the bounding box of the original ramp profile, hence it represents the original thickness and length of the ramp profile but not its exact shape. The solution geometry shown is the deformation, which was retrieved from the velocity vector in Paraview and was exaggerated by the factor 5.

Figure 4 . 1 :

 41 Figure 4.1: Model comparison by velocity between the Shallow Shelf Approximation ('SSA' -top), viscous ('v' -middle) and viscoelastic ('ve' -bottom) models for a the ramp geometry with L = 10 km, H gl = 400 m and H cf = 200 m. The colouring of the four images represent the vector norm of the velocity (|u|) in m a -1 , with darker colours representing a faster flow. The length of the geometry was scaled by 0.02, and the deformation, which is based on the velocity vector, was enhanced for better visibility by 5.

Fig. 4 .

 4 Fig. 4.2 presents the same solutions from Fig. 4.1 separated in its velocity components and shown in its undeformed state, while Fig. 4.3 shows plots of the velocity components at the top and bottom surface, as well as at sea water level.

Figure 4 . 2 :

 42 Figure 4.2: Model comparison by velocity components (2D profile) between the Shallow Shelf Approximation ('SSA' -top), viscous ('v' -middle) and viscoelastic ('ve' -bottom) model for a the ramp geometry with L = 10 km, H gl = 400 m and H cf = 200 m. Shown are the solutions of the horizontal (u x -left) and vertical (u y -right) velocity component. The length of the geometry was scaled by 0.02 for better visibility.

Figure 4 . 3 :

 43 Figure 4.3: Model comparison by velocity components (graph) between the Shallow Shelf Approximation ('SSA'), viscous ('v') and viscoelastic ('ve') model. Shown are the solutions of the horizontal (u x -left) and vertical (u y -right) velocity components per model at the 'top surface', at 'sea water level' and at the 'bottom boundary'.

  4.1). The viscoelastic solution shows a distinct damped oscillation of the vertical velocity. This oscillation is the strongest at the calving front and decreases in amplitude towards the grounding line. It appears that the vertical velocity oscillates around the viscous linear velocity.

Fig. 4 .

 4 Fig. 4.4 and Fig. 4.5 show the difference between two solutions, which is presented as their velocity components. In Fig. 4.4 negative values are shown in blue and positive values in red.The top rows in both figures show the difference between the viscous and the SSA solution. Overall, the velocities computed with SSA are either higher or equal to the viscous solution. The horizontal velocity difference is zero at the grounding line, as prescribed as inflow velocity, and increases with distance to the grounding line. The largest difference can be found at the bottom corner at the calving front, where the SSA solution overestimates the velocity by 2.4 m a -1 (∼ 2 %). The horizontal velocity difference between the two solutions appears to change equally over the shelf's thickness and only diverges within a half kilometre from the calving front.The comparison between the viscous and viscoelastic solutions is shown in the bottom row of both figures (Fig.4.4 and Fig. 4.5). The difference in vertical velocity is the strongest near the calving front, and shows a damped oscillating profile with distance to the grounding line. The oscillating velocity difference centres around zero and is equal over shelf thickness, which suggests that this is the difference induced by the elastic part. The horizontal velocity also shows a damped oscillating difference in velocity, but smaller in amplitude (maximum range 5.7 m a -1 ) compared to the vertical difference (maximum range 12.7 m a -1 ). The horizontal velocity difference shows a negative trend towards the calving front, suggesting an overestimation by the viscous model by ∼ 1 -2 m a -1 . Interestingly, the velocity difference between the top and bottom surface have opposing phases. This can be related to the bending moment considered in the viscoelastic model, which causes

Figure 4 . 4 :

 44 Figure 4.4: Model comparison by velocity component difference (2D profile) between two models per graph, following Fig. 4.2. Shown are the differences of the horizontal (∆u x -left) and vertical (∆u y -right) velocity components per model comparison.

Figure 4 . 5 :

 45 Figure 4.5: Model comparison by velocity component difference (graph) between two models per graph, following Fig. 4.3. Shown are the differences of the horizontal (∆u x -left) and vertical (∆u y -right) velocity components per model comparison at the 'top surface', at 'sea water level' and at the 'bottom boundary'.

  chapter, and kept the model parameters constant. The material properties used here are: the viscosity η = 10 14 Pa s and the elastic modulus G = 10 9 Pa representing a median value from Tab. 2.1. The densities were set to ρ i = 910 kg m -3 for ice and ρ w = 1028 kg m -3 for water and gravity g = 9.81 m s -2 were defined with widely used values. The boundary conditions as described in subsection 3.3.2 were used with the Robin boundary condition on the water boundary (eq. 3.47 + 3.50).

A

  comparison of the velocity profiles shows a close to linear increase of the horizontal velocity component (u x -red in Fig. 4.6) over the length of the ice shelf, starting at 0 km a -1 on the domain inflow (Dirichlet boundary condition). The vertical velocity component (u y -blue in Fig. 4.6) is constant with an oscillation in the frontal ∼ 2500 m. This oscillation can be attributed to the bending moment caused by the imbalance of hydrostatic and lithostatic pressure at the calving front (Fig. 2.4).

Figure 4 . 6 :

 46 Figure 4.6: Sensitivity to the mesh resolution: Evolution of the horizontal (u x -red) and vertical (u y -blue) velocity components along an horizontal line at sea water level for different mesh densities MD. See Tab. 4.1 for details.

Fig. 4 .

 4 Fig. 4.8 shows the evolution at sea water level of the horizontal (u x -red) and vertical (u y -blue) velocity over the length of a 10 km ice shelf.

Figure 4 . 8 :

 48 Figure 4.8: Sensitivity to ice shelf thickness at the calving front: Evolution of the horizontal (u x -red) and vertical (u y -blue) velocity components along an horizontal line at sea water level for different heights at calving front (H cf ). See Tab. 4.3 for details.

  maximum velocities are recorded at the calving front and show a linear increase of velocity with calving front thickness, where the ice ramp with H cf = 100 m experiences the slowest velocities and the shelf with H cf = H gl = 400 m the fastest velocities. Here, the vertical velocity has a stronger increase than the horizontal velocity, proportionally.

Figure 4 . 9 :

 49 Figure 4.9: The maximum absolute velocity per ice shelf thickness at the calving front, presented for the vertical velocity component (u x -red) and horizontal velocity component (u y -blue). The shelf length and thickness at grounding line was kept constant at 10 km and 400 m respectively.

  The horizontal velocity component (u x -red in Fig.4.10) increases with distance to the grounding line starting from the prescribed value of 0 km a -1 and showing a continuous damping of the velocity increase over the whole length, as shown previously (subsection 4.2.2). Additionally, a steeper increase at the last ∼ 1000 m near the calving front can be observed for each dataset. This can be related to the bending moment at the calving front, which causes the upper frontal part of the self front to move faster.

Figure 4 .

 4 Figure 4.10: Sensitivity to ice shelf length: Evolution of the horizontal (u x -red) and vertical (u y -blue) velocity components along an horizontal line at sea water level for different lengths. See Tab. 4.4 for details.

Fig. 4 .

 4 Fig. 4.11 shows the Maxwell time computed by eq. 2.2 for all investigated shear modulus (G) and viscosity (η) parameter combinations, with the Maxwell time ranging from 0.01 s to ∼ 317 years. For a representative shear modulus and viscosity for ice of G = ∼ 1.8 × 10 9 Pa and η = 10 14 Pa s respectively, we get a Maxwell time of λ = ∼ 15 hours.

Figure 4 .

 4 Figure 4.11: Maxwell time and model output for constant forcing. The table shows the Maxwell time for all investigated shear modulus (G) and viscosity (η) combinations. The color scheme indicates the successfulness of the model run, where each cell represents 2 model setups with varying parameter values for height at calving front (H cf = 200/400 m).In case of the grey cells, no successful model run was recorded, while parameter combinations of the white cells indicate some successful model results. For green cells all parameter combinations produced result with constant forcing.

Figure 4 .

 4 Figure 4.12: Velocity components (u x ) and (u y ) as function of the shear modulus (G) and viscosity (η) on the calving front at sea water level. Each point represents one successful model solution from Fig. 4.11 and the colours represent the 2 geometries with different height at calving front (H cf = 200/400 m).

Fig. 4 .

 4 Fig. 4.13 and Fig. 4.14 present the solution velocity profiles for the vertical component (u y ) and horizontal component (u x ), respectively. The profiles for all investigated viscosities and geometries are shown, for a shear modulus of G = 10 12 Pa. The solutions are only presented for one shear modulus, because there are

Figure 4 .

 4 Figure 4.13: Vertical velocity profiles (u y ) for all investigated viscosities (η) with a shear modulus of G = 10 12 Pa. The velocity profiles shown are at different scales and tables with all velocity component profiles can be found in Appendix A.

Figure 4 . 1 :

 41 Figure 4.14: Horizontal velocity profiles (u x ) for all investigated viscosities (η) and a shear modulus of G = 10 12 Pa. The velocity profiles shown are at different scales and tables with all velocity component profiles can be found in Appendix A.

Fig. 5 .

 5 Fig. 5.1 shows the simulated tidal signals on the top row separated by tidal cycle, which were computed with eq. 5.1 and the introduced tidal parameters in Tab. 5.1. The bottom row shows the related vertical water level velocity ('wlv') as a function of time, computed by the first derivative of eq. 5.1. It shows that the peak vertical water velocities are shifted by one-quarters of a wavelength from the tidal peaks and the semi-diurnal cycle (blue) results in around double the velocities compared to the diurnal cycle (purple).

Figure 5 . 1 :

 51 Figure 5.1: Simulated tidal signals (top) and their corresponding vertical water level velocities (bottom) separated by the tidal cycle. Detail on the tidal parameters can be found in Tab. 5.1.

  5.1) on an ice shelf ramp of L = 10 km, H gl = 400 m, H cf = 200 m and a mesh density of 2. The remaining material parameters were set to: ρ i = 910 kg m -3 , ρ w = 1028 kg m -3 , g = 9.81 m s -2 , η = 10 13 Pa s and G = 10 9 Pa. The values for viscosity and shear modulus were chosen as close as possible to the representative values of ice as mentioned in literature (see subsection 2.3.3), but unfortunately, the model did not give results for all parameter combinations (see next section 5.3). The results are presented in Fig. 5.2 showing the velocity responses over time, and in Fig. 5.3 showing the spacial distribution of the maximum velocities. The velocity component results in Fig. 5.2.c-f show a strong wiggling in the very first tidal quarter, this is model ramp up behaviour where it adjusts to the forcing. This part should be generally ignored for any analysis. For the temporal statistics analysis shown in Fig. 5.3, we cut the first three-quarters of the tidal cycle, to ensure to get a clean signal without these ramp-up disturbances. The vertical velocity (u y ; Fig. 5.2.e + f) responds at the calving front ('(uswl)') with all 4 cases showing a clean sinusoidal curve that aligns well with the water level velocity (wlv; Fig. 5.1.c + d). Both the wlv and u y of each tidal configuration show the same phase shift -one-quarters of a wavelength -to the tidal signal, and hence have a near zero phase shift between one other, which indicates a purely elastic response. Additionally, the vertical velocity response and wlv show the same velocity amplitudes per tidal configuration. The shape of the vertical velocity range ('(u -range)') -the minimum and maximum velocities over

Figure 5 . 2 :

 52 Figure 5.2: Velocity response to different sinusoidal tidal forcings. The range of the velocity field are shown as well as the velocity at the intersection of sea water level equilibrium and the calving front for an ice shelf ramp (H cf = 200 m) with material parameters G = 10 9 Pa and η = 10 13 Pa s.

Figure 5 . 3 :

 53 Figure 5.3: Temporal statistics showing the max velocity response over several tidal cycles for the examples in Fig. 5.2. The first cycle has been removed for clean signal. The velocity profiles shown are at different scales.

Figure 5 . 4 :

 54 Figure 5.4: Evolution of the velocity distribution over one tidal cycle for the maximum amplitude case with a semi-diurnal tidal cycle (Fig. 5.3.f+h) on the 'ramp' geometry. The velocity profiles shown are at different scales and the timestamp of the snapshot are indicated in the tidal signal graph on top.

5. 3

 3 Fig. 5.5 shows, as in subsection 4.3.2 -Fig. 4.11, the Maxwell time and model success rate. But this time the figure shows the model success rate for the tidal forcing of 4 different configurations: a combination of the 2 tidal periods and the 2 heights at calving front. Dark grey cells did not have a steady state result (see Fig. 4.11), hence no model run with tidal forcing could be attempted. Light grey cells, did mostly have a steady state result, but the model did not produce a solution for any of the 4 configurations. I have tried different time-step configurations, but unfortunately none have worked. The white to dark green colours indicate the amount of successful model configurations, with white indicating one solution and dark green representing parameter combinations that always produced results. There is a visible diagonal transition from the bottom left to the top right in terms of model success, and only the furthest right column (G = 10 12 Pa) produces model solutions for all investigated viscosities. It also appears that the model was only successful for Maxwell times lower than 1 day, with exceptions to the lowest and highest investigated shear modulus. For values that are widely used for ice as discussed in subsection 2.3.3 (G = ∼ 1.8 × 10 9 Pa and η = 10 14 Pa s) unfortunately, only one solution was found.

Figure 5 . 5 :

 55 Figure 5.5: Maxwell time and model output for tidal forcing. The table shows the Maxwell time for all investigated shear modulus (G) and viscosity (η) combinations. The color scheme indicates the successfulness of the model run, where each cell represents 4 model setups with an tidal amplitude of 5 m and varying parameter values for: tidal periods (12.82/25.82 h) and height at calving front (H cf = 200/400 m). In case of the grey cells, no successful model run was recorded, while parameter combinations of the white to green cells indicate how many successful model results were recorded.

Fig. 5 .

 5 Fig. 5.6 shows the maximum velocity distribution over time in the ramp geometry for the horizontal component (u x -top) and vertical component (u y -bottom) and all investigated viscosities (η -rows) and a selection of shear modulus (Gcolumns).

Figure 5 . 6 :

 56 Figure 5.6: Temporal Statistics profiles of the horizontal (u x ) and vertical (u y ) velocity field depending on shear modulus (G) and viscosity (η) -case: semi-diurnal, H cf = 200 m. The velocity profiles shown are at different scales. Only a selection of the investigated shear modulus are shown, tables with temporal statistics profiles for all parameter combination can be found in Appendix B.

Figure 5 . 7 :

 57 Figure 5.7: Horizontal velocity profiles (u x ) for all investigated viscosities (η) and a shear modulus of G = 10 12 Pa -case: H cf = 200 m ('ramp'). The velocities at sea water level at different distances to the grounding line: 0 km (orange) to 10 km (dark red) are shown. Tables with velocity component profiles for all parameter combination can be found in Appendix C.

Figure 5 . 8 :

 58 Figure 5.8: Horizontal velocity profiles (u x ) for all investigated viscosities (η) and a shear modulus of G = 10 12 Pa -case: H cf = 400 m ('rectangle'). The velocities at sea water level at different distances to the grounding line: 0 km (orange) to 10 km (dark red) are shown. Tables with velocity component profiles for all parameter combination can be found in Appendix C.

Figure 5 . 9 :

 59 Figure 5.9: Vertical velocity profiles (u y ) for all investigated viscosities (η) and a shear modulus of G = 10 12 Pa -case: H cf = 200 m ('ramp'). The velocities at sea water level at different distances to the grounding line: 0 km (light blue) to 10 km (dark blue) are shown. Tables with velocity component profiles for all parameter combination can be found in Appendix C.

Figure 5 .

 5 Figure 5.10: Vertical velocity profiles (u y ) for all investigated viscosities (η) and a shear modulus of G = 10 12 Pa -case: H cf = 200 m ('rectangle'). The velocities at sea water level at different distances to the grounding line: 0 km (light blue) to 10 km (dark blue) are shown. Tables with velocity component profiles for all parameter combination can be found in Appendix C.

Fig. 5 .

 5 Fig. 5.11 compares the total shear stress (τ [Pa]) to the viscous shear stress (2ηD(u) [Pa]) for 3 of the investigated viscosities (η = 10 10 , 10 14 , 10 18 Pa s), a shear modulus of G = 10 12 Pa, the 'ramp' geometry and the tidal configuration: maximum amplitude and 'semi-diurnal' cycle. The results are shown over time and are separated by normal stress components (horizontal -orange/red; vertical -blue) at several locations along the sea water level at equilibrium line. The locations picked are: 0 km, 2 km, 4 km, 6 km, 8 km and 10 km from the grounding line, where the most faded colour represents the stress component over time at the grounding line and the darkest colour the stress component at the calving front. Since the stresses vary cyclically over time with the tide, they form a circle. If the total

Figure 5 .

 5 Figure 5.11: Shear stress comparison for 3 of the investigated viscosities (η = 10 10 , 10 14 , 10 18 Pa s) and a shear modulus of G = 10 12 Pa -case: H cf = 200 m, maximum amplitude, semi-diurnal. Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx -orange/red) and vertical ( yy -blue) component at different distances to the grounding line: 0 km (light) to 10 km (dark). Tables with the stress comparison graphs for all parameter combination can be found in Appendix D.

  3, I initiated an investigation over a much longer time of over 4.75 years. A 'rectangle' shelf geometry with the specification: L = 10 km, H cf = H gl = 400 m and MD=2 was used with a diurnal maximum tidal amplitude ('Max (d)') case and material parameters of G = 10 9 Pa and η = 10 12 Pa s.

Fig. 5 .

 5 Fig.5.12 shows the resulting horizontal and vertical velocity components over the whole time span, while Fig.5.13 is a detailed view of the first 105 days.

Figure 5 .

 5 Figure 5.12: Long time study of the linear viscoelastic model, performed over 4.75 years for a viscosity of η = 10 12 Pa s and a shear modulus of G = 10 12 Pa with the diurnal maximum tidal amplitude ('Max (d)') case on a 'rectangle' geometry. The velocities are taken at sea water level on the calving front .

Figure 5 .

 5 Figure 5.13: Long time study of the linear viscoelastic model -detail Fig.5.12 zoom into the first 105 days.

  2 I introduced a simplified daily tidal forcing to the model and compared 4 different tidal configurations. Two realistic tidal forcings based on two locations in Antarctica [Padman et al. 2018], the Weddell Sea with the diurnal dominant frequency and the Ross Sea with the semi-diurnal frequency. Since both locations had very different spring tide amplitudes (1.5 m and 0.7 m, respectively), I chose to compute both tidal frequencies with the same arbitrary tidal amplitude of 5 m to allow for a tidal cycle comparison. The value of 5 m was picked near the maximum possible tidal amplitude found on earth, to enhance the influence this has on the ice shelf.

  3). The shear modulus appeared to have no noticeable effect on the solution and only impacted the success rate of the model. Therefore, a study of the influence of the Young modulus and Poisson ratio was not considered to bring any new insights, as the Poisson ratio only reduces the Young modulus by 1 -0.5 orders of magnitude, for the simplified shear modulus used, and I investigated the shear modulus over a range of several orders of magnitude. It raises the question if the neglected bulk modulus would change the impact of the elasticity on the solution. It appears that the model was mainly successful for shear moduli larger than 3 × 10 8 Pa and under tidal forcing additionally for Maxwell times lower than 1 day, with small exceptions.

2 . 5 )

 25 from the Rutford ice stream (Weddell sea, semi-diurnal) as well as GPS data from the Bindschadler ice stream (Ross sea, diurnal cycle). The measurements used on Rutford ice stream are 20 km downstream of the grounding line and on Bindschadler ice stream are 15 km downstream, but as being part of large ice shelves, both measurements are very far from the calving front. They compare the tidal displacement with the horizontal displacement of the shelf, and found the displacement on the Rutford ice stream is in phase (0 °) with the tide, while at Bindschadler ice stream the horizontal displacement is out of phase (180 °).[START_REF] Podrasky | Quantifying velocity response to ocean tides and calving near the terminus of Jakobshavn Isbrae, Greenland[END_REF] identified that the horizontal displacement just upstream of Jakobshaven Isbrae is 90 °out of phase to the semi-diurnal tidal signal. These measurements are difficult to compare, as they are in the middel of an ice shelf or near a terminus of a glacier with different tidal cycles at different places in the world, while the simulations done in this study contain unrealistic assumptions, like the fixed grounding line, but it highlights that in nature as well as in this study we find varying phase shifts between the tidal signal and the ice's response. Further possible explanations for the phase shift are either in relation to the viscosity or a introduced delay from the Robin boundary condition, which calculates the displacement correction based on the previous time steps velocity. Or maybe a combination of both. For a better understanding of the phase shift and the related classification of the flow regime, further investigations are needed.

D. 3

 3 VShear stress comparison at tidal forcing for shelf height at calving front H cf = 200 m -Weddell Sea example . . . . . . . . . . . . . . . . D.4 Shear stress comparison at tidal forcing for shelf height at calving front H cf = 200 m -Maximum amplitude, semi-diurnal case . . . . . D.5 Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Ross Sea example . . . . . . . . . . . . . . . . . . D.6 Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, diurnal case . . . . . . . . D.7 Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Weddell Sea example . . . . . . . . . . . . . . . . D.8 Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, semi-diurnal case . . . . .

Figure A. 1 :

 1 Figure A.1: Velocity fields (u [m a -1 ]) at constant forcing for shelf height at calving front H cf = 200 m . Shown are the horizontal (u x -top) and vertical (u y -bottom) velocity components after 6000 model iterations for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2. The velocity fields shown are at different scales.

Figure A. 2 :

 2 Figure A.2: Velocity fields (u [m a -1 ]) at constant forcing for shelf height at calving front H cf = 400 m . Shown are the horizontal (u x -top) and vertical (u y -bottom) velocity components after 6000 model iterations for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2. The velocity fields shown are at different scales.

Figure B. 1 :

 1 Figure B.1: Velocity fields (u max [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 200 m -Ross Sea example . Shown are the maximum velocity for its horizontal (u x -top) and vertical (u y -bottom) components derived from temporal statistics ( where the first 1-2 tidal cycles were removed) for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2. The velocity fields shown are at different scales.

Figure B. 2 :

 2 Figure B.2: Velocity fields (u max [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 200 m -Maximum amplitude, diurnal case . Shown are the maximum velocity for its horizontal (u x -top) and vertical (u y -bottom) components derived from temporal statistics ( where the first 1-2 tidal cycles were removed) for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2. The velocity fields shown are at different scales.

Figure B. 6 :

 6 Figure B.6: Velocity fields (u max [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, diurnal case . Shown are the maximum velocity for its horizontal (u x -top) and vertical (u y -bottom) components derived from temporal statistics ( where the first 1-2 tidal cycles were removed) for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2. The velocity fields shown are at different scales.

Figure C. 1 :

 1 Figure C.1: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 200 m -Ross Sea example . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 2 :

 2 Figure C.2: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 200 m -Maximum amplitude, diurnal case . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 3 :

 3 Figure C.3: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 200 m -Weddell Sea example . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 4 :

 4 Figure C.4: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 200 m -Maximum amplitude, semi-diurnal case . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 5 :

 5 Figure C.5: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 400 m -Ross Sea example . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 6 :

 6 Figure C.6: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, diurnal case . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 7 :

 7 Figure C.7: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 400 m -Weddell Sea example . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure C. 8 :

 8 Figure C.8: Velocity profiles (u [m a -1 ]) at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, semi-diurnal case . Shown are the velocity profiles over time for its horizontal (u x -top) and vertical (u y -bottom) components at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 1 :

 1 Figure D.1: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 200 m -Ross Sea example . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 2 :

 2 Figure D.2: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 200 m -Maximum amplitude, diurnal case . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 3 :

 3 Figure D.3: VShear stress comparison at tidal forcing for shelf height at calving front H cf = 200 m -Weddell Sea example . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 4 :

 4 Figure D.4: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 200 m -Maximum amplitude, semi-diurnal case . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 5 :

 5 Figure D.5: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Ross Sea example . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 6 :

 6 Figure D.6: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, diurnal case . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 7 :

 7 Figure D.7: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Weddell Sea example . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Figure D. 8 :

 8 Figure D.8: Shear stress comparison at tidal forcing for shelf height at calving front H cf = 400 m -Maximum amplitude, semi-diurnal case . Shown are the total shear stress (τ [Pa]) over the viscous shear stress (2ηD(u) [Pa]) over time for its normal horizontal ( xx ) and vertical ( yy ) component at the following locations along sea water level (measured in distance from the grounding line): 0 km (grounding line), 2 km, 4 km, 6 km, 8 km and 10 km (calving front); for the varying input parameters: shear modulus (G [Pa]) and viscosity (η [Pa s]). The remaining geometry parameters are: shelf length L = 10 km, height at grounding line H gl = 400 m and a mesh density of 2.

Table of Content

 of Front matter Declaration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Résumé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v

	Abbreviations & Acronyms Symbols & Notations
	div	AASW divergence operator Antarctic Surface Water
	grad	AIS gradient operator Antarctic Ice Sheet
		AP	Antarctic Peninsula
	u	BAS velocity vector British Antarctic Survey	m s -1
	BC D(u) deformation rate tensor (symmetric part of the velocity gradient) s -1 Boundary Condition
	CDW W (u) vorticity tensor (asymmetric part of the velocity gradient) Circumpolar Deep Water	s -1
	σ	DBC stress tensor	Dirichlet Boundary Condition	P a
	τ	DInSAR differential satellite synthetic aperture radar interferometry deviatoric stress tensor (dev(σ)) P a
	ϵ	EAIS strain tensor	East Antarctic Ice Sheet	-
	ε	FEM deviatoric strain tensor (dev(ϵ)) Finite Element Method	-
	ε	GL strain rate tensor Grounding Line	s -1
	ε	GPS deviatoric strain rate tensor (dev( ε)) Global Positioning System	s -1
	ρ i	ICESat density of ice	Ice, Cloud and land Elevation Satellite	kg m -3
	ρ w	IG density of (sea)water Intragravity waves	kg m -3
	E	IPCC Young's modulus Intergovernmental Panel on Climate Change	P a
	G	LIMA shear modulus Landsat Image Mosaic if Antarctica	P a
	η	MD viscosity	Mesh Density	P a s
	ν	MISI Poisson's ratio Marine Ice Sheet Instability	-
		NBC	Neumann Boundary Condition
		PDE	Partial Differential Equation
		RBC	Robin Boundary Condition
		SAR	Synthetic Aperture Radar
		SSA	Shallow Shelf Approximation
		VLP	very long period gravity waves
		WAIS	West Antarctic Ice Sheet
			ix

x Chapter 1

  West Antarctica, separated from East Antarctica by the Transantarctic Mountains, contains Antarctica's highest mountain: Mount Vinson, measuring 4892 m above sea water level and Antarctica's lowest bed elevation: the Denman canyon, reaching down to 3500 m below sea level. East Antarctica, being a sig-

	nificantly larger area than West Antarctica, is covered by the thickest land ice,
	measuring up to 4776 m [BAS 2005; Morlighem et al. 2020].
	The Antarctic Ice Sheet (AIS), grounded land ice, which is shown as white in
	Fig. 2.1, covers approximately 98 % of Antarctica and is often split into West
	Antarctic Ice Sheet (WAIS) and East Antarctic Ice Sheet (EAIS). While the
	bedrock EAIS is resting on is mainly above sea level, the WAIS bedrock-ice in-

2.1)

. The Antarctic Peninsula is a narrow stretch of land, covered with mountain ranges, pointing towards the southern tip of South America.

Table 2 . 1 :

 21 Elastic properties of ice, specifically the instantaneous Young modulus (E) and Poisson's ratio (ν). The values were either measured in a 'lab' (with 1 ultrasonic pulse technique or 2 sound waves) or an elastic ('e') or viscoelastic ('ve') model was fitted

		ν [-]	method	Reference
	0.001 -0.01 -	e (tp ; -; -)	Wagner et al. [2016]
	0.009 -12	0.33 -(possibly lab) Schulson [1999]
	0.09 -0.9	0.3	e (tp ; -; -)	Le Meur et al. [2014]
	0.88 ± 0.35	-	e (tp ; -; -)	Vaughan [1995]
	1.62 ± 0.69	0.4	ve (tp ; M ; l)	Wild et al. [2017]
	3.0 ± 0.2	0.3	e (tp ; -; -)	Rignot [1996]
	3.0 -10.0	0.35 lab 1	Gold [1977]
	3.2	0.3	ve (fs ; M ; nl) Rosier et al. [2017]
	4.8	0.41 ve (fs ; M ; nl) Gudmundsson [2011]
	9.3	-	ve (tp ; B ; l)	Reeh et al. [2003]
	9.3	-	lab 2	Petrenko and Whitworth [1999]
	10	0.3	ve (tp ; M ; nl) MacAyeal et al. [2015]

3 -0.4

are mostly assumed, see

Tab. 2.1. E [GPa]

to in situ data. The models used either the thin-beam/plate ('tp') or full stokes ('fs') approach and the viscoelastic model was either of Maxwell ('M') or Burgers ('B') type with a linear ('l') stress-strain-rate relationship or non-linear ('nl').

  3.25) in the viscoelastic Stokes equation (eq. 3.31). Here, the problem has to be computed over time and the solution varies at each time step, the so-called transient state. In case the input parameters

do not change over time (constant forcing) the solution approaches a near constant state with time. This (near) steady state is reached if the difference between the solution of two succeeding time steps tends to zero, or in other words if the residual term converges.

Table 4 . 1 :

 41 km, a thickness at grounding line (H gl ) of 400 m and at calving front (H cf ) of 200 m was simulated Mesh configuration for a ramp geometry with varying mesh densities (MD), where α is the surface slope, l gl is the element length at grounding line and l cf is the element length at the calving front.

	with mesh densities (MD) of 1, 2, 5, 7 and 10 elements per shelf thickness, see
	Tab. 4.1 for mesh details.
						Element length	Mesh
	L	H gl	H cf	α	MD	l gl	l cf	nodes elements
					1	400 m	200 m	56	72
					2	200 m	100 m	111	146
	5 km 400 m 200 m ∼ 0.0046 rad	5	80 m	40 m	688	1190
					7 ∼ 57.1 m ∼ 28.6 m	1323	2386
					10	40 m	20 m	2395	4420

Table 4 . 3

 43 4.2.2 Sensitivity to ice shelf thickness at calving frontAn ice shelf can be several kilometres long but usually is only several hundreds of meters thick, resulting in a difference between shelf length and thickness of 1 to 3 orders of magnitude. For this reason, ice shelves are sometimes modelled as a rectangle, where the shelf thickness at the inflow of the domain (in this case H gl ) and calving front (H cf ) are the same.In this section, different shelf thicknesses at the calving front are investigated for a constant shelf length L = 10 km, thickness at grounding line H gl = 400 m and mesh density MD = 2. The shelf thicknesses at calving front range from 100 m to 400 m, see Tab. 4.3 for more details. The later representing the shelf as a rectangle shape (H gl = H cf ).

	Element length	Mesh

: Mesh configuration for a ramp geometry with varying heights at calving front (H cf ), where α is the surface slope, l gl is the element length at grounding line and l cf is the element length at the calving front.

Table 4 . 4 :

 44 Mesh configuration for a ramp geometry with varying lengths (L), and MD= 2, where α is the surface slope, l gl is the element length at grounding line and l cf is the element length at the calving front.

	Element length	Mesh

Acknowledgements Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Table of Content . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Abbrevations & Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix

Viscoelastic model with constant forcing

In this chapter the previously introduced viscoelastic model is investigated under constant forcing. At first the viscoelastic model is compared to a viscous model as well as the analytical SSA solution with the aim to highlight the differences between the models. Further, a sensitivity study of the viscoelastic model with regards to geometry dimensions, mesh resolution and material parameters is presented.

Model comparison

In this section I compare the analytical linear Shallow Shelf Approximation ('SSA' -subsection 2.3.4), the linear viscous Stokes model ('v') and the introduced linear viscoelastic Maxwell model ('ve' -section 3.3).

The viscous model for ice is mathematically described by the Full-Stokes equation (see eq. 3.28). For the implementation in Rheolef (see subsection 3.2.2) I have adapted 'The Stokes problem' example described in Saramito [2020, ch. 2.1.4] for ice with the same boundary conditions (Robin) as used for the viscoelastic model (see subsection 3.3.2).

The Shallow Shelf Approximation was computed in Paraview on the loaded geometry, using the equations described in Greve and Blatter [2009, eq. 6.77, eq. 6.81]. The boundary conditions were defined as described in subsection 3.3.2, where size (number of nodes). The difference of the MD5 and MD7 to the finest mesh (MD10) is ∼ 0.1 m a -1 or less. The MD2 shows a minimal differences at the calving front which represents an error compared to the finest mesh of ∼ 1 m a -1 , while the MD1 shows the strongest differences (∼ 5 m a -1 ). The error to the finest mesh, presented in Fig. 4.7, decreases with the number of nodes (increase of mesh density Tab. 4.1), following a fit between 1 st and 2 nd order convergence as expected by the code, shown as dotted and dashed line respectively. Tab. 4.2 presents the duration of each model run, performed on the same PC under the same overall workload. While the two smaller mesh densities were quite similar and under 10 minutes, the MD10 took over 2.5 hours. A strong increase in computational time can be seen with increasing MD. This is related to the number of nodes per mesh, which defines the stiffness matrix and is the basis for the computation. Following these results and depending on the size of the phenomena to be investigated, e.g. melt ponds, a mesh density of 2 would be reasonable, because it keeps the computational cost low. For more accuracy a mesh with a higher density might be chosen. For the reminder of this thesis a MD2 will be used.

Model sensitivity to material parameters under constant forcing

Following the previous section 4.2, this section investigates the effect of the material parameters: shear modulus (G) and viscosity (η) on the model solution under constant forcing.

Parameter study configuration

The viscosity (η) of ice ranges between 10 13 -10 17 Pa s, with a value of 10 13 -10 15 Pa s being most reasonable (see subsection 2.

Summary

In this section we studied the impact of the material parameters: shear modulus (G) and viscosity (η) on the linear viscoelastic model outcome.

The shear modulus (G) appears to have an impact on the general success rate (convergence) of the model (Fig. 4.11), but no actual effect on the solution, here velocity (Fig. 4.12.a + c). It was shown that the model did not or only partially find a solution for lower shear modulus (G < 8 × 10 8 Pa), which also depended on the shape of the geometry, as we were able to produce more results for lower shear modulus with the 'rectangle' geometry than with the 'ramp' geometry. A study of the actual, underlaying material parameters of the shear modulus: Young modulus (E) and Poisson's ratio (ν) were dismissed, as no new insights were expected.

The viscosity (η) seems to have only a minor, secondary impact on the model success rate (Fig. 4.11), but largely affects the model solutions. The results have

shown that the viscosity is the main driver of the magnitude of the horizontal as well as vertical velocity (Fig. 4.12.b + d), with larger velocities at lower viscosities, which is to be expected since the viscosity is the resistance to deformation of a fluid. This also appears to drive the impact of the bending moment induced by the buoyancy force at the calving front, which showed a larger oscillation wavelength with larger viscosity.

Chapter 5

Viscoelastic model with tidal forcing

This chapter continues the investigation of the linear viscoelastic model by applying a tidal forcing. After introducing the simplified tidal change pattern for sea water level used by this model, we investigate the model's sensitivity to tidal parameters and the tidal impact under different material parameters.

Simplified tidal study configuration

Tides are an oscillation of sea water level, as already introduced in subsection 2.2.5. These oscillations show a daily pattern, low and high tide, and a monthly pattern, variation of the tidal range. In this study we only look at the short term, daily tidal variation and its impact on the ice shelf.

The daily tidal oscillations are mathematically prescribed by a sinusoidal signal:

where Amp is the maximum tidal amplitude, f req the tidal frequency defined as the reciprocal of the tidal period and t is the time. The tidal signal was implemented in the model so that it is computed before each iteration and added to the sea water level at equilibrium. ). The geometry is a 10 km long shelf with a height at grounding line (H gl ) of 400 m and a mesh density of 2.

• Shelf height at calving front

The following pages are not suitable for viewing on A4 print, best viewed in digital.

Appendix B

2D velocity field -tidal forcing

The following graphs show the maximum velocity (u [m a -1 ]) for its horizontal (u The following pages are not suitable for viewing on A4 print, best viewed in digital.

Appendix C

Velocity profiles along sea water level -tidal forcing

The following graphs show the velocity profiles (u [m a -1 ]) over time for its hori- The following pages are not suitable for viewing on A4 print, best viewed in digital.

Appendix D

Shear stress comparison -tidal forcing

The The following pages are not suitable for viewing on A4 print, best viewed in digital.