Dynamics of the Antarctic ice sheet withcoupled ice shelves
南極氷床変動における棚氷の影響
Tatsuru Sato
PhD Dissertation
Division of Earth System Science
Graduate School of Environmental Science
Hokkaido University

May 15, 2012

Contents
Abstractiii

Acknowledgementsv

1 Introduction
1.1 The Antarctic ice sheet . . . . . . . . . . . . . . . . . . . . . .
1.2 Ice shelves . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Global temperature and sea level change . . . . . . . . . . . .
1.4 Changes in the Antarctic ice sheet and ice shelves . . . . . . .
1.5 Ice sheet model SICOPOLIS . . . . . . . . . . . . . . . . . . .
1.6 Purpose and tasks of the study . . . . . . . . . . . . . . . . .

1
1
2
4
6
8
9

2 Theory
2.1 Dynamics of large ice bodies . . . . . . . . . . . . . . . . . . .
2.2 Constitutive equations . . . . . . . . . . . . . . . . . . . . . .
2.3 Stress balance equations in general coordinates . . . . . . . . .
2.4 Large-scale ice sheet dynamics . . . . . . . . . . . . . . . . . .
2.5 Hydrostatic approximation . . . . . . . . . . . . . . . . . . . .
2.6 Shallow shelf approximation . . . . . . . . . . . . . . . . . . .
2.6.1 Formulation in orthogonal curvilinear coordinates . . .
2.6.2 Boundary conditions at the ice surface . . . . . . . . .
2.6.3 Boundary conditions at the ice base . . . . . . . . . . .
2.6.4 Vertical integration . . . . . . . . . . . . . . . . . . . .
2.7 Shallow ice approximation . . . . . . . . . . . . . . . . . . . .
2.8 Temperature evolution of ice . . . . . . . . . . . . . . . . . . .
2.9 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . .
2.9.1 Basal sliding . . . . . . . . . . . . . . . . . . . . . . . .
2.9.2 Calving front . . . . . . . . . . . . . . . . . . . . . . .
2.9.3 Surface temperature and mass balance . . . . . . . . .
2.9.4 Basal temperature and sub-marine mass balance . . . .
2.10 Topographic evolution . . . . . . . . . . . . . . . . . . . . . .

24
24
24
26
28
28
31
31
31
32
32
35
36
36
36
37
38
39
39

i

2.10.1 Ice thickness equation . . . . . . . . . . . . . . . . . . 39
2.10.2 Ice surface equation . . . . . . . . . . . . . . . . . . . . 40
3 Numerical methods
3.1 Ice sheet model SICOPOLIS . . . . . . . . . . . . . . . . . . .
3.2 Details of the numerical grid . . . . . . . . . . . . . . . . . . .
3.3 Coupling of the ice shelf model to the ice sheet model . . . . .
3.4 Shallow shelf dynamic equations . . . . . . . . . . . . . . . . .
3.4.1 Field equations . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Calving front conditions . . . . . . . . . . . . . . . . .
3.4.3 Grounding zone . . . . . . . . . . . . . . . . . . . . . .
3.5 Temperature evolution . . . . . . . . . . . . . . . . . . . . . .
3.6 Topography evolution . . . . . . . . . . . . . . . . . . . . . . .
3.6.1 Ice surface equation . . . . . . . . . . . . . . . . . . . .
3.6.2 Ice thickness equation . . . . . . . . . . . . . . . . . .

46
46
46
48
48
48
49
50
52
52
52
53

4 Numerical tests
57
4.1 Ice shelf ramp . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.2 Shelfy stream experiments . . . . . . . . . . . . . . . . . . . . 58
5 SeaRISE experiments
5.1 SeaRISE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Application of the model to the SeaRISE set-up . . . . . . . .
5.3 Paleoclimatic spin-up . . . . . . . . . . . . . . . . . . . . . . .
5.4 Sensitivities of the ice volume to different datasets . . . . . . .
5.5 Future climate experiments . . . . . . . . . . . . . . . . . . .
5.6 Sensitivities on varied sub-ice-shelf melting rates . . . . . . . .
5.7 Results of whole set of SeaRISE experiments . . . . . . . . . .
5.8 Relaxation of the initial shock . . . . . . . . . . . . . . . . . .

63
63
63
64
65
66
67
69
70

6 Discussion
6.1 Comparison to the present-day ice sheet . . . . . . . . . . . .
6.2 The dynamics of ice sheet and ice shelves . . . . . . . . . . . .
6.3 Grounding line migration . . . . . . . . . . . . . . . . . . . . .
6.4 Calving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90
90
91
92
94

7 Summary

101

Bibliography

103

ii

Abstract
Antarctica is Earth’s southernmost continent, situated almost entirely southof the Antarctic Circle. Except for some coastal regions and high mountains,most of the continent is covered by ice that constitutes the Antarctic ice
sheet. The Antarctic ice sheet is the largest land ice mass on the present-day
Earth. According to the Fourth Assessment Report of the United Nations
Intergovernmental Panel on Climate Change (IPCC AR4), the total volumeof the ice sheet amounts to approximately 25×106 km3 or 56.6 m SLE (metersof sea level equivalent).
Ice shelves are floating ice bodies that are attached to the margins of anice sheet. Half of the coastline of the whole Antarctic ice sheet and almost
the entire West Antarctic ice sheet are surrounded by ice shelves. Most of theice flow towards the coast is discharged through them. Recent observations
have led to strong concerns that ice-dynamical processes (loss of buttressing
from ice shelves, speed-up of ice streams and outlet glaciers) may boost themass loss and thus lead to an additional contribution to sea level rise. The
importance of understanding these processes to assess future sea level change
was explicitly recognized in the IPCC AR4.
SICOPOLIS is a dynamic/thermodynamic ice sheet model which has been
applied to large ice sheets such as the Greenland, Antarctic, Laurentide andFennoscandian ice sheets. Although it was applied to the Antarctic ice sheet,special treatment of ice shelf dynamics has not been included so far. In thisstudy, ice shelf dynamics is considered and coupled with the ice sheet model.
The model is applied to the Antarctic ice sheet to study instabilities of theice sheet related to ice sheet/ice shelf interactions.
Shallow shelf approximation is applied for the ice shelf model. The performance of the model is tested by using an analytical solution for a twodimensional ice shelf ramp. The result shows that the model can reproduce
the analytical solution with high accuracy. The ice shelf model is then coupled to the ice sheet model SICOPOLIS and applied to the Antarctic ice
sheet. Input data are adopted from the SeaRISE project, an international
community effort to estimate the contribution of ice sheets to future sea leveliii

rise in the next 100s of years. In order to obtain a suitable present-day configuration of the Antarctic ice sheet, it is desirable to carry out a paleoclimatic
spin-up over at least a full glacial cycle until the present. The paleoclimatic
spin-up is conducted with a fixed-topography approach (the geometry is kept
fixed over time, thus not allowed to evolve), while the model is forced by pasttemperature and accumulation changes estimated by proxies from the Vostok
ice core record.
Future climate experiments over 500 years are conducted by using theresult of the paleoclimatic spin-up as initial condition. Three types of sensitivity experiments (surface climate, sub-ice-shelf melt rate and basal sliding)are conducted to investigate the key elements for the evolution of the Antarctic ice sheet in a warming climate. Results show that the grounded ice sheet
volume reacts most strongly to changes in sub-ice-shelf melt rates. It also
decreases with increased basal sliding, although the range is not as large as
for sub-ice-shelf melt. The smallest influence is that of direct climate forcing
by surface temperature and precipitation. Combination of all the components causes the largest decrease of the ice sheet. The range of grounded
ice volume change across all experiments is ∼1.3 m SLE after 100 years and∼3.3 m SLE after 500 years. These results show that the contribution ofchanges of the Antarctic ice sheet to global sea level during the next 100
years can be larger than projected by the IPCC AR4.

iv

Acknowledgements
I would like to express my gratitude to my supervisor Ralf Greve who supported my study during the doctor course. I also thank Takeo Hondoh,
Takayuki Shiraiwa, Shin Sugiyama and Yoshinori Iizuka, who helped me not
only during the doctor course, but also during the previous master course.
Nina Kirchner and her colleagues supported me during my visit at Stockholm
University. Hakime Seddik helped me to use the Linux server in the laboratory. I thank Thorben Dunse and Fuyuki Saito for many fruitful discussions.
Thomas Gölles helped me to correct the thesis. I appreciate David Holland’s
invitation to the NYU Abu Dhabi workshop in March 2012. I would liketo thank the students and former students in our group: Toshimitsu Sakurai, Shun Tsutaki, Daisuke Nishimura, Akira Tsuchimoto, Takehiro Fukuda,
Takeshi Yoshizawa, Kenta Tone, Ikumi Oyabu, Daiki Sakakibara, Mikako
Nakata, Ayaka Okamoto, Satoshi Matsuno and Masahiro Minowa. Language
editing was done by Andrew Davis of English Experience Language Services(Jena, Germany). Ayako Abe, Keiichiro Ohshima and Takayuki Shiraiwa
kindly agreed to be examiners. Finally, I owe a big thanks to my family.
This study was supported by a Grant-in-Aid for Scientific Research A(No. 22244058) from the Japan Society for the Promotion of Science JSPS,
and by an exchange grant for a research visit from the Scandinavia-Japan
Sasakawa Foundation.

v

Chapter 1
Introduction
1.1

The Antarctic ice sheet

Antarctica is a continent lying in the Antarctic region. The climate is very
cold throughout the region; mean annual temperature is lower than the melting point of ice [Comiso, 2000]. Most of the continent is covered by ice, except
for some coastal regions and high mountains. This ice covering is called theAntarctic ice sheet (Fig. 1.1).
The Antarctic ice sheet is by far the largest land ice mass on Earth.
The total area of the ice sheet is 12.3 × 106 km2 . Its total ice volume is24.7 × 106 km3 [Solomon et al., 2007] and it contains about 90% of the ice
and 70% of the fresh water on Earth. It is approximately 56.6 m SLE (metersof sea level equivalent), i.e., it would produce a sea level rise of about 55 m if
it melted completely. Ice shelves are among the components of the ice sheet.
The Antarctic ice sheet is divided into two regions by the Transantarctic
Mountains (Fig. 1.1). One is the East Antarctic ice sheet (EAIS). Mostof the EAIS is located east of the Greenwich meridian. The other is theWest Antarctic ice sheet (WAIS). Most of the WAIS is located west of theGreenwich meridian. The EAIS is larger in area than the WAIS. The ice
there is also thicker and has a greater mass than that of the WAIS. Despite
its considerable thickness, the base of the EAIS lies mostly above sea level.
Many ice cores are drilled in the EAIS since it contains old ice, particularly
in inland regions.
Since most of the WAIS lies on bedrock below sea level, it is called amarine ice sheet. The bed lies sufficiently below sea level for ice shelves tobe formed often around the margins (Fig. 1.2). This condition may make the
WAIS more unstable than the EAIS, because a small retreat of its margins
can destabilize most of the WAIS. This then leads to rapid disintegration of1

large regions of ice.
The mass balance of the ice sheet is driven by accumulation and ice flow.
The Antarctic ice sheet is considered a desert in terms of accumulation rates.
The annual accumulation rate is hundreds of millimeters along the coast andtens of millimeters around the dome. Although accumulation rates are low,
the low velocity of the inland ice maintains its large volume.
Although the surface is cold, the base of the ice sheet is relatively warm
and melts in some regions due to geothermal heat and strain heating of ice.
Such melt water lubricates the ice sheet and allows basal sliding. This process
produces fast flowing channels in the ice sheet called ice streams. They are
the key component of the Antarctic ice sheet’s ice mass transport. Although
outlet glaciers and ice streams comprise only 13 of the Antarctic ice sheet
margins, they may drain as much as 90% of the accumulation in the interior[Morgan et al., 1982].

1.2

Ice shelves

The thickness of the ice sheet declines towards the coast. Where it becomesthin enough, it starts floating because the density of ice is less than that ofsea water. Ice shelves are thus formed. There are large areas of such floating
ice attached to glaciers or to the ice sheet (Fig. 1.3). The boundary between
the floating ice shelf and grounded ice, where the ice starts floating, is called
the grounding line or grounding zone. Ice is lost from the ice shelf to theocean by melting under the shelf and iceberg calving at the front of the ice
shelf. The boundary of the ice shelf and the ocean is called the calving front.
Ice shelves surround much of Antarctica. They are also formed in the
Arctic in some fjords and bays on the Greenland and Ellesmere Island coasts.
Alpine glaciers which terminate in a fjord can have a floating tongue which
is similar to an ice shelf. A total of 44% of the Antarctic coastline has ice
shelves attached. The total area is 1.5 × 106 km2 [Solomon et al., 2007].
Most of the WAIS is surrounded by ice shelves [Cuffey and Paterson, 2010].
Through these shelves 80% of the ice flow discharges towards the coastline.
The Ross Ice Shelf is the largest shelf with an area of 472, 960 km2 , whichis larger than Japan (377, 944 km2 ) and twice the area of Honshu island. The
shelf is about 800 km across. Its thickness varies from over 1000 m aroundthe grounding line to about 250 m at the calving front. Typical velocities
are a few hundred meters per year, increasing to more than 1 km/a at thefront [Rignot et al., 2011a].
The second largest is the Ronne-Filchner ice shelf. Its area is 422, 420 km2 ,which is almost as large as the Ross ice shelf.
2Other large ice shelves are Amery (62, 620 km2 ), Larsen C (48, 600 km2 ),Riiser-Larsen (48, 180 km2 ), Fimbul (41, 060 km2 ) and Shackleton (33, 820 km2 )(see Fig. 1.1). Most of the large ice shelves are located in West Antarctica.
Ice shelves are nourished by ice flow from the ice sheet, snow accumulatingon their surface, and ice accumulating on their base as sea water freezes
there. There is typically an input velocity of hundreds of metres per year
at the grounding region. It can be several thousand meters per year at thecalving front [Doake, 2001].
Shelves lose mass through ablation processes, iceberg calving and melting
at the base. The basal mass balance is controlled by the circulation of thesea below the ice shelf. The melt rate is often high near both the grounding
line and the ice front [Joughin and Padman, 2003]. Sea water sometimes
freezes in the central part and this is called marine ice. Grounding over anappreciable area produces a dome-shaped ice rise with its own radial flow
pattern. Less extensive grounding results in surface irregularities called ice
rumples, which move with the shelf.
The essential characteristics of an ice shelf are the dynamics of the ice
body, grounding line dynamics, calving, and basal mass balance. Basal friction does not act on floating ice as it does on grounded ice due to the fairly
low dynamic friction at the ice-water boundary. The flow of grounded ice
sheets is dominated by vertical shear stresses, while ice shelf flow is dominated by longitudinal stretching and lateral shearing. An ice shelf behaves
like a large membrane which spreads by horizontal stretching, while vertical
shearing is unimportant [Schoof, 2007b]. Calving and the basal mass balance
are important processes for the overall mass balance of the ice sheet because
they control the mass loss.
The grounding line divides the dynamics of the ice shelves from those of
the grounded ice. Grounding line dynamics affect the evolution of both. The
horizontal extent of the transition zone is the distance over which longitudinal
stresses from the ice shelf propagate into the ice sheet. A slight retreat in
the grounding line position will lead to an increase in ice thickness and hence
ice discharge at the grounding line. This represents a positive feedback as anincrease in ice discharge leads to further shrinkage of the ice sheet and hence
grounding line retreat. This then leads to a further increase in ice discharge.
The back force (the buttressing effect) is the total force transmitted upstream by forces, other than the pressure of sea water at the ice front, that
oppose the spreading of the ice shelf. The drag exerted by the side walls andice rise are examples of such forces.
The collapse of the Larsen ice shelves has been an impressive phenomenonof recent years. On the Antarctic Peninsula, ice shelves have lost ∼ 300 km2
in a year. Larsen A ice shelf (1, 600 km2 ) broke away in 1995 and Larsen B
3

ice shelf (3, 500 km2 ) in 2002. Such collapses may occur in other areas under
future global warming conditions.
By using geological proxy data, Scherer et al. [1998] conjectured that the
upstream part of the Ross Ice Shelf collapsed at least once in the Pleistocene.
The collapse could have occurred during past interglacial periods. Thus, ice
shelf dynamics has probably played an important role in ice sheet evolution
during glacial-interglacial cycles.

1.3

Global temperature and sea level change

It is a major topic of discussion that the Earth’s surface climate is changing.
One cause of change is considered to be human activities, such as burning
fossil fuels and land management practices leading to as deforestation. The
effects have become evident in recent decades. In consequence, people havestarted to recognize the importance of scientific and social studies related tothis environmental problem.
The Intergovernmental Panel on Climate Change (IPCC) is the leadinginternational community for the assessment of climate change. It was established by the United Nations Environment Programme (UNEP) with theWorld Meteorological Organization (WMO) to provide the world with a clear
scientific view of the current state of knowledge on climate change and its
potential environmental and socio-economic impacts.
Researchers have made progress in understanding the climate changemechanism in space and in time through improvements and extensions of
numerous datasets and data analyses (Fig. 1.4). Increasingly comprehensive
observations are available for glaciers and snow cover going back to the 1960s,
and for sea level and ice sheets over about the last 10 years, although data
coverage is still limited. The IPCC AR4 (fourth assessment report) shows
unequivocally that the climate has warmed over the last 100 years and is
still warming. Observations show an increase in global average air and ocean
temperatures, widespread melting of snow and ice, and rising mean sea level[Solomon et al., 2007].
Future climate change can also be estimated by climate model results.
Many scenarios are considered when estimating the likely range of future
warming. The important consequence is that global warming is inevitableeven if the concentrations of all greenhouse gases and aerosols are kept constant at the levels of 2000 (Fig. 1.5). It is predicted that a further warmingof about 0.1◦ C per decade is to be be expected even this stabilization isachieved. The following scenarios have been considered:• A1: characterized by rapid economic growth. Global population reaches
4

9 billion in 2050 and then gradually declines. International interactions
are extensive. There are several types of experiments; A1F1 emphasizes fossil fuels, A1B has a balanced emphasis on all energy sources,
and in A1T non-fossil energy sources dominate.
• A2: inter-regional interactions are limited and regions are self-reliant.
Population increases continuously.
• B1: the world is more integrated and more ecologically friendly. Economic growth and population trends are the same as in A1.
• B2: the world is more divided than in A2, but more ecologically friendly
than in B1. The population increases continuously but more slowly
than in A2.
There are variations between experiments but all scenarios show an upward
trend. Six scenarios have been used to estimate the global average sea level
rise by the end of the 21st century (Fig. 1.6). Thermal expansion contributes
70-75% to the best estimate for each scenario. Glacier and ice sheet melting
account for the rest.
In some parts of West Antarctica, large accelerations of ice flow haverecently occurred which may have been caused by thinning of ice shelves
due to ocean warming. The limited understanding of such processes was
explicitly spelled out in IPCC AR4: “Dynamical processes related to ice flow
not included in current models but suggested by recent observations could
increase the vulnerability of the ice sheets to warming, increasing future
sea level rise. Understanding of these processes is limited and there is noconsensus on their magnitude.”
Although this has not been formally attributed to anthropogenic climatechange due to greenhouse gases, it suggests that future warming could cause
faster mass loss and greater sea level rise. Quantitative projections of thiseffect cannot be made with confidence. If recently observed increases in ice
discharge rates from the Greenland and Antarctic Ice Sheets were to increaselinearly with global average temperature change, they would add 0.1 to 0.2 m
to the upper bound of sea level rise. Understanding of these effects is, as yet,
too limited to assess their likelihood or to give a good estimate.

5

1.4

Changes in the Antarctic ice sheet andice shelves

The typical scales of the ice sheet are
typical horizontal extent [L] = 1000 km;
typical horizontal velocity [U ] = 10 − 100 m/a;
typical time scale [t] = [L]/[U ] = 104 − 105 a.

(1.1)

Since the time scale of an ice sheet is quite long, it has long hysteresis andchanges slowly. Ice cores drilled at Dome C and Dome Fuji in the EAIS show
eight glacial cycles [Augustin et al., 2004, Goto-Azuma, 2008]. The temperatures of past environments can be determined from measurements on ice
cores taken from boreholes. These measurements suggest that the Antarctic
ice sheets have not changed significantly in 100-1000 years. However, thisview is now changing.
The mass balance of Antarctica and its changes have only recently become
known. The observational record of Antarctic glaciers is short and limited
because humans have seldom visited Antarctic glaciers. But the advent ofsatellite missions has recently revolutionized our knowledge of ice dynamics
in Antarctica (Fig. 1.7). The observations suggest the current mass balanceof the ice sheet is most likely negative, with an accelerating trend [Rignot
et al., 2011b] (Fig. 1.8). Owing to the very low surface temperatures, the ice
sheet is not very susceptible to surface melting even under moderate global
warming scenarios. However, recent observations have led to strong concernsthat dynamic ice processes (loss of buttressing from ice shelves, speed-up
of ice streams and outlet glaciers) may boost mass loss and thus make an
additional contribution to sea level rise.
Recent observations [e.g., Joughin et al., 2003] have led to strong concern about ice-dynamic processes. One of the most impressive phenomena
affecting ice shelves in recent years is the collapse of the Larsen ice shelves.
Very small changes in the ice shelves can have large effects on regions of ice
upstream.
Ice shelves are among the most active parts of the ice sheet system. Theyinteract with the inland ice sheet, glaciers and ice streams which flow into
them, and with the sea into which they eventually melt, either directly or as
icebergs. The typical timescale of an ice shelf is shorter than that of an ice

6sheet:typical horizontal extent [L] = 100 km;
typical horizontal velocity [U ] = 100 − 1000 m/a;
typical time scale [t] = [L]/[U ] = 100 − 1000 a.

(1.2)

Ice shelves therefore change more rapidly than ice sheets and they also affectthe upstream inland ice.
Rignot [2006] summarized the effect of dynamic changes in ice sheet margins. The surface velocity of glaciers around the Larsen ice shelf has changed
considerably (Fig. 1.9). The glaciers have accelerated markedly and left bitsof glacier ice hanging on cliff walls 40 m above the new glacier surface after
the collapse of Larsen A in 1995 [De Angelis and Skvarca, 2003]. Drygalski Glacier flowed three times faster in 2000 than before the collapse [Rottet al., 2002]. Glacier acceleration continued for more than 10 years after theremoval of the ice shelf [Rignot, 2006]. Following the collapse of the Larsen
B ice shelf, Hektoria/Green/Evans flowed eight times as fast, and Crane andJorum twice as fast as they had done before [Scambos et al., 2004].
Pine Island Glacier (Fig. 1.10) is an important drainage of the West
Antarctic ice sheet. Its grounding line is retreating rapidly, at 1000 m/a.
The ice shelf below the Pine Island glacier has the highest rates of bottom
melting in Antarctica. These are are 58 ± 8 m/a near the grounding line
and 24 ± 4 m/a on average for the whole ice shelf [Rignot, 2006]. This is
two orders of magnitude greater than the melting underneath the Ross andRonne ice shelves. This suggests that thermal forcing from the ocean has
a great influence on the ice shelf mass budget [Jenkins et al., 1997]. Interaction between oceans and ice shelf is therefore considered to be the causeof grounding line retreat. The glacier is accelerating and this acceleration
is sufficient to explain the observed thinning [Rignot et al., 2002]. The estimated velocity field and changes from the 1970s to the present are shownin Fig. 1.11. Glacier velocity increased 38% from 1974-2005. The acceleration is detected not only at the grounding line but also more than 100 km
upstream. Changes of flow in ice shelves and ice streams were also foundin East Antarctica. All ice sheet margins have been changing their flow inrecent decades.
The transition zone between grounded and floating ice divides the ice
shelf and the ice sheet. In most cases, ice streams are formed beneath ice
shelves and flow into them. The discontinuity between the dynamics of theice sheet and those of the ice shelf in ice sheet models causes a problem in
representing its dynamics [Pattyn et al., 2006]. A suitable treatment of thetransition zone is important for the dynamics of the coupled system, and for7marine ice sheet dynamics in particular [Weertman, 1994, Schoof, 2007a].
The boundary layer theory [Schoof, 2007a] is applied to the groundingzone for discussing the role of grounding line dynamics in marine ice sheet
evolution. The boundary layer theory states that ice flux through the grounding line in a two-dimensional sheet-shelf system increases sharply with ice
thickness at the grounding line. This result is then applied to the large-scale
dynamics of a marine ice sheet.
The results are that steady grounding lines cannot be stable on reverse
bed slopes [Vieli and Payne, 2005]. They can also undergo hysteresis whenmaterial parameters (such as ice viscosity or basal sliding) or external forcing(accumulation or sea level) vary. The hysteresis mechanism, which is a uniquecharacteristic of the marine ice sheet, is driven by outflow of ice through thegrounding line and requires overdeepening. It is also suggested that thisbehaviour can, in principle explain the retreat of the WAIS following the lastglacial maximum, and may play a role in the dynamics of Heinrich events[Schoof, 2007a], which are the quasi-periodic episodes of massive ice discharge
from the Laurentide ice sheet [Heinrich, 1988].

1.5

Ice sheet model SICOPOLIS

Ice sheet models use quantitative methods to simulate the evolution, dynamics and thermodynamics of ice sheets. The gravity-driven activity of ice
sheets is controlled by two main variables: the strength of the ice sheet baseand the temperature. Glacier flow is controlled by internal deformation of
ice and by basal sliding. In an ice sheet model, the momentum balance iscalculated to estimate velocities of ice with continuum mechanics. Ice sheet
evolution is simulated using calculated velocities and the mass balance atboundaries. Jenssen [1977] constructed a three-dimensional ice sheet modeland applied it to the Greenland ice sheet (Fig. 1.12). The glacial-interglacial
evolution of the Antarctic ice sheet with coupled ice shelves was simulatedby [Huybrechts, 1990]. Some ice sheet models have been constructed more
recently [e.g., Payne, 1999, Ritz et al., 2001, Saito et al., 2003, Bueler and
Brown, 2009]. Pollard and DeConto [2009] simulated five million years of
evolution of the ice sheet. Although most ice sheet models use the finite
difference method due to its computational efficiency, some ice sheet models
use the finite element method [Seddik et al., 2012, Larour et al., 2012]. Such
finite element models use higher order (less approximation) or full stokes field
equations. However, it takes a lot of computational effort to compute theentire ice sheet in higher order dynamics.
SICOPOLIS (http://sicopolis.greveweb.net/) is a three-dimensional dy8

namic/thermodynamic ice sheet model. It can be used to simulate the evolution of large ice sheets. It was originally created as part of the doctoral thesis
of Greve [1995] in a version for the Greenland Ice Sheet. The model solves
the thermomechanically coupled field equations: the equations of mass, momentum and energy, along with the constitutive equation, Glen’s flow law.
Boundary conditions for the free surface and the bed are also included. Isostatic displacement of the bed due to ice loss is also calculated for an ice
sheet. It includes a detailed treatment of basal temperate layers. Within
the temperate layers, water content is calculated, and its influence on the ice
viscosity is taken into account [Greve, 1997b]. The model has been applied tolarge ice sheets such as the Greenland, Antarctic, Laurentide and Fennoscandian ice sheets [Greve, 1997a, Forsström and Greve, 2004] (Fig. 1.13).
Greve [1997a] applied the model to the Greenland ice sheet. This model
was applied here to obtain the steady state under present climate conditions
with different parameters. Transient simulations were also conducted. Past
glacial experiment and future climate runs were carried out.
Calov et al. [1998] applied the model to the Antarctic ice sheet. The
results are presented of a transient paleoclimatic simulation which covers
two entire glacial and interglacial cycles, with special emphasis on Dronning
Maud Land, to study ice coring regions and proposed possible drilling sites byusing depth-age relationships derived from the model. Although it has been
applied to the Antarctic ice sheet, special treatment of ice shelf dynamics
has not yet been included.

1.6

Purpose and tasks of the study

To study Antarctic ice sheet dynamics and their changes, a model of ice shelf
dynamics was constructed and coupled to the ice sheet model SICOPOLIS.
Ice shelf dynamics are considered and coupled with the ice sheet model.
First, an ice shelf model was made without the coupled ice sheet model.
The model’s performance was then checked. It was then coupled to the ice
sheet model. Not only differences in the dynamics, but also differences of flow
regimes and stresses working in the ice shelves cause differences in heat flow
and ice thickness evolution. Future evolution of the ice sheet was estimatedusing contributions from this project and newly-derived ice sheet data. Sealevel Response to Ice Sheet Evolution (SeaRISE) is a community-organizedeffort to estimate the upper bound of ice sheet contributions to sea level overthe next few centuries. The coupled model is applied to the Antarctic ice
sheet using the forcings given by the project.
As stated before, the ice sheet has changed significantly during glacial9

interglacial climate cycles. Modeling the paleo-ice sheet allows comparisonof past and present ice sheet conditions. Paleoclimatic records such as icecore records help in understanding past ice sheet and climatic conditions.
Past ice sheet evolution was simulated using information gained from such
paleoclimatic records. How the Antarctic ice sheet will change in the future is
an important question relating to future global warming. The mass balanceof Antarctica and its changes are not very clear even today. One reason forthis is the uncertainty arising from ice shelf dynamics. The effects of the
Antarctic ice sheet was studied in cooperation with SeaRISE.

10

Figure 1.1: A map of the Antarctic ice sheet with ice shelves [Williams et al.,2005].

11

Figure 1.2: Typical conditions of the marine ice sheet with ice shelves (modified from Schoof [2007b]). The horizontal arrow is the direction of ice flow.
xg is the grounding line and xc is the calving front. Although ice thickness
h(x) is greater than floating shelves, ice base elevation b(x) is still lower thansea level.

12Figure 1.3: Ice shelf geometry in the horizontal map plane [Greve and Blatter,2009]. Incoming flow comes from the upstream ice sheet. The boundary of
the ice shelf and ice sheet is the grounding line. Iceberg calving occurs infront of the ice shelf then ice disperses in the sea. The boundary of the ice
shelf and the sea is the calving front.

13

Figure 1.4: Observed changes for March-April in (a) global average surface
temperature, (b) global average sea level from tide gauge (blue) and satellite(red) data and (c) Northern Hemisphere snow cover for [Solomon et al., 2007].
All these changes are relative to the corresponding averages for the period
1961-1990. Smoothed curves represent decadal average values. White circles
show yearly values. The shaded areas are the uncertainty intervals estimatedfrom a comprehensive analysis of known uncertainties (a and b) and fromthe time series (c).

14

Figure 1.5: Multi-model means of surface warming (compared to the 19801999 base period) for each of the scenarios A2 (red), A1B (green) and B1(blue), shown as continuations of the 20th-century simulation [Solomon et al.,2007]. The A1B and B1 scenarios are continued beyond the year 2100 with
forcing kept constant. An additional experiment, in which the forcing is
kept at the 2000 level is also shown (orange). Shading means ±1σ of themulti-model means.

15Figure 1.6: Projected global average sea level rise and its components in 2090
to 2099 (relative to 1980-1999) for the six scenarios [Solomon et al., 2007].

16

Figure 1.7: The velocity field of the Antarctic ice sheet [Rignot et al., 2011a].
It is derived from ALOS PALSAR, Envisat ASAR, RADARSAT-2 and ERS1/2 satellite radar interferometry. Thick black lines delineate major ice divides. Thin black lines outline subglacial lakes discussed in the text.

17

Figure 1.8: Total of the Antarctic ice sheet mass balance, dM/dt, between
1992 and 2009 in Gt/yr. The black circle is derived from the mass budget method (MBM). The red triangle is derived from GRACE time-variablegravity. The error bar and the trend (black (MBM) and red lines (GRACE))
are also shown [Rignot et al., 2011b].

18

Figure 1.9: Glacier acceleration after the collapse of Larsen A and B for (a)Drygalski, (b) Hektoria/Green/Evanst [Rignot, 2006]. Black square is 1996,
red square is 2000, blue triangle is 2003, yellow triangle is 2003, green triangle
is 2003 and purple diamond is 2005.

19

Figure 1.10: The velocity field of Pine Island bay on the WAIS [Rignot, 2006].
It is derived from ERS-1/2 1996 interferometry, overlaid on a RADARSAT-1
map of radar brightness. The grounding line is thin black.

20

Figure 1.11: Acceleration of the Pine Island glacier from 1974, 1987, 1992,
1996 (ERS-1/2) and 2001-2005 (RADARSAT-1) versus distance along theglacier flow line [Rignot, 2006].

21

Figure 1.12: The domains of the ice sheet model (black line) and a maximum
time step (year) for the simulation. The horizontal grid spacing is 100 km
for west-east and 200 km for north-south. [Jenssen, 1977].

22

Figure 1.13: Examples of the application of the ice sheet model SICOPOLIS
to (a) the steady state velocity of Greenland ice sheet for the present climate[Greve, 1997a].

23

Chapter 2Theory2.1

Dynamics of large ice bodies

The typical geometry (cross-section) of an ice sheet consists of a grounded
ice sheet with a floating ice shelf attached (Fig. 2.1) (The ice shelf will betreated, as well as its interactions with the atmosphere (snowfall, melting),
the lithosphere (geothermal heat flux, isostasy) and the ocean (melting andcalving)).
Ice flow is slow enough for the acceleration term in Newton’s Second Law
to be neglected. This means that a glacier is in static equilibrium even when
it is moving. The forces working in a glacier and an ice sheet are in balance.
Gravity pulls a glacier vertically downward. The resultant driving forces
mainly depends on the slope of the ice surface. Consider an ice body withthickness H and surface slope α. The driving force, τd , isτd = ρgH sin α,(2.1)

where g is the gravitational constant, ρ is the density of ice; which is considered as the constant in this study. The largest forces resisting this stress
are the basal shear stresses, τb . Other components are the side wall effects(τw ,the side drag) and stretching and compression ice (τL , the longitudinaldrag). Therefore, the total momentum balance isτd = τb + τw + τL .

(2.2)

This is the force balance acting in an ice sheet.

2.2

Constitutive equations

Ice crystals have lattice defects. This is because the system has higher entropy with defects than without. Dislocations, a type of defect movement,
24causes the deformation of ice crystals (Fig. 2.2). There is a relationship between dislocation theory and the constitutive equation [e.g., Gudmundsson,2007]. Consider a crystal with height h, width w, and thickness l. Burger’s
vector of a dislocation is b. When n dislocations with length l, moved a
distance x in the crystal, the shearing deformation e of the crystal is
nb x
.
h w

e=

(2.3)

Dislocation density ρ is defined as the sum of dislocations in a unit cell. Sincea dislocation is a line defect, its density is
nl
n
=
.
hwl
hw

(2.4)

with the result that the time differentiation of the shear deformation
equation isė = ρbẋ = ρbv.
(2.5)
Dislocation density ρ is inversely proportional to the square of dislocation
distance d,ρ ∝ d−2 .
(2.6)
d is related to internal stress σ by Hooke’s law,
σ∝

µb
,
d

(2.7)

where µ is the rigidity modulus.
Those two equations imply the Bailey-Hirsch relationshipσ∝

σ2.
µ 2 b2(2.8)

because dislocation velocity is proportional to the stress,
v = Cσ.

(2.9)

With (7), (8) and (4), the strain rates are proportional to the third power ofthe internal stress,ė = Cσ 3 .
(2.10)
The constitutive equation for a viscous material is
ė =

1
.
2η(T, p, σe )25

(2.11)

The inverse viscosity is called fluidity, and its dependence on the temperature
T, the pressure p and the absolute value of σe can be factorized as
1= 2A(T, p)f (σe ) = 2A(T, p)σen−1
η(T, p, |τ |)(2.12)

where A(T,p) is called the rate factor and f (|τ |) is the creep function. Thisequation and equation (9) give Glen’s flow law,ė = A(T, p)σe2 t.

(2.13)

Glen’s flow law yields an infinite viscosity when the effective stress approaches zero. This is not a problem physically because the strain rate issmall when the effective stress is small. It therefore contributes little to large
scale flow. However, depending on the mathematical solution procedure, theinfinite viscosity limit may introduce a singularity in the equations for thevelocity field if the effective stress is very small, which occurs at ice divides
and ice margins. In order to avoid this problem, a regularization has been
proposed in which the power law is replaced by a polynomial relationship,
f (σe ) = σen−1 + σ0n−1 ,(2.14)

where the residual stress σ0 is a small positive constant. This is the “regularized Glen’s flow law” and viscosity is thenη(T 0 , σe ) =1,
2A(T, p)(σen−1 + σ0n−1 )(2.15)

which has finite viscosity at the limit σe → 0,
ησe →0 =2.3

1.
2A(T, p)σ0n−1(2.16)

Stress balance equations in general coordinates

The general coordinate invariance is one of the important principles of generalrelativity. The classical field theory gives
Gij + Λg ij =8πG ijT ,c4(2.17)

where Gij is the Einstein tensor, Λ is the cosmological constant, g ij is themetric tensor, G is Newton’s gravitational constant, c is speed of light in a26

vacuum, and T ij is the stress energy momentum tensor. It uses Einstein’s
convention (sum over the coordinates, i = 1 to 4). It is also possible tounderstand Navier Stokes and Stokes equations as a form of the classicalfield theory, the field equation. The energy and momentum balance equationisT ij ;j = 0.
(2.18)(This equation is valid in a local system. It is also necessary to consider
a term for the energy momentum caused by gravity [Landau and Lifshitz,
1980].) The space components (i = 1 to 3) of the energy momentum tensor
represent stresses. Note that ; means the covariant derivative, , is the ordinary derivative and Γ is the Christoffel symbol. How the energy momentum
tensor should be formulated depends on the problem. In a viscous fluid in
Cartesian coordinates, T ij isT ij = −pg ij + η(ρ + p)ui uj + τ ij .

(2.19)

The form of τ ij in the dissipative system is given by the second law of thermodynamics [cf. Landau and Lifshitz, 1987]. In Cartesian coordinates,
 i
 l
 
i
∂u
∂uj
∂u ij
2
ij
l ∂u
t = −η
+ i + uj u
(g + ui uk ).
(2.20)− ζ− η
j
l
l
∂x
∂x
∂x
3
∂x
In a perfect fluid,
T ij = (ρ + p)ui uj + g ij p,(2.21)

where ρ means the mass energy density and p the pressure. The zero component is the mass balance,


1 ∂ρ
0j
j
+ ρu j .
T,j =(2.22)c ∂tThis provides the mass and energy balance.
Since an ice sheet is relatively slow, it is possible to neglect the acceleration term of momentum balance. This produces Stokes equation which, in
general coordinates, can be formulated as
ij
j ik
i kj
j i
tij
;j = t,j + Γkj t + Γkj t = f δj ,(2.23)

where f j δji is the external force. In this study, the external force is gravity.
The stress balance in large-scale ice sheet dynamics is discussed in this section. The constitutive equation for the general ice sheet model is shown in
the figure. The figure is shown in this section. The external gravitational
force vector is
f i = g ij fj = (0, 0, ρG).
(2.24)
27

2.4

Large-scale ice sheet dynamics

In most parts of a large ice sheet, the flow regime is essentially simple, bedparallel shear. The geometry of real ice sheets is such that they are shallow,
i.e. the typical horizontal length scale exceeds the typical ice-thickness scaleby several orders of magnitude. The slopes of the free surface and the ice baseare small. The typical thickness of the ice sheet is [H] = 1000 m, Therefore,
the typical scale of the surface slope is[H]103= 6 = 10−3 .
[L]
10(2.25)

D
D
Under this condition, the normal stress deviations tD
xx , tyy and tzz are negligible compared to vertical shear stresses,[v][η] [L][H]tD
xx== 10−3 .
=
[v]
D
txz
[L]
[η]

(2.26)[H]

Since the shape of the Earth or a planet in general is close to a sphere,
the atmospheric and ocean general circulation models consider the effect of
the curvature. One way to treat the effect is to apply orthogonal curvilinear coordinates. Consider an ice sheet on a planet surface. In most cases,
ice thickness is very small compared to the distance from the surface to thecore. The difference of the curvature in “the vertical direction”, which isthe direction of the gravitational force, is negligible. This allows the use oforthogonal curvilinear coordinates which have a gravitational force in thevertical direction of γ, and two horizontal coordinates which are perpendicular to the vertical direction, α and β. This represents a shallow ice body ona plane which has curvilinear horizontal coordinates.
In orthogonal curvilinear coordinates, the metric tensor only has diagonal
elements. As a consequence of the above, the vertical component of thecoordinate can be set to g3 = 1. Since the curvature of time is not consideredin this situation, it is unnecessary to consider the zero dimension of the metrictensor. Therefore, the metric tensor is
g11 0 0gij =  0 g22 0 
0
0 1

2.5

Hydrostatic approximation

The scale of the shear stresses (less than 100 kPa) is smaller than the vertical normal stress (10 MPa). Therefore, the stress balance for the vertical
28

coordinate can be approximated by
∂i T

γi

1 ∂tγ γ
= ρG.
=√
g ∂γ

(2.27)

Since atmospheric pressure is negligible compared to the pressure in the ice,the vertical normal stress on the surface, h, istγ γ |γ=h = 0.

(2.28)The balance equation in the vertical direction can be integrated from thedepth, γ, to the surface, h,tγ γ = −ρG(h − γ).

(2.29)

Therefore, the vertical normal stress can be approximated to the hydrostaticpressure in the orthogonal curvilinear coordinate system as in Cartesian coordinates. Pressure is formulated as
p = ρG(h − γ) − (tD )α α − (tD )β β .

(2.30)

This results in the horizontal normal stress changing totα α = 2(tD )α α + (tD )β β − ρG(h − γ),tβ β = 2(tD )β β + (tD )α α − ρG(h − γ).

(2.31)

Then the horizontal stress balance equations are expressed as∂(tD )α α ∂(tD )β β ∂tβ α ∂tγ α
1 ∂g2 D α2++++((t ) α − (tD )β β )
∂α
∂α
∂β
∂γ
2g2 ∂α1 ∂g1 α1 ∂g2 β∂(h − γ)+t β+t α = ρG,
2g2 ∂α2g2 ∂β∂α∂(tD )β β ∂(tD )α α ∂tαβ
∂tγβ
1 ∂g1 D β2++++((t ) β − (tD )α α )
∂β
∂β
∂α
∂γ
2g1 ∂β1 ∂g2 β1 ∂g1 α∂(h − γ)+t α+t β = ρG.
2g1 ∂β2g1 ∂α∂β(2.32)(2.33)(2.34)(2.35)

This is the first order hydrostatic approximation of the stress balance in thiscoordinate system. The vertical velocity is estimated with the incompressible
mass balance equation,
T 0j ;j = 0.
29

(2.36)

On the viscous ice(tD )ij = 2ηDij .

(2.37)The relationship of strain rates to velocities is
1
Di j = (ui ;j + gkj g il uk ;l ).
2(2.38)

It is also possible to apply further approximations. The horizontal difference
between the vertical velocities is negligible [Greve and Blatter, 2009],
uγ ;α = 0,uγ ;β = 0.

(2.39)

Finally, the stress balance equation is
2∂(tD )α α ∂(tD )β β ∂tβ α ∂tγ α
1 ∂g2 D α++++((t ) α − (tD )β β )
∂α
∂α
∂β
∂γ
2g2 ∂α1 ∂g2 β∂(h − γ)1 ∂g1 αt β+t α = ρG,+
2g2 ∂α2g2 ∂β∂α∂(tD )β β ∂(tD )α α ∂tαβ
∂tγβ
1 ∂g1 D β2++++((t ) β − (tD )α α )
∂β
∂β
∂α
∂γ
2g1 ∂β1 ∂g2 β1 ∂g1 α∂(h − γ)+t α+t β = ρG.
2g1 ∂β2g1 ∂α∂β(2.40)(2.41)(2.42)(2.43)

which is the first order approximation. It is important to set the verticalaxis to set the direction of gravitational force. Otherwise, the hydrostatic
approximation (first order approximation) does not work as in the Cartesian
coordinates. The reason is that in that case the gravitational force vectors
are not zero in the horizontal coordinates. If the distance is not very great,tγγ is larger than the other part of the stress balance. Although the gravitational term in the horizontal coordinate needs to be considered, the conceptof the hydrostatic approximation may not be compromised if the verticalnormal stress is sufficiently large. However, if the horizontal distance is not
negligible, the validity of the theory itself is damaged. In other words, thefirst order approximation loses its validity when the space is sufficiently large
and Cartesian coordinates alone are used. It is therefore undesirable to use
Cartesian coordinates in this type of approximation for continental or global
scale conditions. An alternative strategy is to consider the directions of thegravitational force in the force balance. In this case, it is preferable to estimate the ratio of shear stress to normal stress and the gravitational term inhorizontal directions to check validity.

30

2.6

Shallow shelf approximation

2.6.1

Formulation in orthogonal curvilinear coordinates

Since there is no basal drag beneath ice shelves, there is no basal shear
stress. The effects of the side wall force and longitudinal resistant force arecomponents of the total local force balance,
τd = τw + τL .

(2.44)The difference from Eq. (2.2) is that it does not have τb . Forces acting at
the shelf edges and sides influence the flow even at locations far from theseboundaries since the local basal drag is negligible. To construct meaningful analyses of shelf dynamics, the force balance relationship needs to be
integrated over the distance from a point in the shelf to the calving front.
The typical situation is that of plug flow, whereby the horizontal velocities
are essentially constant over depth (Fig. 2.4). Naturally, in the immediate
vicinity of the grounding line, a transitional state between these two limits
will develop. Mathematically, the regime of plug flow can be defined by
∂uα
≈ 0,
∂γ
∂uβ
≈ 0.
∂γ

(2.45)

In other words, horizontal velocities depend only on the horizontal coordinates and on time,uα = uα (α, β, t), uβ = uβ (α, β, t).

(2.46)In orthogonal curvilinear coordinates, basal friction is negligible as in
Cartesian coordinates. Moreover, the direction of the gravitational force
is parallel to the vertical direction, the horizontal velocities do not dependon the vertical direction. Mathematically speaking, the Jacobian is 1 sinceg3 = 1. It is therefore possible to apply the shallow shelf approximation, thevertically integrated stress balance, in orthogonal curvilinear coordinates.

2.6.2

Boundary conditions at the ice surface

The surface function, Fs , is defined as
Fs (α, β, t) = γ − h(α, β, t) = 0.
31

(2.47)

The normalized vector gradient is
Fs;i
1
ni =
=
|Fs;i |
Ns




∂h ∂h
− ,− ,1 ,∂α ∂β(2.48)

where Ns is
Ns = |Fs;i |.

(2.49)

tj i gj = 0.

(2.50)The stress free condition is

2.6.3

Boundary conditions at the ice base

The momentum jump condition is(tjsea )i gj − tj i gj + gij ρab [uj ] = 0,(2.51)(tjsea )i gj = −psea gi + (τsea )j i gj
= −ρG(γsl − b)gi + Cwi ρsw |usea |2 (gτ )i ,(2.52)

where

where (gτ )i is the tangent vector. Since the shear stress from ocean currents
is negligible, this equation changes to
tj i gj = ρG(γsl − b)gi .

2.6.4(2.53)

Vertical integration

This characteristic of ice shelf dynamics can be described well by the shallow
shelf approximation. Since an ice shelf has no basal friction, the vertical
shear stress is negligible. Therefore, horizontal velocities do not dependon the vertical coordinate. To make use of this property, the momentum
balance equations are integrated over depth. In the orthogonal curvilinear
coordinates with shallow conditions, it is possible to treat vertical integration32

as the Cartesian coordinate. The vertically integrated momentum balance is
1 ∂g2 α∂N α α ∂N β β ∂N β α ∂N γ α++++(N α − N β β )
2
∂α
∂α
∂β
∂γ
2g2 ∂α


1
∂g1 α
∂g2 β
∂h
∂b
+
N β+
N α − 2(tD )α α |γ=h
+ 2(tD )α α |γ=b2g2 ∂α
∂β
∂α
∂α
∂b
∂h
∂b
∂h
+ tβ α |γ=b
− (tD )β β |γ=h
− (tD )β β |γ=b−tβ α |γ=h
∂β
∂β
∂α
∂α
∂(h − γ)+tγ α |γ=h − tγ α |γ=b = ρG,∂α∂N β β ∂N α α ∂N α β ∂N α β
1 ∂g1 β2++++
(N β − N α α )
∂β
∂β
∂α
∂α
2g1 ∂β


∂g2 β
1
∂g1 α
∂h
∂b
+
N α+
N β − 2(tD )β β |γ=h
+ 2(tD )β β |γ=b2g1 ∂β
∂α
∂β
∂β
∂h
∂b
∂h
∂b
−tα β |γ=h
+ tα β |γ=b
− (tD )α α |γ=h
− (tD )α α |γ=b
∂α
∂α
∂β
∂β
∂(h
−
γ)
+tγ β |γ=h − tγ β |γ=b = ρG∂βThe membrane stress N , is defined as
Z h
ij
N =
(tD )ij dγ.

(2.54)(2.55)b

The boundary condition on the surface and the base can be considered as in
the Cartesian coordinate. The surface condition is∂h∂h∂h− (tD )β β |γ=h
− tβ α |γ=h
+ tγ α |γ=h = 0,−2(tD )α α |γ=h
∂α
∂α
∂β
∂h
∂h
∂h
−2(tD )β β |γ=h
− (tD )β β |γ=h
− tβ α |γ=h
+ tγ β |γ=h = 0,(2.56)∂β∂β∂αtγ γ |γ=h = 0.
The boundary condition on the ice base is∂b∂b∂b−2(tD )α α |γ=b
− (tD )β β |γ=b
− tβ α |γ=b
+ tγ α |γ=b
∂α
∂α
∂β
ρsw G
∂b
=
(γsl − b) ,
g1
∂α
∂b
∂b
∂b
−2(tD )β β |γ=b
− (tD )β β |γ=b
− tβ α |γ=b
+ tγ β |γ=b
∂β
∂β
∂α
ρsw G
∂b
=
(γsl − b) ,g2∂βγt γ |γ=b = ρsw G(γ − b).
33(2.57)

It is possible to replace ρsw G(γsl − b) to ρGH by using the hydrostatic condition,
ρH = ρsw (γsl − b).
(2.58)
Thus, the boundary condition at the ice base is∂b∂b∂b− (tD )β β |γ=b
− tβ α |γ=b
+ tγ α |γ=b = 0,
∂α
∂α
∂β
∂b
∂b
∂b
−2(tD )β β |γ=b
− (tD )β β |γ=b
− tβ α |γ=b
+ tγ β |γ=b = 0,∂β∂β∂αtγ γ |γ=b = ρH.

−2(tD )α α |γ=b(2.59)

With the condition of 2.56 and 2.59, 2.54 becomes∂N α α ∂N β β ∂N β α ∂N γ α
1 ∂g2 α2++++(N α − N β β )
∂α
∂α
∂β
∂γ
2g2 ∂α1 ∂g2 β∂(h − γ)1 ∂g1 αt β+t α = ρG,+
2g2 ∂α2g2 ∂β∂α∂N β β ∂N α α ∂N α β ∂N α β
1 ∂g1 β2++++
(N β − N α α )
∂β
∂β
∂α
∂α
2g1 ∂β1 ∂g1 α∂(h − γ)1 ∂g2 βt α+t β = ρG.
+
2g1 ∂β2g1 ∂α∂β(2.60)(2.61)(2.62)(2.63)

The constitutive equation for viscous ice when vertically integrated is
N i j = 2η̄Di j ,where η is the vertically integrated viscosity,ZZ h
1 −(1−1/n) hη̄ =ηdγ = deB(T 0 )dγ,2bb(2.64)(2.65)and D is the strain rate tensor,
1
Di j = (ui ;j + gkj g il uk ;l ).
2(2.66)The momentum balance equation in Cartesian coordinates is




 

∂
∂vx
∂
∂vy
∂
∂vx ∂vy
∂h
4
η̄
+2
η̄
+
η̄
+
= ρgH ,
∂x
∂x
∂x
∂y
∂y
∂y
∂x
∂x




 

∂
∂vy
∂
∂vx
∂
∂vx ∂vy
∂h
4
η̄
+2
η̄
+
η̄
+
= ρgH . (2.67)
∂y
∂y
∂y
∂x
∂x
∂y
∂x
∂y
34

This is called shallow shelf approximation [Greve and Blatter, 2009].
The vertical velocity isZ z∂vx ∂vy+.
vz = vz |z=zsl −
∂x
∂y
zsl


∂vx ∂vy= vz |z=zsl − (z − zsl )+.
(2.68)∂x∂y
The difference between this and the shallow ice approximation is the different
reference horizons. The reference for ice shelves is sea level but for ice sheets
it is the base of the ice sheet. Moreover, precise integration is unnecessary
because horizontal velocities do not depend on the vertical direction.

2.7

Shallow ice approximation

The driving force balances basal shear stresses in most parts of a large ice
sheetτd = τb .
(2.69)The local stress balance does not depend on non-local effects, wall drag andlongitudinal drag. This is a characteristic of large ice sheets.
The total hydrostatic pressure approximates top = ρgH = ρg(h − z).

(2.70)and the horizontal components of the momentum balance simplify to
∂p
∂h
∂txz
=
= ρg ,
∂x
∂x
∂x
∂tyz
∂p
∂h
=
= ρg .
∂y∂x∂y(2.71)

Giving the shallow ice approximation [Greve and Blatter, 2009].
With Glen’s flow law and the above concept, the horizontal velocities are
Z z
3
2 ∂h
vx = vbx − 2(ρg) |∇h|A(T, p)(h − z)3 dz,
∂x b
Z z
3
2 ∂h
vy = vby − 2(ρg) |∇h|A(T, p)(h − z)3 dz.
(2.72)∂y b
The vertical velocity can be computed by integrating the continuity equation from z = b to z,Z z∂vx ∂vyvz = vz |z=b −+.
(2.73)
∂x
∂y
b
35

2.8

Temperature evolution of ice

Because viscosity depends on temperature, a thermo-mechanically coupledproblem applies, so complete formulation requires an evolution equation for
the temperature field. The heat equation in orthogonal curvilinear coordinates is


∂T
∂ 1
√
κ ggij j + 4ηde 2 .
(2.74)
ρc∂i T =
√
i
∂x g
∂x
Note that ∂i means sum over four vectors and de is the invariant. Sincehorizontal diffusion is smaller than other parts,





∂T
uα ∂T
uβ ∂T
1
∂
√ ∂T
γ ∂T
ρc
+
+
+u
=√
κ g
+ 4ηde 2
∂t
g1 ∂α
g2 ∂β
∂γ
g ∂γ
∂γ
(2.75)[Greve and Blatter, 2009], where c is specific heat and κ is heat conductivity.
Since the ice temperature must not exceed the pressure melting point, thesolution of the equation is subject to the secondary condition T < Tm .

2.9

Boundary conditions

Ice sheets and glaciers do not only change centrally. Changes at their boundaries are also important processes, including basal sliding, accumulation andablation, and calving front processes. These processes are summarized here.

2.9.1

Basal sliding

There are three components for ice sheet and glacier motion. These areplastic deformation of the ice, sliding of ice over its bed, and deformationof the bed itself. The basal motion of the ice consists of sliding of ice anddeformation of the bed. They are the important components of ice flow.
The bed can deform significantly only if it consists of sediments saturated
with water at a pressure close to that exerted by the overlying ice. Sliding
is expected where the basal ice is close to its melting point. Sliding and
bed deformation are not mutually exclusive, displacement can occur at theinterface between ice and a deforming bed. Indeed, this interface may be
hard to define when the basal ice contains a lot of debris.
The basal processes under glaciers are still under discussion due to thedifficulties of measuring them. Therefore, improving understanding of themis one of the current major problems in glacier physics. Some theories on thesebasal processes have been elaborated based on the available, but limited,36

data and on physical considerations. The “sliding law” is the relationshipbetween basal velocity, shear and normal stresses, and water pressure and thecharacteristics of the glacier’s bed. It is an important boundary condition
for the analysis of flow in glaciers.
Two mechanisms are considered for ice flow over a rough, hard bed(Weertman, 1957). Ice can move over bed obstacles by the same processof creep deformation which gives rise to ice flux when there is no sliding. Alternatively, it can melt upstream of a bed obstacle due to increased pressure
at that point. If refreezing is ignored and the ice is everywhere in contact
with the bed, Fowler’s sliding law is
ub = Cb τbn ,(2.76)

where ub is the basal velocity and τb is basal shear stress. The constant, Cb ,depends on the geometry of bed roughness and n is the exponent of Glen’sflow law.
Since hydrostatic pressure and water pressure affect basal flow, a theory was elaborated incorporating these components (e.g., Hindmarsh and Le
Meur 2001). The basal sliding law used in this study is the sub-melt slidinglawτp0(2.77)
ub = Cb eTb γ bq .
Nb
The sliding law has been modified to allow for sub-melt sliding by adding
the coefficient γ = 1◦ C. T 0 b is the temperature relative to pressure melting,p and q are exponents; p = 3 and q = 2 are used in this study. There are
two options for the effective pressure N . One isN = ρgH,(2.78)

where H is local ice thickness. This option considers only hydrostatic pressureof ice. Another is
N = ρgH − pw .
(2.79)
This option considers the reduction of the effective pressure through waterpressure.

2.9.2

Calving front

The calving front of an ice shelf faces the ocean. Therefore, it is affectedby the ocean’s hydrostatic pressure (Fig. 2.5). The hydrostatic pressure
distribution of sea water provides a dynamic boundary condition for thecalving front and this is(t|cf )j i gj = −psw gi ,(2.80)
37

where psw is
psw = 0,for γ > γsw= ρsw G(γsw − γ) for γ ≤ γsw(2.81)

With the hydrostatic approximation and vertical integration, this gives
ρ GH 2)gα ,ρsw 2
ρ GH 2)gβ .
(2N β β + N α α )cf gβ + N α β |cf gα = ρ(1 −ρsw 2(2N α α + N β β )cf gα + N β α |cf gβ = ρ(1 −(2.82)

The calving front condition with Cartesian coordinates is


ρ
gH 2
∂vx
∂vy
∂vx ∂vy4η
ny =(ρsw − ρ)|cf nx + 2η
|cf nx + η
+
nx ,
∂x
∂y
∂y
∂x cf
ρsw
2


ρ
∂vy
∂vx
∂vx ∂vy
gH 2
4η
nx =
|cf ny + 2η
|cf ny + η+(ρsw − ρ)
ny .
∂y
∂x
∂y
∂x cf
ρsw
2
(2.83)The evolution of the calving front is the other condition. Nye’s calving
formulation is sometimes used in glacier models. But it may not apply to ice
shelves. Although several applications of the dynamic calving law have been
proposed (e.g., Alley et al. 2008), it is still under construction. Because of
the lack of a better criterion, it has been assumed in some modelling studies
that the ice shelf ends at the place where the ice thickness reaches a certain
value, such as 250 m (Paterson 1994).

2.9.3

Surface temperature and mass balance

For temperature evolution, a further thermodynamic boundary condition
must be given. This can be done simply by specifying the surface temperature
Ts ,T = Ts .
(2.84)
Since the temperature of ice cannot exceed its melting point, a finite temperature below melting point is applied.
Surface mass balance, as is given as an input from the atmosphere. If thesurface temperature is above melting point, melting is considered.

38

2.9.4

Basal temperature and sub-marine mass balance

The base of an ice shelf rests in water by definition. At the interface between
liquid water and solid ice, it should be on the pressure melting point of theocean water. The basal temperature at the ice shelf base, Tb , is given as
Tb = Tm − 2◦ C.

(2.85)The mean melting point depression of sea water, 2◦ C, is applied. This isapplied for the ice base as the boundary condition for temperature.
Sub-ice-shelf melting, mm , is parameterized by a method inherited fromthe GRISLI model (Ritz et al., 2001),
 (mm )gl in the vicinity of the grounding line,(mm )cs over the continental shelf, zl > −1600 m ,mm =(2.86)
(mm )ao over the abyssal ocean, zl ≤ −1600 m .
Standard values were used for the three parameters Mgl = 1 m ice equiv. a−1 ,Mcs = 0.1 m ice equiv. a−1 and Mao = 10 m ice equiv. a−1 .

2.10

Topographic evolution

Ice sheet changes its shape by accumulation, ablation at the base and thesurface, and by ice flux. The vertical integral of the continuity equation
then gives the height function of the ice sheet. Two different approaches are
introduced here; the ice surface equation and the ice thickness equation.

2.10.1

Ice thickness equation

The ice thickness equation is∂H= ∇ · Q + as − ab ,(2.87)∂twhere Q is the horizontal ice flux (the vertical integral of horizontal velocity),
as is surface mass balance and ab is basal mass balance.
Since the horizontal velocity is vertically constant at the ice shelf, theflux isQx = uα H, Qy = uβ H.
Therefore, the ice thickness equation is∂H= −Hui ;i + as − ab ,
∂t
 √

√
H ∂ guα ∂ guβ= √++ as − ab .
g
∂α
∂β
39

(2.88)(2.89)(2.90)

2.10.2

Ice surface equation

The ice thickness equation can be modified,∂H= −divQ + as − ab ,∂td(h − b)= divQ + as − ab ,dt

(2.91)

where h is the ice surface and b is the ice base. The elevation of the ice baseis equal to the elevation of the lithospheric surface below the ice sheet. This
is not the case for ice shelves because these float on the ocean. However, if
sea level zsl and ice surface are known, it is possible to estimate the elevationof the ice base using the floating condition
b=−

ρsw
ρ
h+
zsl .
ρsw − ρρsw − ρ(2.92)

Therefore, the ice surface evolution equation in ice shelves is
∂h
ρ
∂zsl
ρ
= −(1 −)∇ · Q ++ (1 −)(as − ab ).
∂t
ρsw
∂t
ρsw

(2.93)

If the horizontal difference of the sea level is neglected,
∂h
ρ
ρ
∂zsl
= −(1 −)∇ · Q + (1 −)(as − ab ) +
.
∂t
ρsw
ρsw
∂t

(2.94)

Note that ∇ · Q = ∇ · v(h − b) when the sea level is horizontally constant.
There is a problem around the grounding line when using the ice surface
evolution equation. This is because, although it is valid for ice shelves, it is
not clear whether the ice at the transition zone should be treated as ice sheet
or ice shelf. It is therefore difficult to choose the correct calculation of massbalance around grounding lines.

40

Figure 2.1: The typical state of an ice sheet. Arrows show driving forces
and flow of ice, and the resisting force working on the bedrock [Greve andBlatter, 2009].

41

Figure 2.2: The sliding motion of the ice sheet. Where b is Burger’s vector.
d is the distance. σ is the internal force. The defect moves through theA-B line. The three red lines become two straight lines, then molecular
movements [modified from Hondoh, 2000].

42Figure 2.3: Flow regimes in an ice sheet. Simple, bed-parallel shear flow
prevails except in the ice dome and margins [Greve and Blatter, 2009].

43

Figure 2.4: Flow regime at the ice sheet to ice shelf transition. In the
grounded ice sheet shear flow prevails, whereas in the floating ice shelf, plug
flow is present [Greve and Blatter, 2009].

44

Figure 2.5: The condition of the calving front. Hydrostatic pressure acts
below sea level [Greve and Blatter, 2009].

45

Chapter 3
Numerical methods
3.1

Ice sheet model SICOPOLIS

The ice sheet model SICOPOLIS is a 3D dynamic-thermodynamic ice sheet
model constructed by Greve [1995]. The purpose of the model is to simulate
large-scale ice sheet evolution with thermodynamical coupling of the ice sheet
and the lithosphere. The components of the model are the dynamics of
ice and of the lithosphere, the atmospheric (temperature and accumulation
rates) and the oceanological boundary conditions (mean sea level and subice-shelf melt) (Fig. 3.1).
A shallow ice approximation with a Weertman-type sliding law is used
for the dynamics. There are several options for the constitutive equation
and Glen’s flow law is one of them. It is also possible to use a polythermal
approach. Isostatic displacement is also applied to the lithosphere dynamics.
The thermodynamics of the ice and of the lithosphere are coupled. The finite
difference scheme with the Arakawa C grid [Arakawa and Lamb, 1977] is used
for the field equations (Fig. 3.2). Scalar variables such as ice thickness andtemperature are on the grid point and the velocities are on the staggered grid
point. For the vertical layers, a sigma coordinate system, with coordinatesnormalized by local thickness, is applied.

3.2

Details of the numerical grid

The surface slope is a significant factor affecting ice flow. If the discretizations
are done in a straightforward way, the exact positions of the ice surface are
missed. Ice flow can then be over- or underestimated. Sigma coordinates areapplied to the ice sheet model to represent precisely the positions of the ice
surface and ice base (Fig. 3.3). A terrain-following coordinate transformation46

is applied which maps the local ice thickness on to unity,ξ = x,
ϕ = y,z − b(x, y, t)ζ=H(x, y, t)
τ = t.

(3.1)
(3.2)

The (x, y, z, t) coordinate system is mapped to (ξ, ϕ, ζ,t au) coordinates. Thistransformation, often referred to as the sigma transformation, maps the ice
surface h = h(x, y, t) to ζ = 1 and the ice base b = b(x, y, t) to ζ = 0. In thetransformed domain, a regular, rectangular grid with spacings δξ, δϕ and σζ
can easily be defined, such that the uppermost layer of grid points matches
the ice surface and the lowermost layer the ice base. It is not enough justto transform the coordinates as it is also necessary to transform the field
equations,
∂
∂x
∂
∂y
∂
∂z
∂
∂t

∂
∂ζ ∂
+
∂ξ ∂x ∂ζ∂∂ζ ∂=+∂ϕ ∂y ∂ζ
∂ζ ∂
=
,
∂z ∂ζ
∂
∂ζ ∂
=
+
.
∂τ
∂t ∂ζ
=

(3.3)

The respective derivatives are not the same.
Another treatment of the grid condition is also applied to the model. It
is called the Arakawa C grid a type of staggered grid system into which
numerical values can be inserted (Fig. 3.2). Scalar field variables are definedon the regular grid (ice thickness, temperature, lithosphere surface, ice surface and so on), but velocities are inserted at the mid-point of the two grid
points. First, we define the regular grid:ξi = ξ0 + i∆xi,ηj = η0 + j∆y,(ζc )kc = kc /kc,max ,(ζt )kt = kt /kt,max ,(ζr )kr = kr /kr,max ,tn = t0 + n∆t,

(3.4)

(3.5)
47

The velocity components vx , vy and vz are inserted at the mid-point between
two points on each axis.

3.3

Coupling of the ice shelf model to the ice
sheet model

It is necessary to couple the ice shelf model with the ice sheet model tounderstand the dynamics of the ice sheet while taking the ice shelves into
consideration. The ice shelf dynamics model is as follows:
Three different time steps are applied to the ice sheet model. The longest
time step, ∆tl , is used in the dynamics of the lithosphere. The second, longer,
step is used in the temperature evolution equation. The shortest time stepis taken for the ice thickness equation.

3.4

Shallow shelf dynamic equations

The shallow shelf approximation is discretized with the staggered Arakawa
C grid. The basic idea is the discretization with respect to the staggered grid
point (i + 1/2, j).

3.4.1

Field equations

If the staggered grid point (i + 1/2, j) has two ice shelf neighbours, (i, j) and(i + 1, j), then it is considered to be an ice shelf. The shallow shelf equation

48

for vx is discretized to




4
ηi+1,j (vx )i+3/2,j − (vx )i+1/2,j − ηi,j (vx )i+1/2,j − (vx )i−1/2,j
2
(∆x)
+

+

+

=





2
ηi+1,j (vy )i+1,j+1/2 − (vy )i+1,j−1/2 − ηi,j (vy )i,j+1/2 − (vy )i,j−1/2
∆x∆y


1
ηi+1/2,j+1/2 (vy )i+1,j+1/2 − (vy )i,j+1/2
∆x∆y


−ηi+1/2,h−1/2 (vy )i+1,j−1/2 − (vy )i,j−1/2


1
η
(v
)
−
(v
)
x
x
i+1/2,j+1/2
i+1/2,ji+1/2,j(∆y)2−ηi+1/2,j−1/2 (vx )i+1/2,j − (vx )i+1/2,j−1ρgHi+1/2,j (hi+1,j − hi,j ).
∆x

(3.6)
Which is applied for vy when both neighbors of the staggered grid point(i, j + 1/2), (i, j) and (i, j + 1), are ice shelves,




4
η
(v
)
−
(v
)
−
η
(v
)
−
(v
)
i,j+1
y
y
i,j
y
y
i,j+3/2
i,j+1/2
i,j+1/2
i,j−1/2
∆y 2
+

+





2
ηi,j+1 (vx )i+1/2,j+1 − (vx )i−1/2,j+1 − ηi,j (vy )i+1/2,j − (vy )i−1/2,j
∆x∆y


1
ηi+1/2,j+1/2 (vx )i+1/2,j+1 − (vx )i+1/2,j
∆x∆y


−ηi−1/2,h+1/2 (vx )i−1/2,j+1 − (vx )i−1/2,j

+



1
η
(v
)
−
(v
)
y
y
i+1/2,j+1/2
i,j+1/2
i+1/2,j∆x2
−ηi−1/2,j+1/2 (vy )i,j+1/2 − (vy )i−1,j+1/2=

ρgHi,j+1/2 (hi,j+1 − hi,j ).
∆y
(3.7)

3.4.2

Calving front conditions

The shallow shelf approximation for the calving front is applied if the grid ison the ice sheet but the next point is the ocean. The calving front condition49

is applied, 2ηi,j 4ηi,j (vx )i+1/2,j − (vx )i−1/2,j +(vy )i,j+1/2 − (vy )i,j−1/2
∆x
∆x
ρg
2
=
(ρsw − ρ)Hi,j.
2ρsw
 2ηi,j 4ηi,j (vy )i,j+1/2 − (vy )i,j−1/2 +(vx )i+1/2,j − (vx )i−1/2,j
∆y
∆x
ρg
2
(ρsw − ρ)Hi,j
.
=
2ρsw

(3.8)

It is presumed that ice fills the whole of a grid and does not leak; thecalving front is perpendicular to each axis.

3.4.3

Grounding zone

Since the ice sheet model uses a staggered grid (Arakawa C grid), there is no
ice thickness information at the points where velocity is defined. There areseveral possible ways of treating this condition.
One way is for the mid point of the ice sheet and ice shelf to be treatedas an ice shelf. However, this causes a large ice flux into the sea because thevelocity calculated for ice shelves is greater than that calculated for ice sheets.
Moreover, because the slope for the ice sheet is steeper than for the ice shelf,
the velocity simulated with the shallow shelf approximation becomes large.
The rate of grounding line migration then increases. An alternative is for the
point to be treated as an ice sheet. It saves ice shelves from the flux above
but the rate of grounding line migration can be too slow. The third possible
way of treating this condition is to check the conditions of the middle point.
The floating conditions of the staggered grid are linearly estimated by usingneighboring ice thickness data,
Hi+1,j + Hi,j,2
Hi,j+1 + Hi,j
Hi,j+1/2 =.
2

Hi+1/2,j =(3.9)

If the estimated ice thickness is lower than the floating thickness, then it istreated as an ice shelf. Otherwise, it is treated as an ice sheet. Then, one
point further downstream on the grid is the ice shelf.
Pollard and DeConto [2009] compared the simulated flux of the middle
point with the flux estimated from 2D boundary layer theory. If the flux
estimated by the boundary layer theory is larger than the flux calculatedby the large scale ice shelf model, the solution of boundary layer theory is50

applied to that point. Otherwise, the boundary layer flux is imposed on the
point one step downstream. Such treatment can improve the fit around thegrounding zone.
It is difficult to define the basal mass balance at the grounding line, thegrid next to an ice shelf. Although the grounding line point is not floating,it is affected by the sea since it is next to an ice shelf. The simple way todetermine melt rates is just to use the ocean-induced melting when the baseof the ice is lower than sea level,
ab = mm .

(3.10)

It is a good approximation when ocean-induced mass balance is dominant.
In another way, the basal mass balance depends on the number of ice
shelf domains next to the grounding line,
rg =

Ngr
,
Nall

(3.11)

where rg is the ratio of the grounded grid number, Ngr , to the total of grids
next to the grounding line, Nall . The total basal mass balance is
ab = rg mm + (1 − rg )mg ,(3.12)

where mg is melt rate when the point is grounded, estimated by the geothermal flux and strain heating of ice.
There is another option which depends on depth. In reality, the grounding
line itself meets the floating condition. However, the numerical grounding
line point can be deeper or shallower than the ice sheet. If the elevation ofits base is not deep enough to satisfy the floating condition, the effect of theocean is limited. The ocean-induced basal mass balance is determined by its
depth relative to the floating condition;• when the ice base is above sea level,
ab = mg ,(3.13)• when the base is below sea level,
bgl = mm (

H ρsw
H ρsw
) + mg (1 −).
zsl − b ρ
zsl − b ρ

51

(3.14)

3.5

Temperature evolution

The temperature evolution equation is not very different from that used forthe ice sheet. The differences are the strain heating and the basal boundary
condition.
The strain heating term is calculated where viscosity is calculated. Therefore, the value is implemented for the temperature evolution equation. For
the basal boundary condition, a temperature two degrees above the melting
point is applied as the boundary condition.

3.6

Topography evolution

The margin of the glacier changes through flow and accumulation. There areseveral ways to manage ice sheet topography evolution. One of the traditional
ways is the ice surface equation and another is the ice thickness equation.

3.6.1

Ice surface equation

The ice surface evolution equation was applied initially. However, it may notconserve mass balance in the vicinity of the grounding lines. The neighborof a grid point at the grounding line is on the ice shelf. The ice surface
equation for the ice sheet is used for the grid and the equation for the ice
shelf is applied to the neighbors. But it is not known whether the grid will
still be ice sheet or ice shelf in the next time step.
The ice surface equation is




n
hn+1
∂h
∂h
i,j − hi,j
= −
vx
− vy
∆t
∂x up
∂y up
(vx )i+1/2,j − (vx )i−1/2,j∆x(v)− (vy )i,j−1/2
yi,j+1/2n− Hi,j∆yn+1bi,j − bni,j+ as − ab +.
∆tn
− Hi,j(3.15)

where the upstream term for the staggered grid has two parts. One example
for the x direction is






hi+1,j − hi,jhi,j − hi−1,j∂h= − (vx )i+1/2,j+ (vx )i−1/2,j
.
vx
∂x up
∆x
∆x
up
up
(3.16)
52

Then, the right hand side is
hi+1,j − hi,j(vx )i+1/2,j
=
∆x
up



hi,j − hi−1,j(vx )i−1/2,j
=
∆x
up

3.6.2

(vx )i+1/2,j0.

(vx )i−1/2,j0

hi+1,j − hi,j∆x

if (vx )i+1/2,j < 0 ,hi,j − hi−1,j∆x

if (vx )i+1/2,j >= 0 ,if (vx )i+1/2,j >= 0 .
(3.17)

if (vx )i+1/2,j < 0 .
(3.18)

Ice thickness equation

As stated above, the ice thickness equation is good for estimating mass conservation. Therefore, attempts have been made to apply it to ice shelves aswell. The ice thickness equation is




n+1
n
− Hi,j
Hi,j
∂H
∂H
− vy
= −
vx
∆t
∂x up
∂y up
(vx )i+1/2,j − (vx )i−1/2,j∆x(v)− (vy )i,j−1/2
yi,j+1/2n− Hi,j+ as − ab .
∆yn
− Hi,jThe upstream term is the same as the ice surface equation.

53(3.19)

Figure 3.1: Components of the ice sheet model, SICOPOLIS.

54

Figure 3.2: The Arakawa C grid is a numerical grid system defined by
Arakawa and Lamb [1977].

55

Figure 3.3: The sigma coordinate system of the ice sheet [Greve and Blatter,2009].

56

Chapter 4
Numerical tests
4.1

Ice shelf ramp

The performance of the shallow shelf model was tested for the case of an ice
shelf ramp for which an analytical solution exists [Greve and Blatter, 2009].
The analytical solution is given in two-dimensional flow in the x − z plane,
where x is the direction of flow and z is positive upward (Fig. 4.1). The
length of the shelf is L. The grounding line lies at x = 0, the calving front
lies at x = L. Ice thickness declines linearly from H = Hgl at the grounding
line to H = Hcf at the calving front,
H = Hgl −Hgl − Hcf
x.
L

(4.1)

With the floating condition, the ice surface and the ice base are given byh=ρsw − ρρH, b = −
H.
ρsw
ρsw

(4.2)

The rate factor of the constitutive equation is a constant, A(T, p) = A =const. The analytical solution for the ice shelf ramp is ρg n Z x
A
H n (x)dx,(4.3)
vx (x) = vgl +40where n is the exponent of Glen’s flow law.
Since the ice shelf model has both x-y horizontal dimensions, the directionneeds to be set perpendicular to the flow (y direction). In the ice shelf ramp
experiment, ice thickness is constant through y direction, i.e. (∂H/∂y = 0)
and the width was set wide enough to eliminate the effect of the side wall.
Finally, it has a rectangular domain with x (the flow direction) and y (the57transverse direction). The solution through the central flow line will be thesame as the original two-dimensional setting. Mean sea level is z = zsl = 0.
Incoming ice flow velocity at the grounding line is vgl = vgl ex , i.e., theincoming ice flux of the shelf is vx Hgl . The velocity at the side wall (y = 0
and y = ymax) is zero. The ice shelf is set at a steady state, ∂/∂t = 0,
for all variables. The input velocity at the grounding line is vgl = 100 m/a.
Ice thickness at the grounding line and the calving front is Hgl = 400 m,
Hcf = 200 m. The length of the shelf is L = 200 km. The rate factor is
A = 4.9 × 10−25 s−1 P a−3 .
The velocity map of the experiment shows that the field is symmetricalabout the central flow line and is fastest at the central flow line, where theeffect of the side wall is negligible (Fig. 4.2). The velocity through the central
flow line nicely reproduced the analytical solution (Fig. 4.3).

4.2

Shelfy stream experiments

The performance of the shelfy stream model was also tested. The rectangulardomain is also used in this test. The length of the ice shelf is Lsh = 100 km
in the x direction. The length of the stream upstream from the ice shelf is
Lst = 50 km (Fig. 4.4a). The ice thickness is Hin = 400 m at the initial point,
Hgl = 300 m at the grounding line and Hcf = 200 m at the calving front. The
ice thickness profile isHin − Hgl(x ≤ 50 km)
Lst
Hgl − Hcf
H(x)sh = Hgl − x(50 < x ≤ 150 km).
Lsf
H(x)st = Hin − x(4.4)

The basal flow law is set the same as Schoof [2007a]. The result indicatesthat the velocity is relatively fast, more than 1 km/a, even around grounded
ice stream areas. The velocity is accelerated to 3 km/a near the calving front.

58

Figure 4.1: The profile of the ice shelf ramp [Greve and Blatter, 2009].

59

Figure 4.2: The condition of the 3D ice shelf ramp.

60

Figure 4.3: The result for the central flow line of the 3D ice shelf ramp.

61

a

b

Mask

Velocity

Figure 4.4: The setting of the ice stream region and the ice shelf region andthe result of the shelfy stream run.

62

Chapter 5
SeaRISE experiments
5.1

SeaRISEAs stated in Chapter 1, it is unclear how much ice sheet instabilities contribute to sea level rise. It is therefore necessary to quantify how sensitive
ice sheets are to several causes of instability.
Sea-level Response to Ice Sheet Evolution (SeaRISE) is a communityorganized effort to estimate the likely range of ice sheet contributions to sea
level change over the next few hundred years. Common datasets for theice sheet study have been developed in the project. The numerical experiments are conducted by using SeaRISE datasets (Fig. 5.1). Several types ofexperiments are being designed within the project.

5.2

Application of the model to the SeaRISE
set-upAfter the model is validated, it is coupled to the ice sheet model. The model isapplied to the Antarctic ice sheet with the SeaRISE dataset. The geometryof the experiment is developed by the SeaRISE project; ice thickness, thebedrock topography, accumulation rates, ice surface and the mask of floating
ice and grounded ice. The parametrization of current surface temperature istaken from Fortuin and Oerlemans [1990],
 7.405 − 0.014285 h − 0.180 φ , for h > 1500 m ,36.689 − 0.005102 h − 0.725 φ , for 200 m < h ≤ 1500 m ,(5.1)
Tma =
49.642 − 0.943 φ ,for h ≤ 200 m ,Tms = 16.81 − 0.00692 h − 0.27937 φ ,
63

(5.2)

where temperatures are in ◦ C, surface elevations in m AMSL and latitudesin ◦ S (counted positive).
The calving criterion thickness is applied in this study. The thresholdvalue is 200 m in some simulations. There are considerable areas of ice shelfthat are thinner than this value (Fig. 5.2). The Riiser-Larsen and Shackleton
ice shelves, in particular, lose large areas with a larger threshold of 200250 m than many of the smaller ice shelves. Moreover, using a large valueof threshold thickness prevents a new ice shelf grid being made. Therefore,
another threshold value, 50 m, was used in most runs.

5.3

Paleoclimatic spin-up

In order to obtain a suitable present-day configuration of the Antarctic ice
sheet, it is desirable to carry out a paleoclimatic spin-up over at least a full
glacial cycle. The simulated present-day configuration of the ice sheet isused as the initial conditions for simulation of the effects of future climate
described below. Initial testing revealed that it is very difficult to reproduce
the observed geometry (in particular the distribution of grounded vs. floatingice) by an unconstrained, freely evolving simulation. For this reason, we
carried out the spin-up simulation in three steps (Fig. 5.3): (i) an initial
relaxation run with freely evolving grounded ice topography over 100 years(starting from the present-day geometry and isothermal conditions at −10◦ C
everywhere) in order to avoid spurious noise in the computed velocity field[Calov, 1994]; (ii) a steady-state run from 250 ka BP to 125 ka BP withthe topography kept fixed over time and the climatic forcing kept steady at
125 ka BP conditions; (iii) a transient run from 125 ka BP to today withthe topography kept fixed over time in order to enforce a good fit between
the simulated and observed present-day topographies [see also Seddik et al.,2012].
The mean annual and mean summer surface temperatures (Tma and Tms ,respectively) are decomposed into present-day spatial distributions plus apurely time-dependent anomaly ∆T (t),present
Tma (h, φ, t) = Tma(h, φ) + ∆T (t) ,present
Tms (h, φ, t) = Tms(h, φ) + ∆T (t) ,(5.3)

where t is time and φ geographical latitude. The time-dependent anomaly∆T (t) results from the Vostok deuterium isotope ratio record converted totemperature with the relation of Petit et al. [1999]. Precipitation, surface
melting and sea level forcing are irrelevant here due to the fixed topography
approach. The paleoclimatic spin-up is conducted with such conditions.
64

The SeaRISE Antarctica Developmental Data Set provides two different
distributions of geothermal heat flux by Shapiro and Ritzwoller [2004] andFox Maule et al. [2005] (Fig. 5.4). We conducted the spin-up run with bothdatasets with heat flux values capped at 70 mW m−2 following a recommendation from the SeaRISE project.
Although several different resolutions were tried, the results at 20 km
resolution are discussed here. The basal temperature results (Fig. 5.5a) areat the pressure melting point in large areas in West Antarctica, while limited
temperature base areas are found in East Antarctica, in George V Land,
Wilkes Land and several other places near the coast. At the deep ice coresites Vostok, Dome C, Dome F, Kohnen and Byrd, the spin-up run produces
basal temperatures of −11.4◦ C, −11.6◦ C, −18.4◦ C, −11.3◦ C and −2.0◦ C(pressure melting point), respectively. Except for Byrd, these temperaturesare significantly below the pressure melting point, while the observed temperatures all reach the pressure melting point. This discrepancy is probably
due to the shortcomings of the Shapiro and Ritzwoller [2004], geothermalheat flux distribution that was applied which gives values in the range of
45 – 50 mW m−2 at Dome C, Dome F and Vostok. This is probably not sufficient to form a temperate base [Llubes et al., 2006]. Seddik et al. [2011] useda higher geothermal heat flux of 60 mW m−2 to produce a realistic boreholetemperature in the vicinity of Dome F. Therefore, larger geothermal heat
fluxes would be needed to reproduce the melting conditions that have been
observed in East Antarctica. By contrast, the base of the West-Antarctic
Byrd ice core, which receives a geothermal heat flux of 70 mW m−2 in our
simulation, reaches the melting point. The basal temperature is high at EAIS
in the Fox Maule et al. [2005] spin-up. However, the WAIS is colder in theseexperiments.

5.4

Sensitivities of the ice volume to different
datasets

Since there are two different sets of input data of geothermal heat flux andaccumulation rates (Fig. 5.4, Fig. 5.6), it is preferable to estimate the effectof using these different datasets. The sensitivities to the input dataset were
estimated. There are four different combinations,
1. Shapiro and Ritzwoller [2004] geothermal heat flux (SR) and Arthern
et al. [2006] (AE) accumulation rate
2. SR geothermal heat flux and Van de Berg et al. [2006] (BE) accumulation rate65

3. Fox Maule et al. [2005] geothermal heat flux (FE) and AE accumulationrate
4. FE geothermal heat flux and BE accumulation rate.
Two different paleoclimatic spin-ups were carried out to accommodate thetwo different geothermal flux datasets. The four experiments were conductedfrom the end of spin-up and calculated over 500 years. Figure 5.7 showsthe evolution of the grounded ice volume for the future control run withthe two different sets of accumulation data (Section 5.5). In both cases,
the ice volume increases by more than 1.5 m SLE over 500 years, and theaccumulation data from Van de Berg et al. [2006] produce ∼ 35 cm SLE
more ice than those from Arthern et al. [2006].
The reason is that the former accumulation rates are relatively larger thanthe latter, particularly around the margin of ice sheets such as the Antarctic
Peninsula and the margins of East Antarctica. By contrast, different geothermal flux settings (Shapiro and Ritzwoller [2004] vs. Fox Maule et al. [2005])
cause only slight changes. It was finally decided to use the SR geothermal
flux and AE accumulation rates for the control run.

5.5

Future climate experiments

The results of the paleoclimatic spin-up run for the present are shown in
Fig. 5.8. Due to the fixed-topography approach, the surface topography isof course very similar to that observed (not shown). Small differences ofthe order of tens of meters arise as a consequence of the initial 100-year
relaxation. The surface velocity shows the expected distribution with small
velocities (< 10 m a−1 ) in the interior, a general speed-up towards the coastand the largest velocities of 1000 m a−1 and more for the ice shelves. The
simulated velocity patterns nicely reproduce the observed velocity pattern.
The experiment on future climate was started from the end of spin-up(2004) and simulated hundreds of years with varied boundary conditions.
For the future climate experiments we used the following set of SeaRISE
experiments:• Experiment C1: constant climate control run; beginning at present(more precisely, 2004-1-1 0:00, corresponding to t = 0) and runningfor 500 a holding the climate steady at present climate values. We
carried out this run alternately with the surface accumulation data from
Arthern et al. [2006] and Shapiro and Ritzwoller [2004] geothermal flux
(as provided in the Antarctica Developmental Data Set by SeaRISE).
66

• Experiment C2: AR4 climate forcing running for 500 a.
• Experiment C3: 1.5 times AR4 climate forcing running for 500 a.
• Experiment C4: 2.0 times AR4 climate forcing for 500 a.
• Experiment M0: constant climate forcing running for 500 a, AE surface
accumulation rates, zero sub-ice-shelf melting rate.
• Experiment M1: constant climate forcing running for 500 a, AE surface
accumulation rates, uniform sub-ice-shelf melting rate of 2 m ice equiv. a−1 .
• Experiment M2: constant climate forcing running for 500 a, AE surface
accumulation rates, uniform sub-ice-shelf melting rate of 20 m ice equiv. a−1 .
• Experiment M3: constant climate forcing running for 500 a, AE surface
accumulation rates. uniform sub-ice-shelf melting rate of 200 m ice equiv. a−1 .
• Experiment S1: the basal sliding velocity is doubled. Constant climate
forcing running for 500 a, AE surface accumulation rates.
• Experiment S2: the basal sliding velocity is increased by a factor of 2.5.
Constant climate forcing running for 500 a, AE surface accumulation
rates.
• Experiment S3: the basal sliding velocity is increased by a factor of 3.
Constant climate forcing running for 500 a, AE surface accumulation
rates.
• Experiment T1: the basal sliding velocity is doubled. AR4 climate
forcing running running for 500 a, AE surface accumulation rates. A
uniform sub-ice-shelf melting rate of 20 m ice equiv. a−1 .
For all experiments, the geothermal heat flux is that of Shapiro and Ritzwoller[2004], capped at 70 mW m−2 (the same as used for the paleoclimatic spin-up,
see above).

5.6

Sensitivities on varied sub-ice-shelf melting rates

The model was applied to the diagnostic condition. The resolution is 40 km.
The number of nodes is 150 × 150 for the horizontal directions, and 80 for the67

vertical direction. The total number is 150 = 1, 800, 000. The ice temperature was set to −10◦ C. The results show that inland ice volume decreases
when basal melt rates increase (see Fig. 5.9). The basal pressure reduction bybasal water pressure was neglected. This is caused by the reduction of buttressing from ice shelves, as stated above. On the other hand, if ice shelves
disappear, land ice volume subsequently increases. The reason for increased
ice volume is changes in boundary conditions. In this simulation, melting
around the marine ice margin is not included. After ice shelf removal due toenhanced melting, ocean induced melting disappears from the last grounding
line. Therefore, volume loss is smaller than when ice shelves exist. Moreover,velocities in ice shelves are higher than those in the ice sheet since there isno friction under ice shelves. The disappearance of ice shelves reduces theflux of ice away from the land ice. This effect and reduction of buttressing
causes the changes in land ice volume.
Another experiment was carried out with melting marine ice margins.
The basal pressure reduction by basal water pressure is also included (Fig. 5.10).
In this case, ice volume of the M3 experiment is smaller than in M2. But thedifference is less than 0.1 mm sea level equivalent in 100 years and less than
0.2 mm sea level equivalent in 500 years. It suggests there is a limitation ofice volume change by ocean induced melting.
The results of simulations of the effect of future climate change of M0-M3
setting in 20 km resolution are shown in Fig. 5.11. As expected, grounded icevolume declines when basal melting rates are higher. However, the effect isstrongly non-linear and diminishes with increasing melting rates. This agrees
with the findings of Huybrechts and de Wolde [1999] for a similar series ofsimulations with smaller basal melting rates of 0, 1, 3, 5 and 10 m a−1 . The
range of ice volume changes across all experiments is ∼ 0.5 m SLE after 100
years and ∼ 1.5 m SLE after 500 years. We thus see a significant, but notcatastrophic, impact of increased basal melting (up to 200 m a−1 !) on ice
sheet decay.
The changes in the area of grounded and floating ice are shown in Fig. 5.12.
In the control run (basal melting according to Eq. (2.86)), the grounded ice
area grows rapidly by ∼ 0.2 × 106 km2 at the expense of floating ice, andremains relatively constant after that. This increased area of grounded ice
also explains the increase in grounded ice volume discussed above. In theextreme run M3 (200 m a−1 basal melting), current floating ice disappears
almost entirely within a few years. The area of grounded ice remains largely
unaffected but, as seen above, the loss of ice shelf buttressing leads to a significant loss of ice volume compared to the control run. In the M2 experiment(20 m a−1 basal melting), most of the floating ice disappears within 100 years.
However, the remaining area of floating ice stabilizes at ∼ 0.2 × 106 km2 due
68to a balance between inflow from the ice sheet and basal melting.
In general, under the influence of high basal melting, the Ross and Antarctic Peninsula ice shelves decay faster than the Filchner-Ronne and Amery ice
shelves. In the M1 experiment (2 m a−1 basal melting), the Ross ice shelf loses
more than 50% of its area by t = 250 a and almost disappears by t = 500 a,while the Amery ice shelf does not lose much of its area. The vulnerability of the Ross ice shelf results from the fact that it is surrounded by theTransantarctic Mountains so that the incoming flux from the ice sheet issmaller than for other ice shelves. By contrast, the Amery ice shelf is largely
stabilized against increased basal melting by the large inflow from the Lambert Glacier, which is the world’s largest glacier. As expected, shrinkage ofthe WAIS was larger than that of the EAIS due to the bedrock topography(Fig. 5.13).
The grounded ice area is largest in run M0 (no basal melting), while thefloating ice area is largest in the control run (basal melting according to
Eq. (2.86)). This is related to grounding line migration. Since there is nobasal melting in run M0, the grounding lines can advance freely all aroundthe ice sheet, so that grounded ice grows at the expense of floating ice.
By contrast, basal melting favours grounding line retreat, so that floating
ice grows at the expense of grounded ice. The floating ice area of run M1(2 m a−1 basal melting) is larger than that of the control run during the first∼ 50 years, and larger than that of run M0 (no basal melting) during the first∼ 110 years. However, the retreat of the calving fronts due to basal-meltinduced thinning later outweighs this effect, so that the floating ice area ofrun M1is ultimately smaller than that of the control and M0 runs. For runs
M2 and M3 (20 / 200 m a−1 basal melting), basal melting is so large that theloss due to calving rapidly outweighs the gain due to grounding line retreat.

5.7

Results of whole set of SeaRISE experiments

The whole set of experiments was conducted in 10 km resolution. The ice
sheet volume is most sensitive to the M-series experiments (Fig. 5.14). The
next level of sensitivity is to the S-series experiments. The C-series experiments caused less change than the others. Ice volume decrease is largest in
T1. The total ice volume ranges from 1.3 m SLE after 100 years to 3.3 m
SLE after 500 years. In the M-series experiment, the difference in ice volume
is 0.7 m SLE after 100 years and 2.6 m SLE after 500 years. In the S-series
experiment, it is 0.4 m SLE after 100 years and 1.1 m SLE after 500 years.

69

In the C-series experiment, it is 0.15 m SLE after 100 years and 0.2 m SLE
after 500 years.
Ice shelf area is sensitive to basal melting (Fig. 5.15). Since changes
caused by the M1-M3 experiments are large, the timing of total ice volume
change is the same as for the change in ice shelf area. Ice shelf area reaches
a steady state after 10 years (M3) and 70 years (M2). A steady state is notreached within 500 years in the M1 experiment. The ice sheet volume of S2
is larger than S1 in 200-500 years (Fig. 5.16). The grounded ice area of S1 issmaller than the grounded ice sheet area of other experiments. The ice shelfarea increased in S1. The reason why ice volume in S1 is smaller than in S2
is that the loss of grounded ice area means loss of land ice regions. Therefore,
the total outgoing flux from the ice sheet is larger in S1 than S2. Althoughgrounded ice area decreased in the C-series experiment (Fig. 5.17), enhanced
accumulation through surface climate change increased its mass. The areaof grounded ice increased after 100 years due to grounding line migration.
However, ice volumes of C4 and C3 were smaller than in C1 at 500 years.
This was caused by increased surface melting around the ice sheet margins
in C3 and C4 compared to C1 experiments.

5.8

Relaxation of the initial shock

There are abrupt changes of ice shelf and ice sheet area in the first stepsof simulation, after the fixed topography spin-up (cf. Fig. 5.12). There areseveral reasons for this effect. Ice sheet evolution is allowed for the firststep of the spin-up. Although the grounding line and ice shelf thickness are
fixed, a certain amount of ice is dammed around the grounding line during
this period. This dammed ice immediately transfers to the grounding zone
after topographical evolution is allowed. It causes the abrupt grounding line
migration at the first step of the simulation. It changes rapidly in the firstyear. It also shows the fast response of the grounding line change. Although
the change becomes moderate after that, it can continue for decades.
Since such rapid change is not desirable, another process is added before
starting any runs. After the fixed topography spin-up, 20 years of ice sheet
and ice shelf topography evolution was added. The first abrupt changes of
ice sheet and ice shelf area were completed during this period. [Figure of ice
shelf, sheet area in 20 years added]. In this experiment, sub-ice-shelf melting
rates in the control run are doubled except for the outer ocean condition.
The results of the new control run and M-series experiments are shown in
Fig. 5.18. The initial abrupt change in floating ice area and grounded area is
not found in this simulation. The new runs also consider marine ice melting.
70

The comparison of marine ice melting with no marine melting experimentshows the importance of this setting. Since this process was not includedin the simulations, the Antarctic ice sheet does not lose enough mass bygrounding line migration. The sea level equivalent ice volume change isalmost 1.5 times larger (2.2 m SLE) than the experiment without marine ice
melting (1.5 m SLE). This agrees with the spread of 0.6 × 106 km3 (1.5 m
SLE) reported by Huybrechts and de Wolde [1999].
Changes in the area of grounded and floating ice are shown in Fig. 5.18.
In the control run, the grounded ice area drops by ∼ 0.25 × 106 km2 (∼ 2%)during the first 200 years, partly compensated by an increase in floating ice,and remains almost constant after that. In run M1 (2 m a−1 basal melting
everywhere), grounded ice evolves very similarly due to the identical melting
rates in the vicinity of the grounding line, whereas two thirds of the floating
ice disappears over 500 years as a consequence of the 10 times higher melting
rate over the continental shelf. For E1b most of it melts within 100 years,
the remaining area stabilizing at 0.15 × 106 km2 (10% of today’s value) due
to a balance between inflow from the ice sheet and basal melting, and forE1c almost all floating ice disappears within a few years.

71Figure 5.1: Ice thickness data used by the SeaRISE project.

72

3000

2000

y(km)

1000

0

−1000

−2000

−2000

−1000

0
1000
x(km)

2000

3000

Figure 5.2: The distribution of ice shelf regions thinner than 200 m. The
regions colored white are where the ice shelf is thinner than 200 m. The
black color is sea, the light gray is the ice sheet and the dark gray is the ice
shelf.

73

Figure 5.3: The temperature settings used in the paleoclimatic spin-up. The
deviation of temperature from 2004 is applied with present surface climate
data [Arthern et al., 2006] and geothermal heat flux [Shapiro and Ritzwoller,
2004].

74

a
˚W

0˚

30

30
˚E

Shapiro and Ritzwoller (2004)
50

˚S
50

3000

mW/m2
120
˚S

100
˚E
60

2000

70

1000

y (km)

80

60
0

90˚E

55
50

-1000

45
12

-2000

0˚E

-2000

-1000

0

1000

2000

40
30

3000

x (km)

b

30

˚W

0˚

30

˚E

Fox Maule et al. (2005)
50

˚S
50

3000

mW/m2
120
˚S

100
˚E
60

2000

70

1000

y (km)

80

60
0

90˚E

55
50

-1000

45
12

-2000

0˚E

-2000

-1000

0

1000

2000

3000

40
30

x (km)

Figure 5.4: Two geothermal flux dataset by (a) Shapiro and Ritzwoller [2004]
and (b) Fox Maule et al. [2005].

75

a

30

0˚

˚W
30

˚E

Shapiro and Ritzwoller (2004)
50

C

˚S

2000

60

˚S

50

0

3000

˚E

-4

y (km)

1000

0

-8

90˚E

-1000
-12
12

-2000

-2000

-1000

0

1000

2000

0˚E

-16

3000

x (km)

b

30

˚W

0˚

30

˚E

Fox Maule et al. (2005)
50

C

˚S

2000

60

˚S

50

0

3000

˚E

-4

y (km)

1000

0

90˚E

-8

-1000
-12
12

-2000

-2000

-1000

0

1000

2000

3000

0˚E

-16

x (km)

Figure 5.5: The results of paleoclimate spin-up using the two different
geothermal heat flux datasets in 20 km resolution (a) Shapiro and Ritzwoller[2004], SR, (b) Fox Maule et al. [2005], SE.

76

a

30

˚W

0˚

30

˚E

Arthern et al. [2006]
50

˚S

mm w.e./a
1000
50

3000

˚S

˚E
60

2000

800
600
500

y (km)

1000

400

0

90˚E

300
200

-1000

150
100
12

-2000

-2000

-1000

0

1000

2000

0˚E

75
50

3000

x (km)

b

30

˚W

0˚

30

˚E

Van de Berg et al. [2006]
50

˚S

mm w.e./a
1000
˚S
50

3000

˚E
60

2000

800
600
500

y (km)

1000

400

0

90˚E

300
200

-1000

150
100
12

-2000

-2000

-1000

0

1000

2000

3000

0˚E

75
50

x (km)

Figure 5.6: Two datasets of accumulation rates (a) Arthern et al. [2006] and(b) Van de Berg et al. [2006].

77

61.5

(mSLE)
SLE)
VV(m

61

60.5

AE FE
BE FE
AE SR
BE SR

60

59.5
0

100

200

yr

300

400

500

Figure 5.7: Simulated changes in grounded ice volume, V (expressed in
m SLE) for experiment CTL (constant climate control run), carried out alternately with surface accumulation data by Arthern et al. [2006] (marked AE)
and Van de Berg et al. [2006] (marked BE). The two geothermal datasets arealso applied to the paleoclimatic spin-up; Fox Maule et al. [2005] (markedFE) and Shapiro and Ritzwoller [2004] (marked SR). Note that t = 0 corresponds to the year 2004. The resolution is 20 km.

78

a

Ice surface
˚E
30

˚W

0˚

30

50

km

˚S

4.0
50

3000

˚S

3.5
2000

0˚E

6

1
2

3.0

3

1

2.5
1
2

1
3

3

2.0

90˚E

2

0

1

y (km)

1000

1.5

1

2

2

-1000

3

1

2

1.0
12
0

-2000

-2000

-1000

0

1000

2000

0.5

˚E

0.0

3000

x (km)

b

30

0˚

˚W
30

˚E

Velocity
3000

˚S

50

3000

50
˚S

1000
0˚E

2000

6

100

1000

y (km)

300

30
0

90˚E

10

m
/
a

3

-1000

1
12

-2000

0˚E

-2000

-1000

0

1000

2000

3000

0.3
0

x (km)

Figure 5.8: The result of the spin-up of the control run in 20 km resolution.
(a) Ice surface topography, (b) ice surface velocity.

790

V (m SLE)

−0.5
−1
−1.5
−2
−2.5
0

M1
M2
M3
100

200

yr

300

400

500

Figure 5.9: Simulated changes in ice volume relative to the control run CTL for experiments M1 (constant climate, sub-ice-shelf
melting 2 m ice equiv. a−1 ), M2 (constant climate, sub-ice-shelf melting 20 m ice equiv. a−1 ), M3 (constant climate, sub-ice-shelf melting
200 m ice equiv. a−1 ) The resolution is 40 km.

800V (m SLE)

−0.5

−1

M1
M2

−1.5

M3

−2
0

100

200

300

400

500

yr
Figure 5.10: Simulated changes in grounded ice volume at 40 km resolution,
V (expressed in m SLE), to the control run CTL for experiments M1 (constant climate, sub-ice-shelf melting 2 m ice equiv. a−1 ), M2 (constant climate,sub-ice-shelf melting 20 m ice equiv. a−1 ), M3 (constant climate, sub-ice-shelf
melting 200 m ice equiv. a−1 )

81

0.5

V (m SLE)

0

−0.5
CTL
2 m/a
−1

20 m/a200 m/aNo melt

−1.5
0

100

200

t (a)

300

400

500

Figure 5.11: Simulated changes in grounded ice volume, V (expressed in
m SLE), relative to the control run CTL for experiments M1 (constant climate, sub-ice-shelf melting 2 m ice equiv. a−1 ), M2 (constant climate, sub-iceshelf melting 20 m ice equiv. a−1 ), M3 (constant climate, sub-ice-shelf melting
200 m ice equiv. a−1 ) and M0 (constant climate, zero sub-ice-shelf melting).
Note that t = 0 corresponds to the year 2004. The resolution is 20 km.

82aFloating ice
1.6
1.4

A (106 km2)

1.2
1
0.8

CTL
2 m/a

0.6

20 m/a200 m/a0.4

No melt

0.2
0
0

b

100

200

t (a)

300

400

500

400

500

Grounded ice
12.8
CTL
12.7

A (106 km2)

12.6

2 m/a
20 m/a200 m/aNo melt

12.5
12.4
12.3
12.2
12.1
0

100

200

t (a)

300

Figure 5.12: Simulated grounded and floating ice area, calculated in M1-M3
experiments. t = 0 corresponds to the year 2004.

83

˚E
30

˚W

0˚

30

50

˚S

km
3
50

3000

˚S

60

2000

˚E

1

1000

y (km)

2

0

90˚E

0

-1

-1000

-2

12

-2000

0˚E

-2000

-1000

0

1000

2000

3000

-3

x (km)

Figure 5.13: The bedrock topography of the ice base. The dataset is definedin the SeaRISE project with the improved BEDMAP1 dataset [Lythe et al.,2001].

84

0.5

0

−0.5

V (m s.l.e.)

−1
C2
−1.5

C3
C4

−2

M1
M2

−2.5

M3
S1
S2

−3

S3
T1

−3.5
0

100

200

yr

300

400

500

Figure 5.14: Simulated changes in grounded ice volume, V (expressed in
m SLE), relative to the control run CTL for the experiments (1). Surface
climate changes: Blue lines, C2 (solid line, AR4 climate applied from 2004
to 2098 and constant after that) C3 (broken line, 1.5× AR4 climate, C4(chain line, 2× AR4 climate). (2). Ocean warming: Red lines, M1 (solidline, sub-ice-shelf melting, 2 m/a). M2 (broken line, sub-ice-shelf melting,
20 m/a). M3 (chain line, sub-ice-shelf melting, 200 m/a. (3). Sub-glacialchange: Black lines, S1 (solid line, double basal sliding applied), S2 (broken
line, 2.5× basal sliding applied) and S1 (solid line, 3× basal sliding applied).
(4). T1 (magenta line, AR4 climate, double basal sliding and 20 m/a subice-shelf melting is applied). The resolution is 10 km.

85aFloating ice
1.5

A (1012 km2)

C1
1

M1
M2
M3

0.5

0
0

100

b

200

yr

300

400

500

400

500

400

500

Grounded ice

A (1012 km2)

12.5

12

C1
11.5

M1
M2
M3

11
0

c

100

200

yr

300

Ice volume
16

2.65

x 10

C1
M1
2.6

M2

Volume (m3)

M3
2.55

2.5

2.45

2.4
0

100

200

yr

300

Figure 5.15: Results of enhanced sub-ice-shelf melt, M-series, experiments.
(a) Floating ice area. (b) Grounded ice area. (c) Total ice volume. t = 0
corresponds to the year 2004. The resolution is 10 km.

86aFloating ice
1.5

A (1012 km2)

1.4

1.3

C1
1.2

S1
S2
S3

1.1
0

100

b

200

yr

300

400

500

400

500

400

500

Grounded ice
12.6

A (1012 km2)

12.5
12.4
12.3
C1

12.2

S1
S2

12.1

S3
12
0

c

100

200

yr

300

Ice volume
16

2.61

x 10

2.6

C1
S1

Volume (m3)

2.59

S2
S3

2.58
2.57
2.56
2.55
2.54
0

100

200

yr

300

Figure 5.16: Results of enhanced basal sliding, S-series, experiments. (a)Floating ice area. (b) Grounded ice area. (c) Total ice volume. t = 0
corresponds to the year 2004. The resolution is 10 km.

87aFloating ice
1.6
C1

1.55

C2

km )

2

1.4

A (10

1.45

12

1.5

C3
C4

1.35
1.3
1.25
0

100

b

200

yr

300

400

500

400

500

400

500

Grounded ice
12.6
C1
C2

12.5

C4

A (10

12

2

km )

C3
12.4
12.3
12.2

0

c

100

yr

300

Ice volume

16

2.62

200

x 10

C1
C2
C3
C4

3

Volume (m )

2.6

2.58

2.56

2.54
0

100

200

yr

300

Figure 5.17: Results of surface climate change, C-series, experiments. (a)Floating ice area. (b) Grounded ice area. (c) Total ice volume. t = 0
corresponds to the year 2004. The resolution is 10 km.

88aFloating ice
2

A (1012 km2)

1.5

1
CTL
2 m/a
20 m/a0.5

200 m/aNo melt

0
0

100

b

200

t (a)

300

400

500

400

500

400

500

Grounded ice
13

A (1012 km2)

12.5

12
CTL
11.5

2 m/a
20 m/a200 m/aNo melt

0

100

c

200

t (a)

300

Ice volume1

0.5

ΔV (m SLE)

0

−0.5

−1

CTL
2 m/a
20 m/a−1.5

200 m/aNo melt

−2
0

100

200

t (a)

300

Figure 5.18: Results of sub-ice-shelf melting rates change with additional 20
year spin-up, M-series, experiments. (a) Floating ice area. (b) Grounded ice
area. (c) Sea level equivalent ice volume. t = 0 corresponds to the year 2004.
The resolution is 20 km.

89Chapter 6
Discussion
6.1

Comparison to the present-day ice sheet

The result of fixed topography spin-up is compared with a observed record
by Rignot et al. [2011a] (Fig. 6.1). The result shows that the simulated
result can represent the features of the Antarctic ice sheet flow. However,the simulated velocities are slower than the observed velocities. Considering
the fact that the present Antarctic ice sheet mass balance is negative [Rignot
et al., 2011b], ice flux is not enough to reduce ice mass inland. It would becaused by the cold ice base due to low geothermal heat flux beneath the ice
sheet as shown in Fig. 5.5. The geothermal heat flux might be warmer thanestimated one. Another possibility is that the fixed topography treatment
caused the cold ice base. Basal frictional heating depends on the basal stress.
This basal stress also depends on the surface slope. The evolution of the ice
surface may cause the additional heating.
The high sub-ice-shelf melting rate at the West Antarctic Ice sheet margin, particularly at Pine Island Bay, is considered to affect the upward glacierflow [Rignot and Jacobs, 2002]. Sub-ice-shelf melting rate influences grounding line migration directly, and the influx from grounded regions indirectly[Schoof, 2007a]. The great sensitivity of ice volume to sub-ice-shelf melting experiments is thought to be responsible for recent changes in the West
Antarctic Ice Sheet [Rignot, 2006]. However, it is inappropriate to set thesame sub-ice-shelf melting rate for the whole of the Antarctic ice shelves. The
sub-ice-shelf melting rate around the grounding line can be 10 m/a [e.g., Wenet al., 2010]. Hellmer et al. [2011] estimated the sub-ice-shelf mass balanceof ice shelves for 100 years in the future. In this estimation, although ice
shelves of the Amundsen sea region, and some part of the Ronne ice shelves,
have high sub-ice-shelf melting rates, that of the Ross ice shelf area did not90change significantly (melt rates were less than 1 m/a at the most of the shelf)
over the next 100 years.
Figure 6.2 shows ice sheet and ice shelf mass and ice thickness changes in
the C4 run. The ice shelves of the Antarctic peninsula were lost by surface
melting due to the higher surface temperature (Fig. 6.3). A part of the
Amery and the Ross ice shelf was also lost. It is suggested that the collapseof Larsen B ice shelf is caused by surface warming [van den Broeke, 2005].
These results suggest that the ice shelf is sensitive to surface warming inthese regions. Ice shelf area did not decrease in the S3 experiment (Fig. 6.4).
However, ice thickness in the East Antarctic Ice Sheet greatly decreased inthis experiment. Decrease affected not only the margins but also the inland
areas such as the upper part of the Shirase glacier. In both the experiments,
ice thickness increased around Pine Island Glacier, Marie Byrd Land. The
cause is probably low ice flux around these regions and low sub-ice-shelf
melting around the transition zone. The Pine Island Glacier region lost its
mass when ice shelf melt rates were high, such as in the M2 experiment(Fig. 6.5). Although the upper part of PIG increased in thickness evenin this experiment, it may be due to low ice flux, their area lost ice masscompared with other experiments. Most of the ice shelves disappeared in theM2 and M3 experiments (Fig. 6.6). Especially in the M3 experiment, inland
regions of the West Antarctic Ice Sheet lost mass due to extreme sub-ice-shelf
melting. Although sub-ice-shelf melting rates in this experiment were muchhigher than the observed result, it indicated the instability of the WAIS to
Antarctic Ocean warming.

6.2

The dynamics of ice sheet and ice shelves

A shallow ice approximation, shallow shelf approximation or coupling of theboth dynamics was applied to most of the ice sheet models. The shallow
shelf approximation was applied to ice shelves, the shallow ice approximation to ice sheets. The application of the shallow shelf approximation to thegrounding zone is appropriate for representing the dynamic coupling between
ice sheets and ice shelves. The coupling of the shallow ice approximation withthe shallow shelf approximation would be better for the momentum balance.
Huybrechts [1990] constructed an ice sheet model which included ice shelf
dynamics. This model used the shallow ice approximation for the ice sheet
and shallow shelf approximation for the ice shelf. A special treatment for
the grounding line zone was also applied in the ice sheet model [Huybrechts,
1990]; the coupling of the shallow ice approximation with the shallow shelf
approximation [Huybrechts, 1990]. The treatment resembles the dynamic
91

part of the Pennsylvania State University ice sheet model. The difference
was that it was applied to the whole of the ice sheet. The boundary layer
theory [Schoof, 2007b] has also been applied to the ice sheet model to estimate ice flux at the grounding line [Pollard and DeConto, 2009]. In theice sheet model GRISLI [Ritz et al., 2001], shallow shelf approximation wasapplied to the ice stream and the ice shelf and the shallow ice approximation
was applied to the ice sheet. The parallel ice sheet model (PISM) uses theshallow shelf approximation as a sliding law. The total horizontal velocity is aweighted sum of the solution of the ice shelf model with the solution of the ice
sheet model. In PISM-PIK, the horizontal velocity is just the sum of the solution of the shallow shelf approximation and the shallow ice approximation.
Weighting is not applied. Although implementation of these mechanism may
improve the performance of the model, the increasing number of shallow shelf
calculation points and the iterations are the computational approach. It it
is necessary to consider the balance between sophistication and computation
time in each case.
The other possibility is applying higher order treatments. Elmer/Ice[Zwinger et al., 2007, Seddik et al., 2012] and ISSM [Ice Sheet System Model,
Larour et al., 2012] use full stokes ice sheet dynamics. Higher order dynamics are applied to ISSM, Glimmer-CISM [Salinger et al., 2011] and IcIES[Saito et al., 2003]. The problem with this approach is that it takes much
more computation time than do applying the shallow ice approximation orshallow shelf approximation [Larour et al., 2012]. Several iteration schemes
have been proposed to reduce the computational effort required for higher
order dynamics [cf. Baral et al., 2001, Soucek and Martinec, 2008]. Ahlkrona[2011] implements second order ice sheet dynamics [Baral et al., 2001] to theice sheet model SICOPOLIS. Coupling these dynamics with the shallow shelf
approximation for ice shelves will improve the treatment of ice streams andof the transition zone.
The other difficulty for the dynamics is the elasticity of ice [Jellinek andBrill, 1956]. On short time scales, ice behaves as an elastic material. This
produces some interesting effects on glacier and ice sheet dynamics [Lingle
et al., 1981, Bindschadler et al., 2003, Gudmundsson, 2006, Sugiyama et al.,2007]. One of these effects is that the ocean tide might modulate ice shelf
flow [Brunt et al., 2010].

6.3

Grounding line migration

One of the most important controls on the ice volume change is grounding
line migration. The results of the SeaRISE experiments shows that Antarctic
92

ice volume is sensitive to the change of ice shelf or ice sheet area. This implies that grounding line migration processes are enormously important for
Antarctic ice sheet evolution. However, numerical representation of grounding line migration is difficult to incorporate in the large-scale ice sheet model.
The dynamics changes greatly in the transition zone between ice sheet andice shelf. This results in the strong dependency on resolution in simulationsof grounding line migration in the fixed grid model [Vieli and Payne, 2005].
Very high resolution, several hundred metres or a kilometre, is required torepresent grounding line migration [cf. Gladstone et al., 2012]. In in SeaRISE
experiments, however, the highest resolution is 10 km. Since this is ten times
larger than the requirement, improved treatment of the grounding line is required. Two methods have been proposed to achive this. One is the adaptive
mesh treatment (AMR) around the grounding line [Gladstone et al., 2010].
This is a means of solving the numerical equation by discretizing the domains
of interest into a grid of many individual elements. It is possible to discretize
domains around the grounding zone with fine resolution. The other is theapplication of boundary layer theory [Pollard and DeConto, 2009] aroundthe grounding zone. A two dimensional analytical solution of the flux at the
grounding line [Schoof, 2007b] is applied around the grounding zone withseveral conditions. Docquier et al. [2011] shows that application of this theory increases the performance of grounding line migration. However, thisflux at the grounding line is based merely on a two dimensional approach.
It is considered that the horizontal gradient parallel to the grounding line is
not so large in typical 2-D flow [Pollard and DeConto, 2009].
The ice sheet does not retreat from the glacial condition to the present
condition in the non-fixed topography spin up experiment. There are several
possible causes for this problem. As stated above, ice volume continuously
increases over five hundreds years in the control run. [Pollard and DeConto,
2009] assumed steady state ice volume in both glacial and interglacial conditions. The outcome was that ice volume does not change when basal melt
rate is fixed. In particular, the ice sheet does not shrink to the extreme
interglacial condition as in the other cases. This implies that sub-ice-shelf
melting is an important influence on the glacial-interglacial evolution of theice sheet. Sub-ice-shelf melt rates are constant over the glacial cycle. Thiscondition and the low ice flux of the ice sheet, prevented deglaciation during
the termination and the Holocene.
Although the grounding line migration velocity depends directly on massbalance, it also depends greatly on ice dynamics and basal flow. Certainly,
observations of the sub-ice-shelf conditions is not enough. However, there areseveral ways of estimating the sub-ice-shelf ocean conditions. The property ofbasal sliding is not known since it depends on many components (interaction93

of the ice base with the lithosphere base, dynamics of the sub-glacial till,
the basal water channels and presence of sub-glacial lakes). Observations of
the basal flow are also unknown. The model experiment of ice stream/ice
shelf/ocean coupling shows the importance of the basal flow law [Walker
et al., 2008]. Even the basal flow law itself and the methodology to estimate
these conditions have their limitations. The properties of the basal conditionis one of the largest uncertainties in estimating grounding line migration. It
would therefore be desirable to improve observations of ice base processes
and provide validation of basal sliding theories.

6.4

CalvingCalving is the most important component of mass balance for the Antarctic
ice sheet. However, the conditions for calving are still under discussion because of the complex mechanisms involved. A simple calving law is used inthis study and a threshold thickness of 50 m is applied. This value was chosen
to keep today’s ice shelf as the initial condition. However, this is considerably
thinner than the thickness of the major present-day Antarctic ice shelves. A
calving event is a consequence of fracture (crevasse) penetration. Therefore,
it is necessary consider the dynamic conditions of the ice. Moreover, there aredifference between how crevasses behave when they ultimately break. Several
types of calving law have been proposed to study the ice shelf or tidewater
glacier calving. Nye [1957] considered that crevasse depth depends on thebalance of hydrostatic pressure and the shear stresses. Such a crevasse depth
model has been applied to tidewater glaciers. Nick et al. [2010] simulated
the evolution of tidewater glaciers by using the crevasse depth model. The
sub-ice-shelf crevasses are also considered in this model. Otero et al. [2010]
included the crevasse depth model into a three dimensional glacier model.
There are other types of calving laws for glaciers. Pralong and Funk [2005]
proposed continuous damage mechanics as a theory on crevasse formation.
The rheology of damaged ice and the evolution of damage is defined in thistheory. The other difficulty of the calving study is that it does not depend
only on the viscosity of ice but also on the elastic property of ice [e.g., Plateet al., 2012]. It is thus necessary to measure stress and strain rates aroundthe calving front to validate these calving studies. Observations of seismic
waves in glaciers could also help in understanding of the calving processes.

94

Figure 6.1: Scatter map comparing simulated ice velocity (10 km resolution)with the observed dataset [Rignot et al., 2011a].

95

a3000

Mask

Y (km)

1500

0

−1500

−2800
−2800

−1500

b

0
X (km)

1500

3000

W

0

30

30
E

Ice thickness difference
m
1000

3000

500
E
60

1500

200
100

y (km)

50

90E

0

0
-50
-100

-1500

-200
12
0E

-1500

0

1500

3000

-500
-1000

x (km)

Figure 6.2: (a) Simulated ice sheet/ice shelf mask of the last time step and(b) ice thickness change from the initial condition to the last time step,
calculated in C4 experiments. t = 0 corresponds to the year 2004.

96

E

30

30

0

W

C
-10

3000

-20
0E

6

-25

1500

y (km)

-30

90E

0

-35
-40
-45

-1500
12

0E

-1500

0

1500

3000

-50
-55

x (km)

Figure 6.3: Annual surface temperature on the ice sheet at the last time stepof C4 experiment.

97

a

3000

Mask

Y (km)

1500

0

−1500

−2800
−2800

−1500

b

0
X (km)

1500

3000

W

0

30

30
E

Ice thickness difference
m
1000

3000

500
E
60

1500

200
100

y (km)

50

90E

0

0
-50
-100

-1500

-200
12
0E

-1500

0

1500

3000

-500
-1000

x (km)

Figure 6.4: (a) Simulated ice sheet/ice shelf mask of the last time step and(b) ice thickness change from the initial condition to the last time step,
calculated in S3 experiments. t = 0 corresponds to the year 2004.

98

3000

a

Mask

Y (km)

1500

0

−1500

−2800
−2800

−1500

b

0
X (km)

1500

3000

W

0

30

30
E

Ice thickness difference
3000

E
60

y (km)

1500

90E

0

-1500
12
0E

-1500

0

1500

3000

x (km)

Figure 6.5: (a) Simulated ice sheet/ice shelf mask of the last time step and(b) ice thickness change from the initial condition to the last time step,
calculated in M3 experiments. t = 0 corresponds to the year 2004.

99

a

3000

Mask

Y (km)

1500

0

−1500

−2800
−2800

−1500

b

0
X (km)

1500

3000

W

0

30

30
E

Ice thickness difference
m
1000

3000

500
E
60

1500

200
100

y (km)

50

90E

0

0
-50
-100

-1500

-200
12
0E

-1500

0

1500

3000

-500
-1000

x (km)

Figure 6.6: (a) Simulated ice sheet/ice shelf mask of the last time step and(b) ice thickness change from the initial condition to the last time step,
calculated in M3 experiments. t = 0 corresponds to the year 2004.

100

Chapter 7
Summary
A theory of ice sheet dynamics on orthogonal curvilinear coordinates wasconsidered and applied to shallow shelf approximation. The application of
the theory and its coupling to the ice sheet was the next step. A module of ice shelf dynamics was therefore constructed and coupled to the dynamic/thermodynamic ice sheet model SICOPOLIS. The performance of theice shelf model was validated against the analytical solution of the shallow
shelf approximation, the ice shelf ramp. The results show that the model
well represents the analytical solution.
The new version of the model is applied to the Antarctic ice sheet by using the SeaRISE dataset. A paleo-climatic spin-up over 250 kyr was applied
to find today’s condition of the ice sheet. The non-fixed topography spin-up
produced the ice sheet conditions during the glacial period. However, there isa problem involved with deglaciation process that occurred between the lastglacial maximum and the present condition. The reasons underlying thisproblem are low grounded ice flux and constant, low sub-ice shelf melting
rates. A fixed topography spin-up was finally applied to future climate scenarios. The sensitivities were checked of different datasets of geothermal heat
flux [Shapiro and Ritzwoller, 2004, Fox Maule et al., 2005] and accumulation
rates [Arthern et al., 2006, Van de Berg et al., 2006]. The results showedthat sea level equivalent ice volume are not very sensitive to geothermal heat
flux difference but to accumulation rate differences. Sea level equivalent icevolume based on the Arthern et al. [2006] dataset is 0.4 m SLE lower thanthat with the Van de Berg et al. [2006] dataset. The difference is causedby the high accumulation rates around the ice sheet margin in Van de Berget al. [2006] dataset.
Sensitivities to sub-ice-shelf melting rates was evaluated with SeaRISE Mseries experiments, where surface climate forcing was kept steady at presentday conditions while the sub-ice-shelf melting rates were varied over 0 −101

200 m/a. The range of grounded ice volume changed across the M-series
being ∼ 0.5 − 0.7 m SLE after 100 years and ∼ 1.5 − 2.0 m SLE after 500
years in 40 km runs. In 20 km runs, it is ∼ 1.3 m SLE after 100 years and∼ 2.25 m SLE after 500 years. The differences between the results of eachexperiment were caused not only by the settings but also by the resolution.
This was because resolution affects the basal topography of the bedrock andthe grounding line migration dynamics themselves. All of these simulations
showed that ice sheet volume is greatly sensitive to sub-ice-shelf melting rate.
Ice sheet volume correlated strongly with floating area in these experiments.
These results indicate the factors likely to be responsible for recent changes
in the West Antarctic Ice Sheet [e.g., Rignot, 2006] and agree with findings
by Huybrechts [1990].
The range of outcomes of all of SeaRISE experiments in 10 km resolutionwas ∼ 1.3 m SLE after 100 years and ∼ 3.3 m SLE after 500 years. Although
ice volume is sensitive to sub-ice-shelf melting, basal sliding is also important. The changes caused by surface climate experiments were small. The
combination run, T1, is the worst case (the coupling of C2, S1 with M2).
It implies that the combination of several forcing factors would cause additional ice loss. These experiments suggest the Antarctic ice sheet is sensitive
to climate change. The sea level change caused by the Antarctic ice sheet
may not be negligible in the future. Further studies of sub-ice-shelf melting
and basal sliding processes are needed because outcomes are quite sensitive
to these processes.

102

Bibliography
J. Ahlkrona. Implementing higher order dynamics into the ice sheet model
SICOPOLIS. Master’s thesis, Department of Information Technology,
Uppsala University, Sweden, 2011.
A. Arakawa and V. R. Lamb. Computational design of the basic dynamical processes of the UCLA general circulation model. In J. Chang, editor, Methods in Computational Physics Vol. 17, pages 173–265. Academic
Press, New York, NY, USA, 1977.
R. J. Arthern, D. P. Winebrenner, and D. G. Vaughan. Antarctic snow
accumulation mapped using polarization of 4.3-cm wavelength microwave
emission. Journal of Geophysical Research, 111(D6):D06107, 2006.
L. Augustin, C. Barbante, P. R. F. Barnes, J. M. Barnola, M. Bigler,
E. Castellano, O. Cattani, J. Chappellaz, D. Dahl-Jensen, B. Delmonte,
et al. Eight glacial cycles from an Antarctic ice core. Nature, 429(6992):
623–628, 2004.
D. R. Baral, K. Hutter, and R. Greve. Asymptotic theories of large-scale
motion, temperature, and moisture distribution in land-based polythermal
ice sheets: a critical review and new developments. Applied Mechanics
Reviews, 54:215, 2001.
R. A. Bindschadler, M. A. King, R. B. Alley, S. Anandakrishnan, and L. Padman. Tidally controlled stick-slip discharge of a West Antarctic ice. Science, 301(5636):1087–1089, 2003.
K. M. Brunt, M. A. King, H. A. Fricker, and D. R. MacAyeal. Flow of theRoss Ice Shelf, Antarctica, is modulated by the ocean tide. Journal of
Glaciology, 56(195):157–161, 2010.
E. Bueler and J. Brown. Shallow shelf approximation as a “sliding law” ina thermomechanically coupled ice sheet model. Journal of GeophysicalResearch, 114(F03008):F03008, 2009.
103

R. Calov. Das thermomechanische Verhalten des Grönländischen Eisschildes unter der Wirkung verschiedener Klimaszenarien – Antworten eines
theoretisch-numerischen Modells. Doctoral thesis, Department of Mechanics, Darmstadt University of Technology, Germany, 1994.
R. Calov, A. A. Savvin, R. Greve, I. Hansen, and K. Hutter. Simulationof the Antarctic ice sheet with a three-dimensional polythermal ice-sheet
model, in support of the EPICA project. Annals of Glaciology, 27:201–206,
1998.
J. C. Comiso. Variability and trends in Antarctic surface temperatures fromin situ and satellite infrared measurements. Journal of Climate, 13(10):1674–1696, 2000.
K. M. Cuffey and W. S. B. Paterson. The Physics of Glaciers. Elsevier, Amsterdam, The Netherlands etc., 4th edition, 2010. ISBN 978-0-12-369461-4.
H. De Angelis and P. Skvarca. Glacier surge after ice shelf collapse. Science,
299(5612):1560, 2003.
C. M. S. Doake. Ice-shelf stability. In D. G. Ainley, G. J. Divoky, J. Steele,
S. Thorpe, and K. Turekian, editors, Encyclopedia of Ocean Sciences, pages
1282–1290. Academic Press, San Diego, CA, USA, 2001.
D. Docquier, L. Perichon, and F. Pattyn. Representing grounding line dynamics in numerical ice sheet models: Recent advances and outlook. Surveys in Geophysics, 32(4-5):417–435, 2011.
P. L. Forsström and R. Greve. Simulation of the Eurasian ice sheet dynamics
during the last glaciation. Global and Planetary Change, 42(1-4):59–81,
2004.
J. P. F. Fortuin and J. Oerlemans. Parameterization of the annual surface
temperature and mass balance of Antarctica. Annals of Glaciology, 14:
78–84, 1990.
C. Fox Maule, M. E. Purucker, N. Olsen, and K. Mosegaard. Heat flux
anomalies in Antarctica revealed by satellite magnetic data. Science, 309(5733):464–467, 2005.
R. M. Gladstone, V. Lee, A. Vieli, and A. J. Payne. Grounding line migrationin an adaptive mesh ice sheet model. Journal of Geophysical Research, 115(F4):F04014, 2010.

104

R. M. Gladstone, A. J. Payne, and S. L. Cornford. Resolution requirements
for grounding-line modelling: sensitivity to basal drag and ice-shelf buttressing. Annals of Glaciology, 53(60):97, 2012.
K. Goto-Azuma. Millennial-scale climate variability during the past 720,000
years recorded in the Dome Fuji ice core. In Geophysical Research Abstracts, volume 10, 2008.
R. Greve. Thermomechanisches Verhalten polythermer Eisschilde – Theorie,
Analytik, Numerik. Doctoral thesis, Department of Mechanics, Darmstadt University of Technology, Germany, 1995. Berichte aus der Geowissenschaft, Shaker Verlag, Aachen, Germany.
R. Greve. Application of a polythermal three-dimensional ice sheet modelto the Greenland ice sheet: response to steady-state and transient climate
scenarios. Journal of Climate, 10(5):901–918, 1997a.
R. Greve. A continuum-mechanical formulation for shallow polythermal ice
sheets. Philosophical Transactions of the Royal Society of London. Series
A: Mathematical, Physical and Engineering Sciences, 355(1726):921–974,
1997b.
R. Greve and H. Blatter. Dynamics of Ice Sheets and Glaciers. Springer,
Berlin, Germany etc., 2009. ISBN 978-3-642-03414-5. doi: 10.1007/978-3642-03415-2.
G. H. Gudmundsson. Fortnightly variations in the flow velocity of Rutford
Ice Stream, West Antarctica. Nature, 444(7122):1063–1064, 2006.
H. Gudmundsson. Continuum mechanics. Karthaus Summer School Lecture
Notes, 2007.
H. Heinrich. Origin and consequences of cyclic ice rafting in the northeast
atlantic ocean during the past 130,000 years. Quaternary Research, 29(2):142–152, 1988.
H. Hellmer, R. Timmermann, and F. Kauker. Will another ice shelf soon
float in ‘hot’ water? In International Symposium on “Interactions of Ice
Sheets and Glaciers with the Ocean”, 610 June 2011, La Jolla, CA, USA,2011.
T. Hondoh, editor. Nature and behavior of dislocations in ice, Sapporo,
Japan, 2000. Hokkaido University Press.

105

P. Huybrechts. A 3-D model for the Antarctic ice sheet: a sensitivity study
on the glacial-interglacial contrast. Climate Dynamics, 5:79–92, 1990.
P. Huybrechts and J. de Wolde. The dynamic response of the Greenlandand Antarctic ice sheets to multiple-century climatic warming. Journal of
Climate, 12(8):2169–2188, 1999.
H. H. G. Jellinek and R. Brill. Viscoelastic properties of ice. Journal ofApplied Physics, 27(10):1198–1209, 1956.
A. Jenkins, D. G. Vaughan, S. S. Jacobs, H. H. Hellmer, and J. R. Keys.
Glaciological and oceanographic evidence of high melt rates beneath Pine
Island Glacier, West Antarctica. Journal of Glaciology, 43(143):114–121,
1997.
D. Jenssen. A three-dimensional polar ice-sheet model. Journal of Glaciology,
18:373–389, 1977.
I. Joughin and L. Padman. Melting and freezing beneath Filchner-Ronne ice
shelf, Antarctica. Geophysical Research Letters, 30(9):1477, 2003.
I. Joughin, E. Rignot, C. E. Rosanova, B. K. Lucchitta, and J. Bohlander.
Timing of recent accelerations of Pine Island Glacier, Antarctica. Geophysical Research Letters, 30(13):1706, 2003.
L. D. Landau and E. M. Lifshitz. Course of Theoretical Physics – The classical theory of fields. Butterworth-Heinemann, 1980.
L. D. Landau and E. M. Lifshitz. Course of Theoretical Physics – Fluid
Mechanics. Butterworth-Heinemann, 1987.
E. Larour, H. Seroussi, M. Morlighem, and E. Rignot. Continental
scale, high order, high spatial resolution, ice sheet modeling using theice sheet system model (issm). Journal of Geophysical Research, 117:
doi:10.1029/2011JF002140, 2012.
C. S. Lingle, T. J. Hughes, and R. C. Kollmeyer. Tidal flexure of Jakobshavns
Glacier, West Greenland. Journal of Geophysical Research, 86:3960, 1981.
M. Llubes, C. Lanseau, and F. Rémy. Relations between basal condition,
subglacial hydrological networks and geothermal flux in Antarctica. Earthand Planetary Science Letters, 241(3-4):655–662, 2006.

106

M. B. Lythe, D. G. Vaughan, et al. BEDMAP: a new ice thickness and subglacial topographic model of Antarctica. Journal of Geophysical Research,
106(B6):11335–11, 2001.
V. I. Morgan, T. H. Jacka, G. J. Akerman, and A. L. Clarke. Outlet glacierand mass-budget studies in Enderby, Kemp and MacRobertson lands,
Antarctica. Annals of Glaciology, 3:204–210, 1982.
F. M. Nick, C. J. Van der Veen, A. Vieli, and D. I. Benn. A physically based
calving model applied to marine outlet glaciers and implications for theglacier dynamics. Journal of Glaciology, 56(199):781–794, 2010.
J. F. Nye. The distribution of stress and velocity in glaciers and ice-sheets.
Proceedings of the Royal Society of London. Series A. Mathematical andPhysical Sciences, 239(1216):113–133, 1957.
J. Otero, F. J. Navarro, C. Martin, M. L. Cuadrado, and M. I. Corcuera.
A three-dimensional calving model: numerical experiments on Johnsons
Glacier, Livingston Island, Antarctica. Journal of Glaciology, 56(196):
200–214, 2010.
F. Pattyn, A. Huyghe, S. De Brabander, and B. De Smedt. Role of transition
zones in marine ice sheet dynamics. Journal of Geophysical Research, 111(F02004):F02004, 2006.
A. J. Payne. A thermomechanical model of ice flow in West Antarctica.
Climate Dynamics, 15(2):115–125, 1999.
J. R. Petit, J. Jouzel, D. Raynaud, N. I. Barkov, J. M. Barnola, I. Basile,
M. Bender, J. Chappellaz, M. Davis, G. Delaygue, M. Delmotte, V. M.
Kotlyakov, M. Legrand, V. Y. Lipenkov, C. Lorius, L. Pepin, C. Ritz,
E. Saltzman, and M. Stievenard. Climate and atmospheric history of thepast 420,000 years from the Vostok ice core, Antarctica. Nature, 399(6735):
429–436, 1999.
C. Plate, R. Müller, A. Humbert, and D. Gross. Evaluation of the criticalityof cracks in ice shelves using finite element simulations. The Cryosphere
Discuss., 6:469–503, 2012.
D. Pollard and R. M. DeConto. Modelling West Antarctic ice sheet growthand collapse through the past five million years. Nature, 458(7236):329–
332, 2009.

107

A. Pralong and M. Funk. Dynamic damage model of crevasse opening andapplication to glacier calving. Journal of Geophysical Research, 110(B1):
B01309, 2005.
E. Rignot. Changes in ice dynamics and mass balance of the Antarctic ice
sheet. Philosophical Transactions of the Royal Society A: Mathematical,
Physical and Engineering Sciences, 364(1844):1637–1655, 2006.
E. Rignot and S. S. Jacobs. Rapid bottom melting widespread near Antarctic
ice sheet grounding lines. Science, 296(5575):2020, 2002.
E. Rignot, D. G. Vaughan, M. Schmeltz, T. Dupont, and D. MacAyeal. Acceleration of Pine Island and Thwaites Glaciers, West Antarctica. Annals
of Glaciology, 34(1):189–194, 2002.
E. Rignot, J. Mouginot, and B. Scheuchl. Ice flow of the Antarctic ice sheet.
Science, 333(6048):1427–1430, 2011a.
E. Rignot, I. Velicogna, M. R. Van den Broeke, A. Monaghan, andJ. Lenaerts. Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise. Geophysical Research Letters, 38(5):L05503,
2011b.
C. Ritz, V. Rommelaere, and C. Dumas. Modeling the evolution of Antarctic
ice sheet over the last 420,000 years: Implications for altitude changes
in the Vostok region. Journal of Geophysical Research, 106(D23):31943–
31964, 2001.
H. Rott, W. Rack, P. Skvarca, and H. De Angelis. Northern Larsen ice shelf,
Antarctica: further retreat after collapse. Annals of Glaciology, 34(1):277–282, 2002.
F. Saito, A. Abe-Ouchi, and H. Blatter. Effects of first-order stress gradientsin an ice sheet evaluated by a three-dimensional thermomechanical coupledmodel. Annals of Glaciology, 37(1):166–172, 2003.
A. Salinger, K. J. Evans, J. F. Lemieux, D. Holland, T. Payne, S. Price, andD. Knoll. Implementation of the jacobian-free newton-krylov method forsolving the for solving the first-order ice sheet momentum balance. Journalof Computational Physics, 230(17), 2011.
T. A. Scambos, J. A. Bohlander, C. A. Shuman, P. Skvarca, et al. Glacier acceleration and thinning after ice shelf collapse in the Larsen B embayment,
Antarctica. Geophysical Research Letters, 31, 2004.
108

R. P. Scherer, A. Aldahan, S. Tulaczyk, G. Possnert, H. Engelhardt, and
B. Kamb. Pleistocene collapse of the West Antarctic ice sheet. Science,
281(5373):82, 1998.
C. Schoof. Ice sheet grounding line dynamics: Steady states, stability, andhysteresis. Journal of Geophysical Research, 112(F3):F03S28, 2007a.
C. Schoof. Marine ice-sheet dynamics. Part 1. the case of rapid sliding.
Journal of Fluid Mechanics, 573(1):27–55, 2007b.
H. Seddik, R. Greve, T. Zwinger, and L. Placidi. A full Stokes ice flow modelfor the vicinity of Dome Fuji, Antarctica, with induced anisotropy andfabric evolution. The Cryosphere, 5(2):495–508, 2011. doi: 10.5194/tc-5495-2011.
H. Seddik, R. Greve, T. Zwinger, F. Gillet-Chaulet, and O. Gagliardini.
Simulations of the Greenland ice sheet 100 years into the future with thefull Stokes model Elmer/Ice. Journal of Glaciology, 58(209):427, 2012.
N. M. Shapiro and M. H. Ritzwoller. Inferring surface heat flux distributions
guided by a global seismic model: particular application to Antarctica.
Earth and Planetary Science Letters, 223(1-2):213–224, 2004.
S. Solomon, D. Qin, M. Manning, Z. Chen, M. Marquis, K. B. Averyt,
M. Tignor, and H. L. Miller, editors. Contribution of Working Group I to
the Fourth Assessment Report of the Intergovernmental Panel on Climate
Change, 2007. Cambridge University Press, Cambridge, United Kingdom
and New York, NY, USA, 2007.
O. Soucek and Z. Martinec. Iterative improvement of the shallow-ice approximation. Journal of Glaciology, 54(188):812–822, 2008.
S. Sugiyama, A. Bauder, P. Weiss, and M. Funk. Reversal of ice motion during the outburst of a glacier-dammed lake on gornergletscher, switzerland.
Journal of Glaciology, 53(181):172–180, 2007.
W. J. Van de Berg, M. R. Van den Broeke, C. H. Reijmer, and E. van Meijgaard. Reassessment of the Antarctic surface mass balance using calibrated
output of a regional atmospheric climate model. Journal of GeophysicalResearch, 111(D11):D11104, 2006.
M. van den Broeke. Strong surface melting preceded collapse of Antarctic
Peninsula ice shelf. Geophysical Research Letters, 32(12):L12815, 2005.

109

A. Vieli and A. J. Payne. Assessing the ability of numerical ice sheet models
to simulate grounding line migration. Journal of Geophysical Research,
110:F01003, 2005.
R. T. Walker, T. K. Dupont, B. R. Parizek, and R. B. Alley. Effects of basalmelting distribution on the retreat of ice-shelf grounding lines. GeophysicalResearch Letters, 35:L17503, 2008.
J. Weertman. Stability of the junction of an ice sheet and an ice shelf. Journalof Glaciology, 13:3–13, 1994.
J. Wen, Y. Wang, W. Wang, K. C. Jezek, H. Liu, and I. Allison. Basalmelting and freezing under the Amery Ice Shelf, East Antarctica. Journalof Glaciology, 56(195):81–90, 2010.
R. S. Williams, J. G. Ferrigno, and K. M. Foley. Coastal-change and glaciological maps of Antarctica. U.S. Geological Survey Fact Sheet 2005-3055,
2 p., 2005. URL http://pubs.usgs.gov/fs/2005/3055.
T. Zwinger, R. Greve, O. Gagliardini, T. Shiraiwa, and M. Lyly. A full Stokesflow thermo-mechanical model for firn and ice applied to the Gorshkov
crater glacier, Kamchatka. Annals of Glaciology, 45(1):29–37, 2007.

110Reference Papers

北海道の雪氷

No.30（2011）

棚氷の縮小による南極氷床の流動変化

佐藤建（北海道大学大学院環境科学院），Ralf Greve（北海道大学低温科学研究所）

１．はじめに近年の大気、海洋の温暖化による南極の棚氷が崩壊や, 縮小が観測されている. このために棚氷が氷床を抑える効果が弱まり, 上流部の氷の流出が促進されることが観測, 理論の両面から示唆されている 1 ) .
棚氷と氷床の内陸部では流動メカニズムが大きく異なる. 氷床内陸部では, 底面の基盤のとの摩擦のため, 鉛直シア応力が最も主要な応力である. 一方, 棚氷の底面は海洋であり底面摩擦が働かない. このため鉛直シア応力は無視できるほどで, 水平方向の応力が主たる応力の成分となる. 下部の底面摩擦がないために典型的な棚氷の流動速度は 1000 m/a であり, 典型的な変動時間スケールは数百キロの大きな棚氷でも氷
床内陸の 1/10 以下である. したがって棚氷は外部の環境変化に対して内陸部より速く反応する.
南極氷床の変動による将来の海水準変化への影響には大きな不確実性があり, その原因の一つが棚氷の変動によるものである. 棚氷底面の融解量を変動させることで、棚氷の変動による氷床変動への影響を考察する.
２．手法
本研究では氷床モデル SICOPOLIS に棚氷の流動モデルを結合させて氷床全体のシミュレーションを行った. 熱動力学結合氷床モデル SICOPOLIS 2 ) はこれまでグリーンラン
ド、南極、ローレンタイド氷床、火星の氷床等の様々な地域に適用されているが, 棚
氷への扱いがなされていなかった. 筆者は棚氷の流動モデルを作成し、氷床モデルと結合してこれを南極氷床へと適用した.
筆者らは氷床による海面上昇を予測する国際プロジェクト SeaRISE に参画している.
ここでは与えられた氷厚, 地熱, 涵養量, 表面温度, 棚氷-氷床分布を使用してシミュ
レーションを行った. モデルの水平分解能は 40km, 鉛直方向は 90 層である.
棚氷底面の融解量については, 海域により異なる扱いを行った. コントロール実験
ではロス棚氷, ロンヌ棚氷, アメリー棚氷の融解量については, 氷床との境界である
接地線近辺では 2 m/a, それ以外の部分では 0.2m/a, 深度 2000m 以下の外洋では 20m/aとした. 陸氷に囲まれていない棚氷や南極半島の棚氷の領域では高い融解量が考えられるため, その 5 倍の融解量を設定した . 次に海洋での融解量を 2 倍、5 倍、10 倍として融解量変動への鋭敏性を調べた.
３．結果と考察
氷床の流動速度は図-1 のようになっており、棚氷や沿岸部で高い速度をもつこと,
棚氷では沿岸部に向かってより速い速度を持つことが再現されている. 次に棚氷の融解量変動実験を行った. 棚氷の融解による内陸氷床の変動を見積もるために海水準相
当氷床変動を求めた. これにより, 棚氷の融解による体積変化を除き, かつ将来の海- 73 Copyright ○c 2011(社) 日本雪氷学会北海道支部

北海道の雪氷

No.30（2011）

水準変化の影響を求める事ができる.
実験の結果、融解量が大きいほど海水準相当の氷床体積が減少する効果を見積もることができた（図-2）. 融
解量を 2 倍とした実験と参照実験では 500 年間で海水準変動に 70cm 程度の差を引き起こす. 参照実験と融解なしの場合でも同様である. また融
解量 5 倍と 10 倍の差は 500 年間でも10cm 程度である. これは融解量が 5倍の状態でも十分に融解量が大きく,
多くの棚氷が無くなり、融解量を上げても変動がないためだろう.
海洋温暖化の実験 4 ) では 21 世紀後半で

図-1.

南極氷床の速度分布の初期値

ロンヌ棚氷下部における融解量が 18m/a を超える事が示唆されており, 温暖化に鋭敏な地域では数値実験で使用したような高い融解量を持つ可能性がある. 一方ロス棚氷の融解量変動は大きくなく, 海域による違いが重要になる.
４．結論
海域による棚氷底部の融解量の差を考慮した氷床変動の数値実験を行った。
融解量設定の変動により、海水準相当の氷
床体積は 100 年で 70cm、500 年で 2m となった. 底面での融解量の違いに対して,

図-2. 融解量 2 倍, 5 倍, 10 倍, 融解なし

氷床鋭敏に変動する. 今後は棚氷-氷床のの各実験の氷床体積変化の参照実験か

面積変化を考える必要がある.

らの差【参考・引用文献】
1) Rignot, E., (2006): Changes in ice dynamics and mass balance of the
Antarctic ice sheet, Philosophical Transactions A,364, 1844 1637-1655,
2) Greve, R., (1997): Application of a polythermal three-dimensional ice sheet
model to the Greenland ice sheet: response to steadystate and transient
climate scenarios. Journal of Climate, 10(5), 901-918,
3) Pollard, D. (2009): Response of the Antarctic Ice Sheet to increased ice-shelfoceanic melting, WAIS/FRISP Workshop, Washington, USA.
4) Hellmer,H.(2011): Will another ice shelf soon float in ‘hot’ water?
International Glaciological Society Symposium on Interactions of Ice
Sheets and Glaciers with the Ocean, La Jolla, USA.

- 74 Copyright ○c 2011(社) 日本雪氷学会北海道支部

Annals of Glaciology 00 0000

1

1

Sensitivity experiments for the Antarctic ice sheet with

2varied sub-ice-shelf melting rates

3

Tatsuru SATO,1,2 Ralf GREVE2
1

4

Graduate School of Environmental Science, Hokkaido University,
Kita-10, Nishi-5, Kita-ku, Sapporo 060-0810, Japan

5

2

6

Institute of Low Temperature Science, Hokkaido University,

7

Kita-19, Nishi-8, Kita-ku, Sapporo 060-0819, Japan

8

ABSTRACT. Ice sheet modelling is an important tool for predicting the possible response9

of ice sheets to climate change in the past and future. An established ice sheet model is

10

SICOPOLIS (SImulation COde for POLythermal Ice Sheets), and for this study the previ-

11

ously grounded-ice-only model was complemented by an ice shelf module. The new version

12of SICOPOLIS is applied to the Antarctic ice sheet, driven by standard forcings defined by

13

the SeaRISE (Sea-level Response to Ice Sheet Evolution) community effort. A crucial point

14

for simulations into the future is to obtain reasonable initial conditions by a paleoclimatic15

spin-up, which we carry out over 125,000 years from the Eemian until today. We then carry16

out a set of experiments for 500 years into the future, in which the surface temperature and17

precipitation are kept at their present-day distributions, while sub-ice-shelf melting rates be-

18

tween 0 and 200 m a−1 are applied. These simulations show a significant, but not catastrophic

19

sensitivity of the ice sheet. Grounded ice volumes decrease with increasing melting rates,

20

and the spread of the results from the zero to the maximum melting case is ∼ 0.65 metres of21

sea level equivalent after 100 years and ∼ 2.25 metres of sea level equivalent after 500 years.

22

1. INTRODUCTION

23

The Antarctic ice sheet is by far the largest land ice mass on the present-day Earth, and its volume amounts to approximatively

24

56.6 m SLE (metres of sea level equivalent) (Lemke and others, 2007). The current mass balance of the ice sheet is most likely

2Sato and Greve: Sensitivity experiments for Antarctica

25

negative with an accelerating trend, even though the uncertainty is still significant (Lemke and others, 2007; Rignot and others,

26

2011; van den Broeke and others, 2011; Zwally and Giovinetto, 2011). Owing to the very low surface temperatures, the ice

27

sheet is not very susceptible to surface melting even under moderate global warming scenarios. However, recent observations

28(e.g., Shepherd and others, 2001; Joughin and others, 2003; Rignot and others, 2004) have led to strong concerns that ice-29

dynamical processes (loss of buttressing from ice shelves, speed-up of ice streams and outlet glaciers) may boost the mass

30

loss and thus lead to an additional contribution to sea level rise. The limited understanding of such processes was explicitly

31

highlighted in the Fourth Assessment Report (AR4) of the United Nations Intergovernmental Panel on Climate Change

32(IPCC): “Dynamical processes related to ice flow not included in current models but suggested by recent observations could

33

increase the vulnerability of the ice sheets to warming, increasing future sea level rise. Understanding of these processes is34

limited and there is no consensus on their magnitude.” (IPCC, 2007).

35

The scientific community has reacted to the need for an improved modelling of ice sheet dynamics, including their ice

36

shelves, ice streams and outlet glaciers. Coordinated research projects have been launched, such as the European-led ice2sea

37

programme funded by the European Union Framework-7 scheme (http://www.ice2sea.eu/), or the US-led SeaRISE effort

38(Sea-level Response to Ice Sheet Evolution; http://tinyurl.com/srise-lanl, http://tinyurl.com/srise-umt). The latter39

is an already well-advanced, community-organized effort to estimate the likely range of ice sheet contributions to sea level rise

40

over the next few hundred years.

41

The Japanese ice sheet modelling community contributes to both ice2sea and SeaRISE as part of several funded research

42

projects. In this study, we present results obtained with the model SICOPOLIS (SImulation COde for POLythermal Ice Sheets)

43

for the Antarctic ice sheet within the framework of the SeaRISE effort. To this end, SICOPOLIS was complemented by an44

ice-shelf module. We present a paleoclimatic spin-up that provides reasonable initial conditions for future climate runs. Based45

on that, we conduct a series of future climate scenarios specified by SeaRISE in which the climatic conditions at the ice surface46

are kept at their present-day conditions, but increased sub-ice-shelf melting as a consequence of increased ocean temperatures

47

is assumed. The impact of these forcings on the mass balance of the ice sheet over the next 500 years is discussed.

48

2. ICE SHEET MODEL SICOPOLIS

49

SICOPOLIS simulates the large-scale dynamics and thermodynamics (ice extent, thickness, velocity, temperature, water con-50

tent and age) of ice sheets three-dimensionally and as a function of time (Greve, 1997; for the latest version 3.0 used here51

see http://sicopolis.greveweb.net/). It is based on the shallow ice approximation (SIA) for grounded ice (Hutter, 1983;

52

Morland, 1984), the shallow shelf approximation (SSA) for floating ice (Morland, 1987; MacAyeal, 1989; Weis and others,

531999) and the rheology of an incompressible, heat-conducting, power law fluid.

Sato and Greve: Sensitivity experiments for Antarctica

54

55

3

In this study, we use the regularized Glen flow law for the ice fluidity (inverse viscosity) 1/η in the form of Greve and Blatter(2009),
1= 2EA(T 0 )f (σe ) ,η(T 0 , σe )

(1)

q

56

where T 0 is the temperature relative to pressure melting1 , A(T 0 ) the rate factor, σe =57

root of the second invariant of the deviatoric stress tensor tD ), f (σe ) the creep function and E the flow enhancement factor.

58

The rate factor and creep function are expressed in the form of an Arrhenius law,A(T 0 ) = A0 e−Q/RT

1
D 2
2 tr (t ) the effective stress (square

0

(2)

59

(A0 : pre-exponential constant, Q: activation energy, R: universal gas constant, T 0 in K), and a power law with additional

60

constant term,f (σe ) = σen−1 + σ0n−1

(3)

61

(n: stress exponent, σ0 : residual stress), respectively. The values for the several parameters follow largely the recent recom-

62

mendations by Cuffey and Paterson (2010) and are listed, among others, in Table 1.

63

A particular feature of the model thermodynamics is that it distinguishes between cold ice with a temperature below the64

pressure melting point and temperate ice with a temperature at the pressure melting point, the latter being considered as a65binary mixture of ice and small amounts of water. In temperate ice, following Lliboutry and Duval (1985), the temperature-66

dependent rate factor A(T 0 ) is replaced by a rate factor depending on the water content W ,At (W ) = A(T 0 = 0◦ C) × (1 + 1.8125 W [%]) .

(4)

67

The interface that separates cold and temperate ice is tracked through the use of Stefan-type energy flux and mass flux

68

matching conditions (this procedure is referred to as the “polythermal mode”).

69

Basal sliding under grounded ice, vb , is described by a Weertman-type sliding law with sub-melt sliding of the form applied70

to the Austfonna ice cap by Dunse and others (2011) (based on earlier work by, e.g., Weertman (1964), Lliboutry (1968),

71

Bindschadler (1983), Hindmarsh and Le Meur (2001)),0vb (Tb0 ) = −Cb eTb /γ

τbp
,
Nbq

(5)

72

where τb is the basal drag (shear stress), Nb the effective pressure and Tb0 the basal temperature relative to pressure melting73(in ◦ C, always ≤ 0◦ C). In the SIA, the effective pressure Nb is equal to the hydrostatic pressure Pb = ρgH (where ρ is the
1 More precisely, following Greve and Blatter (2009), T 0 = T − T + T , where T is the ice temperature (in either K or ◦ C), T = T − βpmm00

the pressure (p) dependent melting point and T0 the melting temperature at standard atmospheric pressure (Table 1).

4Sato and Greve: Sensitivity experiments for Antarctica

74

density of ice, g the gravitational acceleration and H the ice thickness) minus the water pressure Pw . The latter is assumed

75

to be equal to the pressure of a sea water column of density ρsw and thickness Hw = zsl − b, where zsl is the mean sea level76

and b the ice base topography. This yields

Nb = Pb − Pw = ρgH − ρsw gHw .

(6)

77

In order to avoid the Nb → 0 singularity when grounded ice becomes nearly floating close to the grounding line, Nb is78

constrained to be ≥ 0.2 × Pb . The basal drag τb is equal to the basal pressure times the surface slope,τb = ρgH |grad h| ,

(7)

79

where h is the surface elevation and grad the gradient operator in the horizontal plane (Greve and Blatter, 2009). The minus

80

sign in Eq. (5) indicates that the direction of basal sliding is anti-parallel to the basal drag. The parameters have the values

81

Cb = 11.2 m a−1 Pa−1 (sliding coefficient), p = 3, q = 2 (sliding exponents) and γ = 1◦ C (sub-melt-sliding parameter).

82

83

Isostatic depression and rebound of the lithosphere due to changing ice load is modelled by the elastic lithosphere/relaxing
asthenosphere (ELRA) approach (Le Meur and Huybrechts, 1996; Greve, 2001).

84

The shallow shelf module is a new feature of version 3.0 of SICOPOLIS. Its core is the solver for the non-linear elliptical85

differential equations for the horizontal velocities vx and vy (Greve and Blatter, 2009, Eq. (6.55); note that in the SSA vx

86

and vy are constant over depth). This is done by discretizing the equations by central differences on the regular grid of

87

SICOPOLIS, solving the resulting system of linear equations by the biconjugate gradient method with the Library of Iterative

88

Solvers for Linear Systems (Lis, http://www.ssisc.org/lis/), and non-linear iterations between the horizontal velocities and89

the depth-integrated viscosity with a relaxation scheme. The position of the grounding line is determined by checking the90

floating condition,ρsw(zsl − zl ) :ρρswH>(zsl − zl ) :ρH≤floating ice,
(8)
grounded ice

91(where zl is the topography of the lithosphere surface), and the inflow into regions of floating ice results from the SIA solution

92

for the horizontal velocity field of adjacent grounded ice. A more sophisticated treatment of the sheet/shelf transition zone

93

(e.g., Pollard and DeConto, 2007, 2009) will be implemented in future work. The calving front is exposed to the hydrostatic

94

pressure of the adjacent sea water, and its position is determined by the simple, implicit calving criterion that floating ice

95

calves if its thickness is smaller than a threshold value of 50 m. This is considerably thinner than the thickness of the major

96

present-day Antarctic ice shelves; however, with a larger threshold of 200–250 m many of the smaller ice shelves (e.g., Riiser-97

Larsen, Shackleton) tend to erode during a simulation, and new ice shelves are largely prevented from forming. Sub-ice-shelfSato and Greve: Sensitivity experiments for Antarctica

98

5

melting, M , is parameterized by a method adopted from the GRISLI model (Ritz and others, 2001):

M=





Mgl




Mcs






 Mao

in the vicinity of the grounding line (see below),over the continental shelf, zl > −1600 m ,(9)over the abyssal ocean, zl ≤ −1600 m .

99

For the three parameters, we use the standard values Mgl = 2 m ice equiv. a−1 , Mcs = 0.2 m ice equiv. a−1 and Mao =100

10 m ice equiv. a−1 . For Mgl and Mcs , this is 2/3 of the values suggested in GRISLI (A. Quiquet, personal communication

1012010). The high value of Mao prevents ice shelves from extending over the abyssal ocean, where they are unstable due to the102

absence of embayments or islands/shoals that can serve as pinning points.

103

The model domain covers the entire area of Antarctica and the surrounding oceans, projected on a polar stereographic

104

grid with standard parallel 71◦ S and central meridian 0◦ E. Distortions due to this projection are accounted for as metric105

coefficients in all model equations except for the above-mentioned SSA equations for the horizontal velocities. The present

106

geometry (surface and basal topographies, ice thickness, equilibrated bedrock elevation, grounded vs. floating ice) is derived

107

from the “Antarctica Developmental Data Set” (Antarctica 5km dev1.0.nc) provided on the SeaRISE website, resampled to108

a horizontal resolution of 20 km. In the vertical direction, sigma coordinates are used; the cold ice column, the temperate

109

ice layer (if present) and the thermal boundary layer of the lithosphere are mapped separately to [0, 1] intervals. The cold110

ice column is discretized by 81 grid points (concentrated towards the base), and the temperate ice and lithosphere layers are

111

discretized each by 11 equidistant grid points. Time steps are 0.05–0.5 a for the topography and velocity, and 0.25–2 a for the112

thermodynamics.

113

The “vicinity of the grounding line” in Eq. (9) refers to floating-ice grid points that have a grounded-ice neighbour. In

114

addition, grounded-ice grid points with a floating-ice or ocean neighbour are assigned a basal melting rate M averaged over115

the grid point itself and the four neighbours,

M = rgr Mgr + (1 − rgr )Mgl

(10)

116

where rgr is the ratio of grounded-ice points among the five points (rgr = #(grounded-ice points)/5) and Mgr is the computed

117

basal melting rate of the central grounded-ice point. Consequently, the “vicinity of the grounding line” (influence of the118

prescribed parameter Mgl ) extends over twice the horizontal resolution, that is, 40 km. Due to this rather large, resolution-

119

dependent zone of influence, we refrain from choosing a value of the order of tens of metres per year for Mgl , even though such120

melting rates were measured near deep grounding lines (Rignot and Jacobs, 2002).

6

Sato and Greve: Sensitivity experiments for Antarctica

121

3. SEARISE EXPERIMENTS

122

3.1. Paleoclimatic spin-up

123

In order to obtain a suitable present-day configuration of the Antarctic ice sheet, it is desirable to carry out a paleoclimatic124

spin-up over at least a full glacial cycle. However, initial testing revealed that it is very difficult to reproduce the observed

125

geometry (in particular the distribution of grounded vs. floating ice) by an unconstrained, freely evolving simulation. For this126

reason, we carry out the spin-up simulation in four steps:

127(i) An initial relaxation run with freely evolving grounded ice topography over 100 a, starting from the present-day geometry

128

and isothermal conditions at −10◦ C everywhere, in order to avoid spurious noise in the computed velocity field (Calov,1291994). The grounding lines and floating ice topographies are kept fixed. The surface temperature and the sea level are130

those of today (see below), the surface mass balance is set to zero.

131

132

(ii) A steady-state run from 250 ka BP2 until 125 ka BP with the entire topography kept fixed over time. The surface
temperature is that of 125 ka BP, the surface mass balance is unspecified (due to the fixed topography).

133(iii) A transient run from 125 ka BP until today with the entire topography kept fixed over time in order to enforce a good fit

134

between the simulated and observed present-day topographies. The surface temperature varies over time (see below), the135

surface mass balance is unspecified.

136(iv) A transient run over 20 a with unconstrained evolution of the ice topography in order to avoid transition shocks in the137

beginning of the subsequent future climate experiments. The climatic forcing (surface temperature, surface mass balance)138

and the sea level are kept steady at today’s conditions.

139

Apart from the “frozen” geometry of the ice sheet in runs (ii) and (iii), the runs are conducted with the forcings suggested

140

by SeaRISE. The mean annual and mean summer surface temperatures (Tma and Tms , respectively) are decomposed into

141

present-day spatial distributions plus a purely time-dependent anomaly ∆T (t),present
Tma (h, φ, t) = Tma(h, φ) + ∆T (t) ,(11)
present
Tms (h, φ, t) = Tms(h, φ) + ∆T (t) ,2 Here, the notation ka BP means thousand calendar years before present.

Sato and Greve: Sensitivity experiments for Antarctica

7

where t is time and φ geographical latitude. The present-day parameterizations are taken from Fortuin and Oerlemans (1990),

Tma =





7.405 − 0.014285 h − 0.180 φ ,
36.689 − 0.005102 h − 0.725 φ ,






 49.642 − 0.943 φ ,for h > 1500 m ,for 200 m < h ≤ 1500 m ,(12)for h ≤ 200 m ,Tms = 16.81 − 0.00692 h − 0.27937 φ ,

(13)

142

where temperatures are in ◦ C, surface elevations in m AMSL and latitudes in ◦ S (counted positive). The time-dependent143

anomaly ∆T (t) results from the Vostok δD record converted to temperature with the relation by Petit and others (1999).

144

Precipitation, surface melting and sea level forcings are not required for runs (ii) and (iii) due to the “frozen” geometry

145

approach. In fact, in these runs the net surface mass balance (precipitation minus surface melting) is determined by the146

dynamical ice flow via the continuity equation, as a consequence of requiring the geometry to be unchanged. For run (iv),

147

precipitation is equated to the surface accumulation data by Arthern and others (2006), and surface melting is parameterized by

148

Reeh’s (1991) positive degree day (PDD) method, supplemented by the semi-analytical solution for the PDD integral by Calov

149

and Greve (2005). The PDD factors are (in ice equivalents) βice = 8 mm d−1 ◦ C−1 for ice melt and βsnow = 3 mm d−1 ◦ C−1 for150

snow melt (Ritz and others, 2001). Furthermore, the standard deviation of short-term, statistical air temperature fluctuations

151

is σ = 5◦ C, and the saturation factor for the formation of superimposed ice is chosen as Pmax = 0.6 (Reeh, 1991).

152

As for the geothermal heat flux, the SeaRISE “Antarctica Developmental Data Set” provides two different distributions by

153

Shapiro and Ritzwoller (2004) and Fox Maule and others (2005). We conduct the spin-up run with the one by Shapiro and154

Ritzwoller (2004), and, following a recommendation from the SeaRISE website, heat flux values are capped at 70 mW m−2 .

155

The simulated present-day configuration of the ice sheet (result of run (iv)) is used as initial condition for the future climate

156

experiments described in the following.

157

3.2. Future climate experiments

158

For the future climate experiments, we use the following set of SeaRISE experiments:

159

Experiment CTL: Constant climate control run; beginning at present (more precisely, the epoch 2004-1-1 0:00, correspond-160

ing to t = 0) and running for 500 a holding the climate steady to the present climate. We carry out this run twice, with161

the surface accumulation data by (i) Arthern and others (2006) and (ii) van de Berg and others (2006) (both provided in162

the “Antarctica Developmental Data Set” by SeaRISE).

163

Experiment E1a: Constant climate forcing running for 500 a, surface accumulation by Arthern and others (2006), uniform

164

sub-ice-shelf melting rate of Mgl = Mcs = Mao = 2 m ice equiv. a−1 (instead of the parameters given after Eq. (9)).

8Sato and Greve: Sensitivity experiments for Antarctica

165

Experiment E1b: Constant climate forcing running for 500 a, surface accumulation by Arthern and others (2006), uniform

166

sub-ice-shelf melting rate of Mgl = Mcs = Mao = 20 m ice equiv. a−1 (instead of the parameters given after Eq. (9)).

167

Experiment E1c: Constant climate forcing running for 500 a, surface accumulation by Arthern and others (2006), uniform

168

sub-ice-shelf melting rate of Mgl = Mcs = Mao = 200 m ice equiv. a−1 (instead of the parameters given after Eq. (9)).

169

Experiment E1z: Constant climate forcing running for 500 a, surface accumulation by Arthern and others (2006), zero170

sub-ice-shelf melting rate Mgl = Mcs = Mao = 0 (instead of the parameters given after Eq. (9)).

171

For all experiments, surface melting is parameterized by the PDD method (see Sect. 3.1), basal melting at grounded-ice grid172

points with a floating-ice or ocean neighbour is computed by the averaging scheme of Eq. (10), and the geothermal heat flux

173

is that by Shapiro and Ritzwoller (2004), capped at 70 mW m−2 (same as used for the paleoclimatic spin-up, see above).

174

The more recently defined “2011 Sensitivity Experiments” of the SeaRISE group, including investigations on changes of the

175

climate at the upper surface and changes of sliding at the subglacial interface, will be considered in future work.

176

4. RESULTS AND DISCUSSION177

The results of the paleoclimatic spin-up run (see Section 3.1) for the present are shown in Fig. 1. Due to the fixed-topography

178

approach, the surface topography (Fig. 1a) is of course very similar to the observed one (not shown). Small differences of the179

order of tens of meters arise as a consequence of the initial 100-year relaxation and the final 20-year run with freely evolving

180

topography. The surface velocity (Fig. 1b) shows the expected distribution with small velocities (< 10 m a−1 ) in the interior, a181

general speed-up towards the coast and the largest velocities of 1000 m a−1 and more for the ice shelves, while the limitations

182

of our approach (20-km resolution, fixed-topography spin-up with short transient run at the end) do not allow capturing all

183

fine structures adequately.

184

Basal temperatures (Fig. 1c) are at the pressure melting point in large areas in West Antarctica, while limited temperate

185base areas are found in East Antarctica in George V Land, Wilkes Land and several other places near the coast. At the deep186

ice core sites Vostok, Dome C, Dome F, Kohnen and Byrd, the spin-up run produces basal temperatures of −11.5◦ C, −11.6◦ C,

187−18.5◦ C, −11.4◦ C and −2.0◦ C (pressure melting point), respectively. Except for Byrd, these temperatures are significantly

188

below the pressure melting point, while the observed temperatures all reach the pressure melting point (Petit and others, 1999;

189

Motoyama, 2007; Parrenin and others, 2007; Wilhelms and others, 2007).

190

This discrepancy is probably due to shortcomings of the applied geothermal heat flux distribution by Shapiro and Ritzwoller

191(2004), which gives values in the range of 45 – 50 mW m−2 at Dome C, Dome F and Vostok. This is probably not sufficient192

to form a temperate base (Llubes and others, 2006, Fig. 3a). Seddik and others (2011) used a higher geothermal heat flux ofSato and Greve: Sensitivity experiments for Antarctica

9193

60 mW m−2 to produce a realistic borehole temperature in the vicinity of Dome F. Therefore, larger geothermal heat fluxes

194

would be needed to reproduce these observed melting conditions in East Antarctica. By contrast, the base of the West-Antarctic

195

Byrd ice core, which receives a geothermal heat flux of 70 mW m−2 in our simulation, reaches the melting point. A second

196

possible reason for the underprediction of basal temperatures is the fixed-topography approach of the spin-up run. Compared197

to the real, evolving Antarctic ice sheet, our simplified approach certainly influences the englacial heat transport (by both198

advection and conduction) to some extent, which affects basal temperatures.

199

Figure 2 shows the evolution of the grounded ice volume for the future climate control run with the two different accumu-

200

lation data mentioned in Section 3.2. In both cases, the ice volume increases by more than 1 m SLE over 500 years, and the201

accumulation data by van de Berg and others (2006) produce ∼ 0.32 m SLE more ice than those by Arthern and others (2006).

202

The reason is that the former accumulation rates are relatively larger than the latter, particularly around the margin of the ice

203sheet. By contrast, different geothermal flux settings [Shapiro and Ritzwoller (2004) vs. Fox Maule and others (2005)] cause

204

only slight changes (not shown). In order to minimize the imbalance of the ice sheet in the future climate control run, we

205

discard the accumulation data by van de Berg and others (2006) and use only those by Arthern and others (2006) for the future206

climate simulations with varied sub-ice-shelf melting rates (see Sect. 3.2). The imbalance for the initial time (t = 0 in Fig. 2,

207

corresponding to the year 2004) is ∼ 2.5 mm SLE a−1 , which indicates that the drainage of the ice sheet is underpredicted208(probably related to the underpredicted basal temperatures, see above).

209

The results of the future climate simulations with varied sub-ice-shelf melting rates are shown in Figs. 3 and 4. As expected,

210

the grounded ice volume (Fig. 3, shown relative to the control run because of the imbalance discussed above) becomes smaller

211

when the basal melting rates are larger. However, the effect is strongly non-linear and diminishes with increasing melting

212

rates. This agrees with the findings by Huybrechts and de Wolde (1999) for a similar series of simulations with smaller basal

213

melting rates of 0, 1, 3, 5 and 10 m a−1 . The spread of ice volumes across all experiments (from the zero to the maximum

214

melting case) is ∼ 0.65 m SLE after 100 years and ∼ 2.25 m SLE after 500 years. Excluding the most extreme experiment E1c

215

with 200 m a−1 basal melting and thus considering only the range from 0 to 20 m a−1 yields a spread of ∼ 2.0 m SLE after 500

216

years. This agrees well with the spread of ∼ 0.6 × 106 km3 (∼ 1.5 m SLE) reported by Huybrechts and de Wolde (1999) (their

217

Fig. 14) for the range from 0 to 10 m a−1 . So we see a significant, but not catastrophic impact of increased basal melting on

218

ice sheet decay.

219

Changes of the area of grounded and floating ice are shown in Fig. 4. In the control run (basal melting according to Eq. (9)),

220the grounded ice area drops by ∼ 0.25 × 106 km2 (∼ 2%) during the first 200 years, partly compensated by an increase of221

floating ice, and remains almost constant after that. In run E1a (2 m a−1 basal melting everywhere), grounded ice evolves

222

very similarly due to the identical melting rates in the vicinity of the grounding line, whereas two thirds of the floating ice

10

Sato and Greve: Sensitivity experiments for Antarctica

223

disappear over 500 years as a consequence of the 10 times larger melting rate over the continental shelf. In the more extreme

224

runs E1b and E1c (20 and 200 m a−1 basal melting everywhere, respectively), today’s floating ice disappears even stronger:

225

For E1b most of it melts within 100 years, the remaining area stabilizing at ∼ 0.15 × 106 km2 (∼ 10% of today’s value) due

226

to a balance between inflow from the ice sheet and basal melting, and for E1c almost all floating ice disappears within a few

227

years. This has a significant impact on the grounded ice area, which, for both E1b and E1c, shrinks continuously over the228

entire period of 500 years, by when the loss amounts to ∼ 0.5 × 106 km2 (∼ 4%) relative to the control run.

229

In general, under the influence of high basal melting, the Ross and Antarctic Peninsula ice shelves decay faster than the230

Filchner-Rønne and Amery ice shelves. In the E1a experiment (2 m a−1 basal melting), the Ross ice shelf loses more than 50%

231

of its area by t = 250 a and has almost disappeared by t = 500 a, while the Amery ice shelf does not lose much of its area.

232

The vulnerability of the Ross ice shelf results from the fact that it is surrounded by the Transantarctic Mountains, so that233

the incoming flux from the ice sheet is smaller than for other ice shelves. By contrast, the Amery ice shelf is largely stabilized

234

against increased basal melting by the large inflow from Lambert Glacier, which is the world’s largest glacier.

235

The grounded ice area is largest in run E1z (no basal melting), while the floating ice area is largest in the control run (basal

236

melting according to Eq. (9)). This is related to grounding line migration. Since there is no basal melting in run E1z, the237

grounding lines can advance freely all around the ice sheet, so that the grounded ice grows at the expense of the floating ice. By

238

contrast, basal melting favours grounding line retreat, so that in the control run floating ice grows at the expense of grounded

239

ice. Even the floating ice area of run E1a (2 m a−1 basal melting) is larger than that of run E1z (no basal melting) during

240

the first ∼ 80 years. However, later the retreat of the calving fronts due to basal-melt-induced thinning outweighs this effect,241

so that finally the floating ice area of run E1a becomes smaller than that of run E1z. For runs E1b and E1c (20 / 200 m a−1

242

basal melting) basal melting is so large that the loss due to thinning-induced calving rapidly overcompensates the gain due to243grounding line retreat, and thus no initial increase of the floating ice area is observed in the results.

244

5. CONCLUSION

245

A module for floating ice was added to the dynamic/thermodynamic ice sheet model SICOPOLIS. We have applied the new246

version of the model to the Antarctic ice sheet and carried out a paleoclimatic spin-up as well as a set of future climate runs

247

specified by the SeaRISE group. In the latter, the surface climate forcing was kept steady at present-day conditions, while248

the sub-ice-shelf melting rate was varied over a wide range. The simulations have revealed a significant, but not catastrophic

249

sensitivity of the ice sheet, the grounded ice volume showing a spread from the zero to the maximum melting case of ∼250

0.65 m SLE after 100 years and ∼ 2.25 m SLE after 500 years.

251

Some shortcomings of this study must be noted. Due to computing time limitations, we did the simulations of this study at

252

a horizontal resolution of 20 km. This resolution does not resolve well many fine structures in the vicinity of ice streams that

Sato and Greve: Sensitivity experiments for Antarctica

11

253

are important for the dynamics of the ice sheet as a whole. It is thus desirable to improve the resolution to at least 10 km.

254This is already on the way and will be the basis for the authors’ contribution to the SeaRISE community effort (publication

255

in preparation). A second aspect is the current treatment of grounding line dynamics. As explained in Section 2, the inflow

256

into regions of floating ice results from the SIA solution for the horizontal flow field of adjacent grounded ice. A smoother way

257

of dealing with this change of flow regimes would be to consider a transitional zone with shelfy stream dynamics (SSA with258

additional basal drag; MacAyeal (1989)). Following findings by Schoof (2007) and Pollard and DeConto (2009), for the rather

259

coarse resolutions typically applied in large-scale ice sheet models the inflow should be corrected by imposing a mass flux

260

constraint across the grounding line, which parameterizes ice velocities there as a function of local ice thickness and stresses.

261

Further, the simplified parameterization of basal melting (Eq. (9)) that does not distinguish between different regions and262

types of ice shelves, as well as the calving criterion in the form of a simple threshold value for the ice thickness, leave room263

for improvements. These points will be considered in future work.

264

ACKNOWLEDGEMENTS265

We thank R.A. Bindschadler, S. Nowicki and others for their efforts in the management of the SeaRISE project, J.V. Johnson266

and others for compiling and maintaining the SeaRISE datasets, A. Quiquet for information about the parameterization of267

sub-ice-shelf melting in the model GRISLI, and D. Pollard for discussions about the grounding line treatment. Comments of268

two anonymous reviewers and scientific editor D.M. Holland helped considerably to improve the manuscript. This study was

269

supported by a Grant-in-Aid for Scientific Research A (No. 22244058) from the Japan Society for the Promotion of Science

270(JSPS).

271

REFERENCES

272

Arthern, R. J., D. P. Winebrenner and D. G. Vaughan, 2006. Antarctic snow accumulation mapped using polarization of 4.3-cm wavelength

273

274

275

276

277

278

279

microwave emission, J. Geophys. Res., 111(D6), D06107.
van de Berg, W. J., M. R. van den Broeke, C. H. Reijmer and E. van Meijgaard, 2006. Reassessment of the Antarctic surface mass balanceusing calibrated output of a regional atmospheric climate model, J. Geophys. Res., 111(D11), D11104.
Bindschadler, R., 1983. The importance of pressurized subglacial water in separation and sliding at the glacier bed, J. Glaciol., 29(101),
3–19.
van den Broeke, M. R., J. Bamber, J. Lenaerts and E. Rignot, 2011. Ice sheets and sea level: Thinking outside the box, Surv. Geophys.,32(4-5), 495–505.

12Sato and Greve: Sensitivity experiments for Antarctica

280

Calov, R., 1994. Das thermomechanische Verhalten des Grönländischen Eisschildes unter der Wirkung verschiedener Klimaszenarien –281

Antworten eines theoretisch-numerischen Modells, Doctoral thesis, Department of Mechanics, Darmstadt University of Technology,

282

Germany.

283

284

Calov, R. and R. Greve, 2005. A semi-analytical solution for the positive degree-day model with stochastic temperature variations, J.
Glaciol., 51(172), 173–175.

285

Cuffey, K. M. and W. S. B. Paterson, 2010. The Physics of Glaciers, Elsevier, Amsterdam, The Netherlands etc., 4th ed.

286

Dunse, T., R. Greve, T. V. Schuler and J. O. Hagen, 2011. Permanent fast flow versus cyclic surge behaviour: numerical simulations of

287

288

289

290

291

292

293

the Austfonna ice cap, Svalbard, J. Glaciol., 57(202), 247–259.
Fortuin, J. P. F. and J. Oerlemans, 1990. Parameterization of the annual surface temperature and mass balance of Antarctica, Ann.
Glaciol., 14, 78–84.
Fox Maule, C., M. E. Purucker, N. Olsen and K. Mosegaard, 2005. Heat flux anomalies in Antarctica revealed by satellite magnetic data,
Science, 309(5733), 464–467.
Greve, R., 1997. Application of a polythermal three-dimensional ice sheet model to the Greenland ice sheet: Response to steady-state
and transient climate scenarios, J. Climate, 10(5), 901–918.

294

Greve, R., 2001. Glacial isostasy: Models for the response of the Earth to varying ice loads, Straughan, B., R. Greve, H. Ehrentraut and295

Y. Wang, eds., Continuum Mechanics and Applications in Geophysics and the Environment, Springer, Berlin, Germany etc., 307–325.

296

Greve, R. and H. Blatter, 2009. Dynamics of Ice Sheets and Glaciers, Springer, Berlin, Germany etc.

297

Hindmarsh, R. C. A. and E. Le Meur, 2001. Dynamical processes involved in the retreat of marine ice sheets, J. Glaciol., 47(157),

298

299

300

301

302

271–282.
Hutter, K., 1983. Theoretical Glaciology; Material Science of Ice and the Mechanics of Glaciers and Ice Sheets, D. Reidel Publishing
Company, Dordrecht, The Netherlands.
Huybrechts, P. and J. de Wolde, 1999. The dynamic response of the Antarctic and Greenland ice sheets to multiple-century climatic
warming, J. Climate, 12(8), 2169–2188.

303

IPCC, 2007. Summary for policymakers, Solomon, S., D. Qin, M. Manning, Z. Chen, M. Marquis, K. B. Averyt, M. Tignor and H. L.

304

Miller, eds., Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of305

the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK and New York, NY, USA, 1–18.

306

307

308

309

Joughin, I., E. Rignot, C. E. Rosanova, B. K. Lucchitta, J. Bohlander, J. L. Bamber and S. P. Gogenini, 2003. Timing of recent accelerationsof Pine Island Glacier, Antarctica, Geophys. Res. Lett., 30(13), 1706.
Le Meur, E. and P. Huybrechts, 1996. A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice
sheet during the last glacial cycle, Ann. Glaciol., 23, 309–317.

310

Lemke, P., J. Ren, R. B. Alley, I. Allison, J. Carrasco, G. Flato, Y. Fujii, G. Kaser, P. Mote, R. H. Thomas and T. Zhang, 2007.

311

Observations: Changes in snow, ice and frozen ground, Solomon, S., D. Qin, M. Manning, Z. Chen, M. Marquis, K. B. Averyt,

312

M. Tignor and H. L. Miller, eds., Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the FourthSato and Greve: Sensitivity experiments for Antarctica

13

313

Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, UK, and New York,

314

NY, USA, 337–383.

315

Lliboutry, L., 1968. General theory of subglacial cavitation and sliding of temperate glaciers, J. Glaciol., 7(49), 21–58.

316

Lliboutry, L. and P. Duval, 1985. Various isotropic and anisotropic ices found in glaciers and polar ice caps and their corresponding

317

318

319

320

321

rheologies, Ann. Geophys., 3(2), 207–224.
Llubes, M., C. Lanseau and F. Rémy, 2006. Relations between basal condition, subglacial hydrological networks and geothermal flux in
Antarctica, Earth Planet. Sci. Lett., 241(3-4), 655–662.
MacAyeal, D. R., 1989. Large-scale ice flow over a viscous basal sediment: theory and application to ice stream B, Antarctica, J. Geophys.
Res., 94(B4), 4071–4087.

322

Morland, L. W., 1984. Thermomechanical balances of ice sheet flows, Geophys. Astrophys. Fluid Dyn., 29, 237–266.

323

Morland, L. W., 1987. Unconfined ice-shelf flow, van der Veen, C. J. and J. Oerlemans, eds., Dynamics of the West Antarctic Ice Sheet,

324

D. Reidel Publishing Company, Dordrecht, The Netherlands, 99–116.

325

Motoyama, H., 2007. The second deep ice coring project at Dome Fuji, Antarctica, Sci. Drill., 5, 41–43.

326

Parrenin, F., G. Dreyfus, G. Durand, S. Fujita, O. Gagliardini, F. Gillet, J. Jouzel, K. Kawamura, N. Lhomme, V. Masson-Delmotte,

327

C. Ritz, J. Schwander, H. Shoji, R. Uemura, O. Watanabe and N. Yoshida, 2007. 1-D-ice flow modelling at EPICA Dome C and Dome

328

Fuji, East Antarctica, Clim. Past, 3(2), 243–259.

329

Petit, J. R., J. Jouzel, D. Raynaud, N. I. Barkov, J. M. Barnola, I. Basile, M. Bender, J. Chappellaz, M. Davis, G. Delaygue, M. Delmotte,

330

V. M. Kotlyakov, M. Legrand, V. Y. Lipenkov, C. Lorius, L. Pepin, C. Ritz, E. Saltzman and M. Stievenard, 1999. Climate and331

atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica, Nature, 399(6735), 429–436.

332

Pollard, D. and R. M. DeConto, 2007. A coupled ice-sheet/ice-shelf/sediment model applied to a marine-margin flowline: forced and333

unforced variations, Hambrey, M. J., P. Christoffersen, N. F. Glasser and B. Hubbard, eds., Glacial Sedimentary Processes and334

Products, Blackwell Publ., International Association of Sedimentologists Special Publication No. 39, 37–52.

335

336

Pollard, D. and R. M. DeConto, 2009. Modelling West Antarctic ice sheet growth and collapse through the past five million years, Nature,458(7236), 329–332.

337

Reeh, N., 1991. Parameterization of melt rate and surface temperature on the Greenland ice sheet, Polarforsch., 59(3), 113–128.

338

Rignot, E., G. Casassa, P. Gogineni, W. Krabill, A. Rivera and R. Thomas, 2004. Accelerated ice discharge from the Antarctic Peninsula

339

340

341

342

343

344

345

following the collapse of Larsen B ice shelf, Geophys. Res. Lett., 31(18), L18401.
Rignot, E. and S. S. Jacobs, 2002. Rapid bottom melting widespread near Antarctic ice sheet grounding lines, Science, 296(5575),
2020–2023.
Rignot, E., I. Velicogna, M. R. van den Broeke, A. Monaghan and J. Lenaerts, 2011. Acceleration of the contribution of the Greenlandand Antarctic ice sheets to sea level rise, Geophys. Res. Lett., 38, L05503.
Ritz, C., V. Rommelaere and C. Dumas, 2001. Modeling the evolution of Antarctic ice sheet over the last 420,000 years: Implications foraltitude changes in the Vostok region, J. Geophys. Res., 106(D23), 31943–31964.

14

Sato and Greve: Sensitivity experiments for Antarctica

346

Schoof, C., 2007. Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res., 112(F3), F03S28.

347

Seddik, H., R. Greve, T. Zwinger and L. Placidi, 2011. A full Stokes ice flow model for the vicinity of Dome Fuji, Antarctica, with induced

348

349

350

351

352

anisotropy and fabric evolution, The Cryosphere, 5(2), 495–508.
Shapiro, N. M. and M. H. Ritzwoller, 2004. Inferring surface heat flux distributions guided by a global seismic model: particular application
to Antarctica, Earth Planet. Sci. Lett., 223(1-2), 213–224.
Shepherd, A., D. J. Wingham, J. A. D. Mansley and H. F. J. Corr, 2001. Inland thinning of Pine Island Glacier, West Antarctica, Science,
291(5505), 862–864.

353

Weertman, J., 1964. The theory of glacier sliding, J. Glaciol., 5(39), 287–303.

354

Weis, M., R. Greve and K. Hutter, 1999. Theory of shallow ice shelves, Cont. Mech. Thermodyn., 11(1), 15–50.

355

Wilhelms, F., S. G. Sheldon, I. Hamann and S. Kipfstuhl, 2007. Implications for and findings from deep ice core drillings – an example: the356

ultimate tensile strength of ice at high strain rates, Kuhs, W. F., ed., Physics and Chemistry of Ice, The Royal Society of Chemistry,

357

Cambridge, UK, Special Publication No. 311, 635–639.

358

359

Zwally, H. J. and M. B. Giovinetto, 2011. Overview and assessment of Antarctic ice-sheet mass balance estimates: 1992-2009, Surv.
Geophys., 32(4-5), 351–376.

Sato and Greve: Sensitivity experiments for Antarctica

Quantity

Value

Gravitational acceleration, g

9.81 m s−2

Density of ice, ρ

910 kg m−3

Density of sea water, ρsw1028 kg m−3

Power law exponent, n3Residual stress, σ0

10 kPa15

2.847 × 10−13 s−1 Pa−3

(T 0 ≤ − 10◦ C)

2.356 × 10−2 s−1 Pa−3

(T 0 ≥ − 10◦ C)

Pre-exponential constant, A0

60 kJ mol−1(T 0 ≤ − 10◦ C)

115 kJ mol−1(T 0 ≥ − 10◦ C)

Activation energy, Q

Flow enhancement factor, E

5 / 1?

Melting point at atmosphericpressure, T0

273.15 K = 0◦ C

Heat conductivity of ice, κ

9.828 e−0.0057 T [K] W m−1 K−1Specific heat of ice, c(146.3 + 7.253 T [K]) J kg−1 K−1

Latent heat of ice, L

335 kJ kg−1

Clausius-Clapeyron gradient, β

8.7 × 10−4 K m−1

Universal gas constant, R

8.314 J mol−1 K−1Isostatic time lag, τiso3000 aAsthenosphere density, ρa3300 kg m−3

Flexural stiffness of thelithosphere, Kl

1025 N m

Density × specific heat of thelithosphere, ρr cr2000 kJ m−3 K−1

Heat conductivity of thelithosphere, κr

3 W m−1 K−1Thickness of the thermal boundary
layer of the lithosphere, Hr

5 km

Table 1. Physical parameters used for the simulations of this study. ? : E = 5 for grounded ice, E = 1 for floating ice.

Sato and Greve: Sensitivity experiments for Antarctica

50

km

S

4.0

Surface topography

S
50

3000

0

W
30

a

30
E

16

3.5
2000

E
60

1
2

3.0

3

1

1

2.5
2

3

3

2

1

1

0

2.0

90E
1

y (km)

1000

1.5

2

2

3

-1000

2

1.0

-2000

-1000

0

1000

2000

0.5

0E

0.0

3000

m/a
3000

3000

Surface velocity

50

b

12

1

-2000

S

1000

2000

6

0E

300
100

y (km)

1000

30
0

90E

10
3

-1000

1
12

-2000

-2000

0

1000

2000

0.3
0

3000

C
3000

0

Basal temperature

S

50

c

-1000

0E

E

60

2000

-4

y (km)

1000

0

90E

-8

-1000
-12
12

-2000

-2000

-1000

0

1000

2000

3000

0E

-16

x (km)

Fig. 1. Present-day configuration of the Antarctic ice sheet computed by the paleoclimatic spin-up. (a) Surface topography (in km AMSL,
contour spacing is 250 m), (b) surface velocity (in m a−1 ), (c) basal temperature (in ◦ C relative to the pressure melting point).

Sato and Greve: Sensitivity experiments for Antarctica

17

61.5
CTL(AE)V (m SLE)61

CTL(BE)

60.5

60

59.5
0

100

200

t (a)

300

400

500

Fig. 2. Simulated grounded ice volume, V (expressed as SLE), for experiment CTL (constant climate control run), carried out alternatively
with the surface accumulation data by Arthern and others (2006) (marked AE) and van de Berg and others (2006) (marked BE). Notethat t = 0 corresponds to the year 2004.

18

Sato and Greve: Sensitivity experiments for Antarctica1∆V (m SLE)

0.5

0

−0.5

−1

CTL
2 m/a
20 m/a−1.5

200 m/aNo melt

−2
0

100

200

t (a)

300

400

500

Fig. 3. Simulated changes in grounded ice volume, ∆V , relative to the control run CTL (expressed as SLE), for experiments E1a (constantclimate, sub-ice-shelf melting 2 m ice equiv. a−1 ), E1b (constant climate, sub-ice-shelf melting 20 m ice equiv. a−1 ), E1c (constant climate,sub-ice-shelf melting 200 m ice equiv. a−1 ) and E1z (constant climate, zero sub-ice-shelf melting). Surface accumulation by Arthern andothers (2006). Note that t = 0 corresponds to the year 2004.

Sato and Greve: Sensitivity experiments for Antarctica

a 13

19

Grounded ice

A (106 km2)

12.5

12
CTL

11.5

2 m/a
20 m/a200 m/aNo meltb2

Floating ice

A (106 km2)

1.5

1
CTL
2 m/a
0.5

20 m/a200 m/aNo melt

0
0

100

200

t (a)

300

400

500

Fig. 4. Simulated ice area, A, for experiments E1a (constant climate, sub-ice-shelf melting 2 m ice equiv. a−1 ), E1b (constant climate,sub-ice-shelf melting 20 m ice equiv. a−1 ), E1c (constant climate, sub-ice-shelf melting 200 m ice equiv. a−1 ) and E1z (constant climate,zero sub-ice-shelf melting). Panel (a): grounded ice area, panel (b): floating ice area. Note that t = 0 corresponds to the year 2004.