Simulating ice thickness and velocity evolution of Upernavik Isstrøm 1849–2012 by forcing prescribed terminus positions in ISSM

Abstract

Tidewater glacier velocity and mass balance are known to be highly responsive to terminus position change. Yet it remains challenging for ice flow models to reproduce observed ice margin changes. Here, using the Ice Sheet System Model ISSM;, we simulate the ice velocity and thickness changes of Upernavik Isstrøm (north-western Greenland) by prescribing a collection of 27 observed terminus positions spanning 164 years (1849–2012). The simulation shows increased ice velocity during the 1930s, the late 1970s and between 1995 and 2012 when terminus retreat was observed along with negative surface mass balance anomalies. Three distinct mass balance states are evident in the reconstruction: (1849–1932) with near zero mass balance, (1932–1992) with ice mass loss dominated by ice dynamical flow, and (1998–2012), when increased retreat and negative surface mass balance anomalies led to mass loss that was twice that of any earlier period. Over the multi-decadal simulation, mass loss was dominated by thinning and acceleration responsible for 70 $\mathrm{\%}$ of the total mass loss induced by prescribed change in terminus position. The remaining 30 $\mathrm{\%}$ of the total ice mass loss resulted directly from prescribed terminus retreat and decreasing surface mass balance. Although the method can not explain the cause of glacier retreat, it enables the reconstruction of ice flow and geometry during 1849–2012. Given annual or seasonal observed terminus front positions, this method could be a useful tool for evaluating simulations investigating the effect of calving laws.

Introduction

In recent decades, glaciers terminating into the ocean (tidewater glaciers) have exhibited widespread thinning and velocity acceleration e.g.. Increased air and ocean temperatures induce increased surface melt rates and frontal retreat , represented by submarine melt and iceberg calving. The Greenland ice sheet has more than 240 tidewater glacier outlets and its mass balance is highly affected by changes in tidewater glacier discharge . Greenland's ice mass changes have dominated global sea level contributions of the past two decades e.g.. Sea level projections rely on models to estimate discharge and Greenland's contribution to sea level, which are coming more into line with observations . However, accurate simulation of terminus position remains a major challenge IPCC, 2013, Sect. 13.

Tidewater glacier retreat occurs due to calving and submarine melt . Yet, no universal calving law exists and model approaches either (1) focus on the development and performance of a particular calving law e.g.; (2) simplify glacier simulations using flow line or flow band models e.g., neglecting, for example, across-flow stresses or (3) determine glacier terminus changes on ice particle scale and are thereby not well suited for long-term studies .

Upernavik Isstrøm (UI), a set of West Greenland tidewater glaciers, has been the focus of several observational studies. compiled historical records of UI terminus positions between 1849 and 1953, concluding that terminus retreat had increased starting in the 1930s. Observed periods of increased UI terminus retreat in 1931 to 1946, in the late 1990s and in 2005–2009 correlate with elevated air temperatures . Two periods of increased dynamically driven ice loss that took place in 1985–1993 and 2005–2010 were responsible for 79 $\mathrm{\%}$ of the ice mass loss during 1985–2010 .

Previous studies either simulate tidewater glacier retreat with ice flow models or discuss observed terminus changes and its implications for tidewater glaciers. In this study, we combine observations and ice flow models by using observed terminus positions in the Ice Sheet System Model ISSM; to simulate Upernavik's glacial system evolution from 1849, near the end of the Little Ice Age, to 2012. We reconstruct the 1849 ice surface elevation and force ISSM glacier terminus retreat with 27 observed terminus positions.

This study does not aim to simulate physically caused retreat, instead we evaluate the effects of changing termini on UI's ice surface elevation and velocity. We (1) investigate whether prescribed terminus change produces a realistic thinning and velocity history; (2) compare simulated mass loss, surface elevation and velocity changes with 1985–2012 observations; and (3) correlate the calculated dynamic ice loss with observational studies. ISSM produces a weekly reconstruction of UI ice thickness and surface velocity from 1849–2012.

Area and data

UI has a catchment area of $\sim \mathrm{64}$ 700 ${\mathrm{km}}^{\mathrm{2}}$, terminating into several tidewater glaciers. We focus on the three main glaciers and denote them UI-1, UI-2 and UI-3 from north to south (Fig. ). Historically, the three glaciers shared the same terminus between 1849 and 1931 Fig. ;. In the 1930s the glaciers separated in two, UI-1/UI-2 and UI-3. UI-1 and UI-2 decoupled after 1966. Historical front positions (Fig. ) were collected from several sources: 1849–1953 historical records;, 1966–1975 satellite images;, 1985–1996 aerial photographs (1985) and satellite images;, 1999 to 2012 satellite images;.

For initialization and evaluation of the model we use data from different studies, described in Table .

Upernavik Isstrøm's observed margin front positions between 1849 and 2012 (lines) and trimline positions yellow dots;. The background image is from Landsat 8 (September 2013). Inset is the location and shape of the Upernavik catchment (red area), determined by 2008/09 surface velocity from , which define the model domain.

[Figure omitted. See PDF]

Data for initialising and evaluating the simulation.

Datum	Source	Description
Bed topography		Derived with mass conservation approach, extended with bathymetry measurements
Bathymetry measurement		2012 NASA project, led by Eric Rignot and Todd Dupont
Bathymetry measurement	;	NASA project Oceans Melting Greenland OMG
Bathymetry measurement		Ship-based single point echo sounders
Trimline points		Little Ice Age maximum extent (Fig. )
Surface mass balance (SMB)		Monthly data, covering 1840–2012
1985 Digital elevation model(DEM)		Based on aerial photographs, 25 $\mathrm{m}$ resolution
2005 DEM		Greenland Ice Sheet Mapping Project (GIMP), 30 $\mathrm{m}$ resolution
2012 DEM		ArcticDEM, 2–10 $\mathrm{m}$ resolution
Ice surface velocity		Winter 2008/09
Ice surface velocity	http://esa-icesheets-greenland-cci.org/ (described in )	Provided by ESA project Climate Change Initiative (CCI) Greenland Ice Sheet in winters between 1991/92 and 2008/09
Ice surface velocity		Provided by MEaSUREs, in the winters 2000/01, 2007/08 and 2009/10
Ice surface elevation	;	IceBridge ATM; UI-1 in 2009–2012 and UI-3 in 1994, 1999, 2002, 2009, 2010, 2012
Mass change	;	Provided by the Jet Propulsion Laboratory (version: JPL RL05M GRACE mascon solution); suitable for regional (300 $\mathrm{km}$ scale) ice sheet mass change comparisons

Ice flow model

We use the Ice Sheet System Model ISSM;, a finite-element thermomechanical ice flow model. Ice flow is calculated applying the Shelfy-Stream Approximation SSA;, that integrates vertically averaged ice properties (e.g. ice rheology, thickness, velocity) and neglects vertical shear stresses. The SSA is well suited for fast-flowing glaciers like Upernavik, where the ice flow is primarily driven by basal sliding.

Ice viscosity follows Glen's flow law . The initial viscosity is taken from Table 3.4 in p. 75, assuming ice temperature of $-\mathrm{5}$ ${}^{\circ }\mathrm{C}$ and will be refined in Sect. .

A Budd-like friction law is applied to all grounded ice: $${\displaystyle }{\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=-C^{\mathrm{2}}N{v}_{\mathrm{b}},$$ where $v_{\mathrm{b}}$ is the basal velocity. The effective pressure $N$ is assumed to be $N=g\left({\mathit{\rho }}_{\mathrm{i}}H\right(x,y,t)+{\mathit{\rho }}_{\mathrm{w}}b(x,y\left)\right)$ where $H$ is the ice thickness at the current time step, $b$ is the bed elevation with respect to sea level, $g$ is the gravity, ${\mathit{\rho }}_{\mathrm{i}}$ and ${\mathit{\rho }}_{\mathrm{w}}$ are the densities of ice and water respectively. The friction coefficient $C$ is variable in space, but constant in time (Fig. 1, Supplement). Perfect sliding is assumed on floating ice.

Steps for model initialization.

Step	Input	Output
Relaxation 1	GIMP extended to 1849 terminus position,ice viscosity (initial guess),basal friction (initial guess)	Reconstructed 1849 ice thickness and velocity
Thermal	Ice thickness and velocity from relaxation 1	Improved ice viscosity
Inversion	Surface velocity from relaxation 1,ice viscosity from thermal	Inverted basal friction
Relaxation 2	Ice thickness from relaxation 1,ice viscosity from thermal,basal friction from inversion	Steady state ice thickness and velocity

The model domain is set to the Upernavik catchment, which is defined by the flow direction given by the 2008/09 surface velocity from (red area in Fig. ). We use an adaptive mesh that has a resolution varying between 300–800 $\mathrm{m}$ in the area of observed terminus changes and 12 $\mathrm{km}$ near the ice divide, resulting in about 17 000 mesh elements. Resolution increases with larger changes in ice velocity or bedrock topography and decreases stepwise with distance from the front.

We impose hydrostatic pressure at the terminus and keep the ice velocity and surface elevation constant at the inland boundary. No submarine frontal melt or calving rates are applied, since the study aims to simulate ice velocity and thickness changes caused by observational prescribed terminus changes. The ice is allowed to float depending on a hydrostatic criterion .

Model initialization

Since starting the simulation in 1849 extends the present day ice extent by 356 ${\mathrm{km}}^{\mathrm{2}}$, model initialization requires reconstruction of the ice surface elevation in the extended area. To initialise the model we thus reconstruct the 1849 ice surface elevation, as described below. Over the present day ice covered area, the initial ice surface is given by the 2005 ice surface elevation GIMP;. At the 1849 marine terminus given by , the ice surface elevation is set to 70 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ consistent with marine termini in the area, based on IceBridge data . Trimline data points Fig. ; mark the 1849 surface elevation and ice extent on the bedrock along the fjords. In the remaining area, the ice surface elevation is interpolated linearly being constrained to a minimum elevation of 40 $\mathrm{m}$. The ice thickness is set to floatation height or to the maximum thickness, defined through the initialised ice surface elevation and bed topography. The ice surface velocity is resolved performing a stress balance solution.

During the relaxation, the reconstructed glacier area thickened by 50–400 $\mathrm{m}$ from the UI-1/UI-2 1966 terminus position reaching 40 $\mathrm{km}$ upstream, while the ice surface velocity slowed down by a maximum of 2.5 $\mathrm{km}{\mathrm{yr}}^{-\mathrm{1}}$. Along UI-3, the glacier thickness changed by $\pm \mathrm{150}$ $\mathrm{m}$ and ice surface velocity decreased by 500–1500 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ (see Supplement).

As we are interested in determining how the model geometry and velocity react to the prescribed terminus change and not internal model instability, we relax the model prior to the transient run, bringing ice surface elevation and velocity into equilibrium (following ). Equilibrating model geometry and velocity requires constant forcing, i.e. a stable SMB. The SMB at Upernavik is found to be stable in 1854–1900 and 1964–1990. The mean 1854–1900 SMB value is used for equilibrating the model for 1849 conditions and 1964–1990 is set as the SMB reference period to evaluate simulated mass balance.

We perform two relaxation runs stepwise (Table ), keeping SMB constant to the 1854–1900 mean value . The first relaxation provides reconstructed 1849 ice thickness, given the GIMP surface elevation extended to the 1849 terminus. Thus, in the first relaxation basal friction is based on the assumption that driving stress is equal to basal stress at any given point using the initial geometry.

Given computed ice velocity and thickness from the first relaxation, ice viscosity and basal friction can be redefined. The ice viscosity is calculated by extruding the model with 15 layers and solving for the thermal steady state based on forcing the surface with 1854–1900 UI mean surface air temperature . The basal friction coefficient is constant in time, but varies in space, and is calculated by an adjoint-based inversion, following and , given the updated ice viscosity from the thermal steady state simulation.

The second relaxation runs for 5000 years until ice velocity and thickness are equilibrated, provided with ice thickness from the first initialization, simulated ice viscosity and inverted basal friction. The end state of this relaxation provides the initial values of simulated ice surface elevation and surface velocity for the 1849–2012 simulations.

Simulation setup

We run the following two different model simulations: (1) a control run ISSM${}_{\text{control}}$, forced only by monthly SMB using a fixed terminus at the observed 1849 ice margin and (2) a prescribed terminus change simulation ISSM${}_{\text{PT}}$, forced by the same monthly SMB and observed calving front positions. ISSM${}_{\text{control}}$ serves to estimate the ice mass, velocity and thickness changes that are simulated without prescribed terminus change.

The prescribed terminus position change in ISSM${}_{\text{PT}}$ is implemented through a levelset-based method and performed in July of the observation year, according to observed terminus positions (Fig. ). The highest surface air temperatures and melt rates on UI are observed in July , increasing the likelihood of terminus retreat . We introduce 20 additional calving front positions, created through linear interpolation between the observed termini positions and constrained by the mesh resolution. The additional calving fronts are prescribed at the halfway points in time between observations and aim to improve realistic simulation behaviour by splitting 20 large ice area changes induced by the prescribed terminus changes into smaller areas within shorter time periods (dotted lines, Fig. ).

(a) Simulated cumulative ice mass in Gt. ISSM${}_{\text{PT}}$ changes are shown in red; control run changes in purple. (b) The blue curve illustrates simulated annual change in ice mass for ISSM${}_{\text{PT}}$. The black bars indicate the ice mass that is removed due to $\mathrm{\Delta }$SMB and prescribed changes of the terminus position. The area between the blue line and black bars corresponds to ice mass changes caused by changes in the ice dynamics that where not prescribed. The green outline marks the portion of mass change due to $\mathrm{\Delta }$SMB, and the orange outline the share of prescribed terminus change.

[Figure omitted. See PDF]

Within the prescribed ice area, the grounding line is evolving freely and floating tongue formation is thereby allowed.

The simulation evaluation time step is set to 73 $\mathrm{h}$, constrained by the Courant–Friedrich–Lewy condition , ensuring numerical stability when solving the ice flow equations at each time step.

Simulated ice mass changes from anomalies (relative to 1964–1990 mean values) for the simulation ISSM${}_{\text{PT}}$. The background is highlighted in yellow for periods of time where $\mathrm{\Delta }$SMB controlled MB, blue is where ice mass loss is driven by $\mathrm{\Delta }$DIL and red, where $\mathrm{\Delta }$SMB and $\mathrm{\Delta }$DIL have an equally increased influence on the MB. Prescribed termini changes are marked with dashed (observations) and dotted (interpolations) lines.

[Figure omitted. See PDF]

Results and comparison

During the simulation, most of the ice thickness and velocity changes occur near the central flow lines of UI-1, UI-2 and UI-3. Simulated changes in ice thickness and velocity in the majority of the model domain (more than 70 $\mathrm{km}$ inland from the 2012 terminus or 5 $\mathrm{km}$ away from the central flow lines of the three glaciers) are below 5 $\mathrm{\%}$, corresponding to changes of 20 $\mathrm{m}$ and 10 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ over 164 simulation years. Thus, in the following, we present relative and absolute changes in ice velocity and thickness along the central flow lines of UI-1, UI-2 and UI-3 from the 2012 terminus reaching 30 $\mathrm{km}$ upstream (Fig. ).

Model comparison

Between 1849 and 2012, ISSM${}_{\text{control}}$ shows less than 7 $\mathrm{\%}$ thinning and 5 $\mathrm{\%}$ acceleration, simulating a change in velocity less than 120 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ and a thinning less than 30 $\mathrm{m}$ along the central flow lines for the entire period. In contrast, ISSM${}_{\text{PT}}$ produces a thinning between 20 $\mathrm{\%}$ along the flow lines and up to 60 $\mathrm{\%}$ in the area between 2012 terminus and 70 $\mathrm{km}$ upstream in 1849–2012, corresponding to thinning between 100 and 450 $\mathrm{m}$ along the flow lines. The average ice surface velocity increase along UI-1 and UI-2 is 180, and it is 47 $\mathrm{\%}$ on UI-3. Cumulative ice mass loss over the simulation period of the entire model domain (converted from modelled water equivalent assuming 917 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ ice density) was by the end of the model simulation $-$50 $\mathrm{Gt}$ for ISSM${}_{\text{control}}$ and $-$585 $\mathrm{Gt}$ for ISSM${}_{\text{PT}}$ (Fig. ). 99 $\mathrm{\%}$ of simulated ISSM${}_{\text{control}}$ mass loss was prescribed by SMB anomalies while 30 $\mathrm{\%}$ of total ice mass loss simulated by ISSM${}_{\text{PT}}$ was prescribed, with SMB anomalies accounting for 9 $\mathrm{\%}$ ($-$50 $\mathrm{Gt}$) and prescribed terminus position change contributed 21 $\mathrm{\%}$ ($-$121 $\mathrm{Gt}$). Thus, 70 $\mathrm{\%}$ of by ISSM${}_{\text{PT}}$ simulated mass loss is caused by thinning and acceleration. The following subsections describe ISSM${}_{\text{PT}}$ results in more detail.

Mass balance

In the following section we focus on the simulated mass balance (MB) through the model runs (see cumulated mass change in Fig. ). For marine terminating glaciers, mass balance can be attributed to either changes in SMB or changes in dynamic ice loss (DIL). A tidewater glacier is in equilibrium, when SMB and DIL are in balance. Deviations in SMB and DIL change the glacier and its stability is hereafter referred to as anomalies $\mathrm{\Delta }$SMB and $\mathrm{\Delta }$DIL. SMB is a model input and $\mathrm{\Delta }$SMB are calculated relative to the mean value of the stable UI period 1964–1990 SMB. $\mathrm{\Delta }$DIL is calculated as the residual between the simulated MB and $\mathrm{\Delta }$SMB.

The simulated annual MB for the UI catchment (Fig. ) is positive from 1849 to 1920. In this period, the MB from the ISSM${}_{\text{PT}}$ and ISSM${}_{\text{control}}$ are similar due to very few and small terminus changes (Fig. ) and MB is thus dominated by $\mathrm{\Delta }$SMB. Anomalies in DIL (Fig. ) are evident by small ($-$0.5 to $-$4 $\mathrm{Gt}$) peaks that coincide with prescribed terminus retreat. After 1920, the MB becomes negative, except in 1996, when $\mathrm{\Delta }$SMB has a peak (8 $\mathrm{Gt}$), which is attributed to a high winter accumulation . Figure  highlights three periods in MB trends: (1) 1849–1932, when MB is near equilibrium, (2) from 1932 to 1992, when the negative MB is driven by $\mathrm{\Delta }$DIL, and (3) 1998–2012, when SMB and DIL both have high negative anomalies and the total mass loss each year was twice as high as any year before.

Observed vs. simulated ice mass changes (with ISSM${}_{\text{PT}}$).

	${}^{\mathrm{a}}$
	1985–2002/05	2002/05–2010	2000–2005	2006–2008	2009–2011
Total observed ice mass changes (Gt)	$-\mathrm{32}\pm \mathrm{9}$	$-\mathrm{17}\pm \mathrm{10}$	$-\mathrm{6}\pm \mathrm{20}$	$-\mathrm{25}\pm \mathrm{14}$	$-\mathrm{39}\pm \mathrm{17}$
Total simulated ice mass changes (Gt)	$-{\mathrm{37}}^{\mathrm{b}}$	$-{\mathrm{32}}^{\mathrm{b}}$	$-\mathrm{48}$	$-\mathrm{41}$	$-\mathrm{44}$
Observed dynamic ice loss (Gt)	$\mathrm{29}\pm \mathrm{9}(\mathrm{90}\mathrm{\%}{)}^{\mathrm{c}}$	$\mathrm{12}\pm \mathrm{11}(\mathrm{80}\mathrm{\%}{)}^{\mathrm{c}}$	$\mathrm{5}\pm \mathrm{10}(\mathrm{80}\mathrm{\%}{)}^{\mathrm{c}}$	$\mathrm{16}\pm \mathrm{4}(\mathrm{62}\mathrm{\%}{)}^{\mathrm{c}}$	$\mathrm{27}\pm \mathrm{4}(\mathrm{68}\mathrm{\%}{)}^{\mathrm{c}}$
Simulated dynamic ice loss (Gt)	${\mathrm{32}}^{\mathrm{b}}(\mathrm{86}\mathrm{\%}{)}^{\mathrm{d}}$	${\mathrm{26}}^{\mathrm{b}}(\mathrm{81}\mathrm{\%}{)}^{\mathrm{d}}$	$\mathrm{40}(\mathrm{83}\mathrm{\%}{)}^{\mathrm{d}}$	$\mathrm{24}(\mathrm{59}\mathrm{\%}{)}^{\mathrm{d}}$	$\mathrm{28}(\mathrm{64}\mathrm{\%}{)}^{\mathrm{d}}$

${}^{\mathrm{a}}$ converted from ${\mathrm{km}}^{\mathrm{3}}$ to Gt ice equivalent; ${}^{\mathrm{b}}$ results from 2002/05 as mean values of that time; ${}^{\mathrm{c}}$ average percentage of total mass change induced by dynamic ice loss; ${}^{\mathrm{d}}$ percentage of total mass change that is induced by dynamic ice loss

and measure surface elevation changes from aerial photographs, satellites and digital elevation models between 1985 and 2010. These yield a total mass change during different time periods and, congruent to our calculations, $\mathrm{\Delta }$DIL is estimated as the residual of mass change and $\mathrm{\Delta }$SMB. Both studies refer to different areas within the UI catchment. Table  presents a comparison of the observed mass changes and our simulation results, recalculated for the particular areas. Due to sparse data coverage, combine surface elevation measurements acquired between 2002 and 2005 to quantify elevation changes and refer to this period as 2002/05. The average of simulated ice mass loss between 2002 and 2005 is taken for comparison with the 2002/05 observations from .

Simulated total ice mass changes in 1985–2002/05 and 2006–2011 correspond with observed ice mass changes from and (Table ). Additionally, the DIL during 2000 to 2005 makes up 83 $\mathrm{\%}$ of the mass change and in 2006–2011 this percentage is reduced to 64 $\mathrm{\%}$, in agreement with and . In 2000–2005, however, simulated total mass changes are 81 $\mathrm{\%}$ larger than the maximum of what is observed. A comparison with GRACE that measures gravity field variations from which mass change is computed, shows equivalent seasonal mass loss fluctuations in summer and mass gain in winter with an overall negative trend. The simulated mass change rates resemble 98 $\mathrm{\%}$ of GRACE's rate (see Supplement).

Ice thickness

ISSM${}_{\text{PT}}$ simulates 10–80 $\mathrm{\%}$ thinning from 1849 to 2012 over an area reaching 70 $\mathrm{km}$ upstream from the 2012 terminus (see Supplement). Transient surface elevation changes along the central flow lines of UI-1, UI-2 and UI-3 are visualised in the supplementary Fig. 8 and movie01 (Supplement). The model simulation shows increased surface lowering in the time periods 1930/40, 1970/80 and from 2000 onwards.

Observed vs. simulated ice thickness along flight lines IceBridge surface elevation data; over UI-1 and UI-3. Stars mark mean values between 0 and 5 $\mathrm{km}$ from the 2012 terminus, dots refer to mean values 5–10 $\mathrm{km}$ upstream from the 2012 terminus. Flight lines over UI-1 are available for the years 2009, 2010, 2011, 2012 and over UI-3 in the years 1994, 1999, 2002, 2009, 2010, 2012.

[Figure omitted. See PDF]

To evaluate simulated ice thickness, we compare simulation results with the residual ice thickness obtained from observed surface elevation data and the bed topography from , that is used in the simulation setup. We refer to the Supplement for illustrations of spatial comparisons between simulation results and observations. Simulated thickness of the UI catchment in 2005 lies within $\pm \mathrm{30}$ $\mathrm{\%}$ of observed thickness (GIMP), except in the shear margin regions of UI-1, where simulated ice thickness is too high by up to 160 $\mathrm{\%}$ of observations. A comparison of absolute ice thickness in 2005 shows up to 200 $\mathrm{m}$ lower simulated thickness than observed, apart from the shear zones of UI-1, where the ice is up to 200 $\mathrm{m}$ thicker than observed. Differences between the simulated 2012 ice thickness and observations (ArcticDEM) show the same pattern with less difference in the UI-1 shear zone. The 1985 DEM based on aerial photographs covers only the UI coastal area, reaching at most 40 $\mathrm{km}$ inland and covering primarily the UI-3 area. Simulated ice is 20–100 $\mathrm{\%}$ thicker around UI-1 and UI-2 than the 1985 observations and 10 $\mathrm{\%}$ thinner on average along UI-3.

NASA Operation IceBridge provides ice surface elevation along UI-1 (2009–2012) and UI-3 (1994, 1999, 2002, 2009, 2010, 2012). A mean value comparison along the UI-3 flow line illustrates that the simulated ice thickness is on average 10 $\mathrm{\%}$ less than observations (Fig. ). The same comparison on UI-1 shows simulated ice thickness being 104 $\mathrm{\%}$ of observations close to the UI-1 terminus and 93 $\mathrm{\%}$ of observed thickness 5 to 10 $\mathrm{km}$ upstream the 2012 terminus. Observed IceBridge and simulated surface elevation along flow lines 5 $\mathrm{km}$ downstream and 20 $\mathrm{km}$ upstream of the 2012 terminus have high correlation with $R^{\mathrm{2}}$ values of 0.80 for UI-1 and 0.95 for UI-3.

ISSM${}_{\text{PT}}$ simulates the major thinning trends as described by and between 1985 and 2010 on UI-1 and UI-3, though not on UI-2. Note that the observed thinning south of UI-3 between 1985 and 1991 is not reproduced in ISSM${}_{\text{PT}}$.

Ice surface velocity

Observed vs. simulated ice surface velocity along the central flow lines of UI-1, UI-2 and UI-3. Stars mark mean velocity between 0 and 5 $\mathrm{km}$ from the 2012 terminus, dots refer to mean values 5–10 $\mathrm{km}$ upstream. Winter velocity maps for 1991/92, 1992/93, 1994/95, 2002/03, 2005/06 and 2008/09 are produced from data available from http://esa-icesheets-greenland-cci.org/ and described in . Winter velocity maps from 2000/01, 2007/08 and 2009/10 are given by MEaSUREs .

[Figure omitted. See PDF]

By the end of the ISSM${}_{\text{PT}}$ simulation, ice flow velocity has doubled at UI-1 and UI-2 and increased by 55 $\mathrm{\%}$ at UI-3 compared to 1849 (relating the 1849 and 2012 velocity along each flow line between 1 and 30 $\mathrm{km}$ upstream the 2012 terminus position). The simulated ice surface velocity evolution in plan view over the study period can be viewed in movie02 (Supplement). Short-term accelerations coincide with the induced ice mass change due to the prescribed terminus change (see movie01, Supplement). The simulation reproduces seasonal and annual velocity variations due to the SMB forcing in the model. Small (20 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$) annual velocity fluctuations are forced by seasonal SMB fluctuations. Each retreat from the prescribed terminus change is followed by acceleration between 1 and 70 and 5–30 $\mathrm{\%}$ surface lowering, lasting 0.5 to 6 months.

Simulated 2009 ice surface velocity is within $\pm \mathrm{20}$ $\mathrm{\%}$ of observations from , except in the shear margins, where simulated velocities are up to 250 $\mathrm{\%}$ higher than observations. Winter velocity maps between 1991 and 2010 (Table ) are used to evaluate recent changes in simulated velocity. Observed and simulated winter ice surface velocity averaged between 0 and 5 $\mathrm{km}$ and 5 to 10 $\mathrm{km}$ upstream of 2012 terminus (Fig. ) have $R^{\mathrm{2}}$ values of 0.90 on UI-1, 0.88 on UI-2 and 0.92 on UI-3. Observations show 20 $\mathrm{\%}$ velocity increase on UI-1 from 2007/08 to 2008/09, however, this is not captured in ISSM${}_{\text{PT}}$.

Discussion

The comparison of ISSM${}_{\text{PT}}$ and ISSM${}_{\text{control}}$ shows that the ice surface velocity and thickness are significantly affected by the prescribed marginal changes. After each prescribed terminus change, ISSM${}_{\text{PT}}$ simulates short (0.5 to 6 months) periods of faster flow (1–70 $\mathrm{\%}$ acceleration), and the surface elevation lowers up to 30 $\mathrm{\%}$ at the new terminus in response to the ice flow acceleration. These are dynamic readjustments to the instantly reduced terminal flow resistance from the prescribed retreat, which is induced in discrete time steps.

While ISSM${}_{\text{PT}}$ produces maximum velocity and surface elevation changes of 275 and 84 $\mathrm{\%}$ respectively over the simulation period, ISSM${}_{\text{control}}$ simulates minor changes (maximum $\pm \mathrm{7}$ $\mathrm{\%}$) in ice thickness and velocity, representing sole mass changes prescribed by $\mathrm{\Delta }$SMB. This highlights the importance of simulated terminus retreat in order to reproduce a UI glacial system evolution.

In 1985–2002/05, ISSM${}_{\text{PT}}$ simulates mass changes similar to observations . However, observed and simulated mass changes between 2000 and 2008 differ from each other, showing up to 50 $\mathrm{\%}$ more simulated than observed total mass loss (Table ). During this time, the largest area changes are prescribed at UI-1 with increasing retreat rates from 200–500 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ between 1985 and 2005 to 4 $\mathrm{km}{\mathrm{yr}}^{-\mathrm{1}}$ in 2006–2009 along the UI-1 centre flow line. Between 2006 and 2007, UI-1 splits into three calving fronts, followed by the disintegration of UI-1's floating ice tongue in 2008 . However, the simulated UI-1 terminus is grounded between 1990 and 2012, only starting to float above an overdeepening in the bathymetry in 2007 in order to be grounded again after the prescribed retreat in 2008. The misfit between observed and simulated mass change could be justified by the absent floating tongue in the simulation. The UI-1 bathymetry is deeper than 500 $\mathrm{m}$ below sea level. Therefore, prescribed retreat leads to higher ice mass loss retreating over the same distance than an observed disintegration of a floating tongue. Moreover, a floating tongue has potential to stabilize the glacier more, decreasing the glacier velocity and dynamical discharge .

Recent studies suggest dividing mass balance into atmospheric and dynamically driven processes . Our simulation indicates three distinct MB periods when considering $\mathrm{\Delta }$SMB and $\mathrm{\Delta }$DIL. From the simulation start in 1849 to 1932, the total UI MB is the same for ISSM${}_{\text{control}}$ and ISSM${}_{\text{PT}}$, only diverging five times by one to four Gt ${\mathrm{yr}}^{-\mathrm{1}}$ when prescribed retreat is enforced. The increasing $\mathrm{\Delta }$SMB trend leads to a positive MB and thus mass gain. ISSM${}_{\text{PT}}$ velocity starts to differ from ISSM${}_{\text{control}}$ following the first prescribed retreat in 1862, showing a short ($<\mathrm{1}$ month) acceleration. The simulation indicates stable glacier behaviour without dynamically caused acceleration or thinning.

From 1925 onwards, $\mathrm{\Delta }$SMB reveals a negative trend, initiating the negative MB trend that lasts until the simulation end in 2012. Between 1931 and 1992, in two instances (1931–1960 and 1960–1992), 5–7 year periods of sustained less-positive SMB are followed by approximately 20 year long periods of elevated $\mathrm{\Delta }$DIL.

Within the 60 years of simulation, 31 terminus changes are prescribed, each removing 0.4–5 $\mathrm{Gt}$ of ice, which is as much as each of the five terminus changes during the preceding 82 years (Fig. ). The simulated mass loss in this period is therefore highly controlled by the prescribed retreat. $\mathrm{\Delta }$DIL consists of the removed ice mass at a prescribed retreat and of changes in ice mass flux caused by the acceleration of the glacier. We simulate two increased $\mathrm{\Delta }$DIL periods preceded by low $\mathrm{\Delta }$SMB as the result of observed terminus retreat. Induced by the prescribed terminus change in 1960 and 1966, a new period with increased $\mathrm{\Delta }$DIL lasts until 1992.

From 1999 onwards, $\mathrm{\Delta }$DIL and $\mathrm{\Delta }$SMB are roughly equivalent in contribution to the elevated negative MB. The simulation computes elevated dynamic ice loss due to 5.5 $\mathrm{km}$ terminus retreat on UI-1 within 12 years. We can not resolve, whether the increased retreat of UI-1 is due to (1) the change in $\mathrm{\Delta }$SMB from positive to negative values ($+$7 to $-$7 $\mathrm{Gt}$) or whether the glacier itself reaches an unstable position. However, as a result, the retreat causes increased simulated $\mathrm{\Delta }$DIL adding up to the same amount as the increased negative $\mathrm{\Delta }$SMB; UI-2 shows similar behaviour. The result is a negative MB, twice as negative as in any year before. In contrast, UI-3 is nearly stable, retreating $\sim \mathrm{0.5}$ $\mathrm{km}$ between 1999 and 2012 and even advances in some years. It cannot be determined, whether UI-1 and UI-2 also will reach a stable position soon or whether they will continue to retreat and accelerate.

Although we primarily discuss prescribed ice margin retreat, it is worth mentioning that our method also includes advancing observed terminus position changes at UI-1 and UI-2 in summer 2012 and at UI-3 in the summers 2001, 2003 and 2007. When ice margin advance is prescribed, ice mass is advected downstream and extrapolated over the regions that are newly activated.

MEaSUREs data indicate a 20 $\mathrm{\%}$ speed-up on UI-1 from 2007–2008 to 2008–2009, when a large floating ice tongue breaks off . Yet, the observed acceleration is not captured by the simulation and may be related to unresolved loss of buttressing in the simulation.

Simulated frontal changes occur on annual time scales, marked by observation years, and thereby happen less often than in nature. In addition, the magnitude of the removed ice mass is defined by the mesh resolution, set between 300 and 800 $\mathrm{m}$ on the area of frontal retreat. Hence, simulated frontal changes appear more abrupt than in nature and do rather capture changes in glacier velocity and thickness on decadal time scales than a seasonal pattern. The simulation is likely to overestimate the velocity and thickness changes in response to the larger decrease in frontal backstress.

Moreover, the timing of the prescribed terminus changes of ISSM${}_{\text{PT}}$ is not well constrained, given by observations with gaps of up to 13 years after 1900. In ISSM${}_{\text{PT}}$, we include 20 additional calving fronts, to minimize the ice mass removed at each terminus change. A simulation excluding the 20 interpolated terminus positions does not affect the overall simulation results (see Supplement). The total ice mass change shows higher peaks, but results in similar cumulative mass changes. Simulated 2012 ice velocity and thickness of ISSM${}_{\text{PT}}$ and the sensitivity simulation agree as well with $\pm \mathrm{5}$ $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ and $\pm \mathrm{3}$ $\mathrm{m}$ respectively.

However, given more frequent observations could minimize the timing uncertainties and with multiple observed terminus front positions throughout a year, this approach could capture seasonal glacier changes.

Conclusions

Our study shows that prescribing glacier front positions and surface mass balance are necessary to realistically simulate the multi-decadal evolution of ice velocity and thickness at Upernavik Isstrøm. Our simulation suggests that dynamic response caused by prescribed terminus position change is responsible for 70 $\mathrm{\%}$ of the total simulated mass change. Thus, moving terminus positions play an important role for UI's acceleration and thinning. The simulation with prescribed terminus changes reproduces distinct mass loss periods of dynamically driven ice mass loss and extends the periods discussed in and from 1985 to 1932.

Prescribed terminus position change avoids calving and melt rate estimations and reduces simulated retreat uncertainty. Yet, our approach requires knowledge about terminus positions and thus cannot be applied in future projections. However, the simulation results show the importance of calving in order to produce velocity and thickness change of tidewater glaciers. Better physically based calving laws are needed to understand and predict future glacier behaviour and glacier contribution to sea level rise. With an increasing amount of collected observed terminus front positions, the method discussed here will become a progressively more useful tool for evaluating calving laws or determining calving law parameters for hind-cast simulations before they are applied to future simulations. Short-term simulations with prescribed terminus position changes can determine what observations are needed to evaluate and construct new calving laws, by establishing if seasonal terminus position variations are necessary to capture long-term glacier behaviour. Future work could include comparisons with simulations using physically based calving laws e.g. as well as the application of our method to other tidewater glaciers.

Data from simulation results are available upon request from the authors.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-12-1511-2018-supplement.

KH, NJS, EYL, JEB, KHK, KKK and SHL designed the study and setup the model. KH performed the study and data comparison and led the writing of the manuscript, in which she received input and feedback from all authors. MM created the bed geometry from bathymetry data, supported with data from ER and TKD. KKK provided trimline data and observed terminus positions. AMS processed winter ice surface velocity maps from ESA, CCI Greenland.

The authors declare that they have no conflict of interest.

Acknowledgements

The manuscript was improved substantially by the editors and reviewers remarks, and the authors would like to thank Olivier Gagliardini, Jeremy Bassis and the two anonymous reviewers for their constructive comments. This study is part of the project “Multi-millennial ice volume changes of the Greenland Ice Sheet” funded by the Geocenter Danmark. Kristian K. Kjeldsen was supported by the Danish Council for Independent Research (DFF-4090-00151). We thank Brian Vinter and his team at Niels Bohr Institute, Copenhagen University, for generously supplying high performance computing resources. We wish to thank Camilla Snowman Andresen, Geological Survey of Greenland and Denmark, for providing bathymetry measurements. Observed termini between 1999 and 2012 were digitized by Trine S. Jensen and Karina Hansen, Geological Survey of Greenland and Denmark. We acknowledge and thank the Ice Sheet System Model group for producing and making available their model. We also acknowledge the use of the DEMs from GIMP, ArcticDEM and Niels Korsgaard and the velocity data provided by ESA (CCI Greenland) and NASA. Edited by: Olivier Gagliardini Reviewed by: Jeremy Bassis and two anonymous referees