RESEARCH ARTICLE

How long will an Arctic mountain glacier survive? A case study of Austre Lovénbreen, Svalbard

Zemin Wang, Guobiao Lin & Songtao Ai

Chinese Antarctic Center of Surveying and Mapping, Wuhan University, Wuhan, China

Abstract

To study Arctic valley glacier responses to global climate change, the Elmer/Ice ice-flow model was used to investigate long-term changes in Austre Lovénbreen, a typical polythermal glacier in Svalbard. Evolution and features, including volume, area, ice thickness, runoff and time and mode of glacier disappearance, were projected. Firstly, steady-state simulations were performed to determine the best parameters for the ice-flow model, which were then used to simulate glacial dynamics. Based on the 21st-century Arctic warming trend in the fifth assessment report published by the Intergovernmental Panel on Climate Change, the evolution of the glacier was simulated under three hypothetical climatic scenarios: pessimistic, high-probability and optimistic. The results predicted that the glacier will retreat until disappearance under all three scenarios, and its disappearance time will likely be approximately 111 years (by 2120). Under all scenarios, glacier volume and area reductions will be slow at first, then fast and finally slow again at the end. In particular, glacial runoff will increase markedly until 2070 in the high-probability scenario, and the peak runoff will be double the current value.

Keywords
Climate scenarios; Elmer/Ice; basal sliding; peak runoff; glacier disappearance

Abbreviations
AL: Austre Lovénbreen
DEM: digital elevation model
ELA: equilibrium-line altitude
GPS: global positioning system
SMB: surface mass balance
w.e.: water equivalent

Introduction

Glaciers, which result from specific climatic conditions, are very sensitive to climate change, and many of them have retreated to varying degrees in a global warming context (e.g., Huss et al. 2010; Gardner et al. 2011; Małecki 2013; Barnard et al. 2014; Zhang et al. 2014; Pérez et al. 2018). Small mountain glaciers and icecaps have significant roles in climate and sea-level changes at decadal or century scales because of their short response times to climate change (Oerlemans & Fortuin 1992). Small glaciers also have impacts on river discharge and energy exchange within basins (Verbunt et al. 2003). Svalbard (10°–35°E, 74°–81°N) is located in close proximity to the warm West Spitsbergen Current (Matul et al. 2018), and the mass balance of glaciers in the archipelago is very sensitive to fluctuations in this current and the corresponding climate changes (Xu et al. 2007). At present, Svalbard is an international hotspot for glacier research. This archipelago has a total area of approximately 59 250 km², and more than 57% of this area is covered by more than 2100 glaciers (Nuth et al. 2013). Most glaciers in Svalbard are polythermal (Jania et al. 1996), and AL is a typical polythermal valley glacier based on ice temperature measurements in recent years (Sun et al. 2016). AL is located in north-western of Svalbard (Fig. 1). AL is 6.2 km from Yellow River Station, China’s first Arctic research station, built in July 2004 (Ai, Wang, E, Pang et al. 2012). This research station is the base from which glaciologists have continuously conducted field observations of AL since 2005 (Ai et al. 2006).

Fig. 1  Geographic location of the glacier AL, Svalbard.

The glacierized area of Svalbard has decreased by 7% over the past 30 years or so, with an average decrease rate of –80 km² a⁻¹ (Nuth et al. 2013). According to Lang et al. (2016), the total SMB over all of Svalbard was negative (–1.6 Gt a⁻¹) from 1979 to 2013. Möller & Kohler (2018) researched the glacier mass balance for all of Svalbard from 1900 to 2010 and found that the average annual mass balance was –2 mm w.e. a^–1. Based on the time series analysis of high-resolution DEMs, small glaciers in Svalbard have experienced accelerated thinning from 1990 to 2005 (James et al. 2012). Since 1990, most of the small glaciers in central Spitsbergen, the largest island in Svalbard, have been thinning continually at a rate higher than the regional average (Małecki 2016). Considering the impact of superimposed ice, Zwinger & Moore (2009) performed a 53-year simulation of Midtre Lovénbreen, a nearby glacier on the west of AL, and predicted that the glacier tongue would retreat by approximately 500 m from 1977 to 2030. Through comparative analysis of glacier terrain, another study showed that the Pedersenbreen glacier, which is near AL, had experienced a marked retreat in 1936, 1990 and 2009 (Ai et al. 2013). Similar to nearby studied glaciers in the area, AL experienced a pronounced retreat from 1948 to 2013: the total retreat was 1064 ± 20 m in length and 1.82 ± 0.28 km² in area (Marlin et al. 2017). GPS observations indicated that the maximum ice-flow velocity along the glacier centre line of AL was 3.8 m a^–1 from 2005 to 2010 and that the mean velocity of all observed stakes was 2.4 m a⁻¹ during this period (Ai, Wang, E & Yan 2012).

AL has been the subject of numerous studies that have examined the glacier’s mass balance (Xu et al. 2010; Marlin et al. 2017), topography (Ai, Wang, E, Pang et al. 2012; Saintenoy et al. 2013), temperature changes (Sun et al. 2016) and motion (Bernard et al. 2011; Ai, Wang, E & Yan 2012). However, little research has been done on the future evolution of AL, and no scholar has used numerical simulation methods to predict the evolution and features of the glacier in the future.

In this study, we used an ice-flow model—Elmer/Ice—to simulate the future evolution of AL, including the geometric and physical responses of AL to climate changes during its decay under certain climatic scenario to address such questions as: How will AL respond to the future climate warming? Will AL become extinct? If so, when will it eventually disappear?

Materials and methods

Data sources

The data for this study included a surface DEM, a bedrock DEM and the surface ice-flow velocity of AL, all collected in 2009 (Ai, Wang, E, Pang et al. 2012). Field measurements were performed on the glacier using real-time kinematic and ground-penetrating radar (Ai et al. 2014) by a snowmobile carrying measuring instruments. With a Leica GS10 real-time kinematic set, kinematic GPS surveying provided the entire surface topography of AL at centimetre-level accuracy. The surface DEM was obtained from real-time kinematic measurements by natural neighbour interpolation, and the bedrock DEM was provided by subtracting the ice thickness obtained by ground-penetrating radar from the surface topography. The spatial resolutions of both DEMs were 10 m. The surface ice-flow velocity was acquired from repeated measurements of stakes by double-frequency GPS, which has a positioning accuracy better than 1 cm. The surface DEM, bedrock DEM and locations of stakes are shown in Fig. 2. The maximum ice thickness was 155 m in 2009, and the maximum ice velocity (3.8 m a^–1) measured by GPS along the glacier centre line occurred at point C2 (Ai, Wang, E, Pang et al. 2012; Fig. 2a).

Fig. 2  (a) Surface, and (b) bedrock DEMs of the glacier AL. The contour interval for both DEMs is 25 m. The spatial distribution of stakes is also shown in (a).

There were three ice-temperature observation points located on the glacier centreline, and the ice temperature was observed every May and September during 2009–2011 (Sun et al. 2016). Although the thermal structure of AL may be as complex as that of Midtre Lovénbreen (Zwinger & Moore 2009), we only use the in situ ice temperature observations at these three points to approximately estimate the temperature profile of AL. According to the ice temperature observations by Sun et al. (2016), the average ice temperature was approximately –3.0 °C at a 14-m depth in 2009 and increased linearly with increasing depth, reaching a maximum value of 0 °C (consider 0 °C as the pressure melting point) at 100-m depth. Therefore, for the purpose of the ice-flow model, the ice temperature (in Celsius), T, was approximated as follows:

where d is the depth and T₀ is the average ice temperature at lower boundary of the active layer, which is –3.0 °C at a depth of 14 m in 2009. This relationship was assumed to be constant over time in the model.

In the IPCC’s fifth assessment report (Collins et al. 2013), there were four representative concentration pathways for greenhouse gases, referred to as RCP2.6, RCP4.5, RCP6.0 and RCP8.5. The Arctic region is expected to experience varying warming from the end of the 20th century to the end of the 21st century; in the report, warming is projected under different scenarios: under RCP2.6, RCP4.5, RCP6.0 and RCP8.5, the temperature of the region will increase by 2.2, 4.2, 5.2 and 8.3 °C, respectively (Collins et al. 2013). Compared with the upper (RCP8.5) and the lower (RCP2.6) end of emission pathways, the intermediate scenarios RCP4.5 and RCP6.0, representing an average increase of 4.7 °C, are more likely to occur. Accordingly, during this 100-year period, the Arctic is expected to warm by 8.3 °C under the most pessimistic, high-emission scenario (RCP8.5), with an average warming rate of 0.083 °C a⁻¹; most likely, it will warm by 4.7 °C, with an average warming rate of 0.047 °C a⁻¹ (the synthesis of RCP6.0 and RCP4.5). Under the most optimistic scenario (RCP2.6), the Arctic will warm by 2.2 °C, with an average warming rate of 0.022 °C a⁻¹. In consideration of the report, three hypothetical future climate scenarios are considered in this study: (a) pessimistic scenario in which the Arctic will continue to warm rapidly at a rate of 0.083 °C a⁻¹, (b) a high-probability scenario in which the Arctic will continue to warm at a rate of 0.047 °C a⁻¹ and (c) an optimistic scenario in which the Arctic will continue to warm slowly at a rate of 0.022 °C a⁻¹.

Mass balance model

Glacier mass balance, as an important driver of glacier change, is closely related to the climate (Pelto 2018). Mass balance usually refers to SMB and may be approximated using two key parameters: the ELA and the gradient. The ELA is closely related to air temperature (Zhao et al. 2014) and the gradient is the rate of change in SMB with elevation. Generally, the SMB increases with increasing elevation, which can be extracted from previous works (Kohler et al. 2007; Välisuo et al. 2017). Xu et al. (2010) found a strong linear correlation between the SMB of AL and surface elevation during the balance year 2005–06. Zhang (2011) found a similar relationship between the SMB and surface elevation based on measurements over 2005–2010, reporting a mean gradient of 1/200 m^–1 m a^–1 ice equivalent. Hence, the annual SMB at different values of surface elevation, Z_S, can be calculated as follows:

The unit of SMB is m a⁻¹ ice equivalent. This relationship is assumed to remain unchanged over the entire period of glacier simulation.

The ELA of AL in 2009 was ca. 423 m according to the data published by the World Glacier Monitoring Service (2017). The sensitivity of ELA to change in temperature, α, is 90 m/°C, that is, the ELA increases (decreases) 90 m for every 1 °C increase (decrease) in air temperature (Fleming et al. 1997). Since most of the change in ELA is attributed to variations in temperature on the peninsula of Brøggerhalvøya, north-western Svalbard (Kohler et al. 2002), we ignore the effect of precipitation on ELA, and consider only the effect of temperature on ELA in this study. ELA in the ith year starting from 2009 was computed as

where ΔT_i is the net temperature change of the ith year relative to 2009.

Glacier runoff

Glacier runoff is an important hydrological process mainly determined by precipitation and glacial meltwater (Xu et al. 2017; Huss & Hock 2018). Currently, AL has almost no net accumulation; the yearly glacier runoff can also be defined as the sum of annual precipitation and ice loss. The annual precipitation in the AL basin can be calculated by multiplying the total precipitation (427 mm a⁻¹, average for 1981–2010 from the nearby meteorological station) (Førland et al. 2011) by the basin area (10 km²) (Bernard et al. 2013). Considering that annual precipitation will increase in the future by 5.4% per decade according to the linear precipitation trend for 1975–2011 (Førland et al. 2011) and that the reduction in glacier volume occurs entirely through transformation into meltwater, the future basin runoff is estimated as the sum of annual precipitation and meltwater. To coincide with our yearly simulation, the future precipitation trend is simply converted to 0.54% a^–1.

Ice-flow model

The glacier is taken as an incompressible fluid with a constant density, and its flow is governed by the Stokes equations (Zwinger & Moore 2009):

where u is the ice velocity, τ is the deviatoric stress, p is the ice pressure, g is the gravitational acceleration and ρ is the ice density.

Assuming the glacier to be an isotropic material, its constitutive equation satisfies Glen’s law (Réveillet et al. 2015):

where is the strain rate, and the effective viscosity η is defined as:

where is the square of the second invariant of the strain rate, E is the Glen enhancement factor (measured from 0.06 to 10, usually E = 1) and n is the Glen exponent (usually n = 3). A is a rheological parameter depending on the ice temperature (in Kelvin) relative to the pressure melting point, T′(Greve & Blatter 2009):

where A₀ is the rate factor, Q is the creep activation energy and R is the gas constant.

The glacier may be subject to basal sliding due to the presence of sub-glacial meltwater. The basal boundary condition (Réveillet et al. 2015) can be expressed as follows:

where τ_β is the basal shear stress, β is the basal friction parameter and u_b is the basal tangential velocity.

The surface elevation Z_S at any point (x, y) changes with time t following the kinematic boundary condition (Zhao et al. 2014):

where u_x, u_y and u_z are the three components of ice-flow velocity in the three directions x, y and z, respectively.

The parameters of the ice-flow model in this study are shown in Table 1.

Table 1 The parameters of the ice-flow model.

Symbol	Description	Value and unit
ρ	Ice density	910 kg m−3
g	Gravitational acceleration	9.81 kg s−2
n	Glen exponent	3
A0	Rate factor
	When T ≤ –10 °C	2.89 × 10−13 s−1 Pa−3
	when T > –10 °C	2.43 × 10−2 s−1 Pa−3
Q	Creep activation energy
	when T ≤ –10 °C	60 kJ mol–1
	when T > –10 °C	115 kJ mol–1
R	Gas constant	8.31 J kg–1 K–1

Simulation process

The simulation process consisted of two phases. Initially, steady-state simulation (diagnostic simulation) constrained by GPS data was performed; then, transient simulation (prognostic simulation) was performed using the steady-state simulation result as initial condition. Both phases of simulation are described in the following sections.

Ice velocity simulation setup. The basal friction parameter β and the Glen enhancement factor E are two important parameters in the ice-flow model. In this study, the model parameters β and E were obtained by comparison between the measured ice-flow velocities from GPS and simulated values from multiple steady-state simulations.

Firstly, a glacier finite element grid was constructed. The boundary coordinates of AL in 2009 were input into an open-source finite element generator (Gmsh) (Geuzaine et al. 2009) using the Northern Hemisphere 33° projection of the Universal Transverse Mercator projection (UTM33N), and the glacier’s two-dimensional plane was divided into a network of irregular triangles with side lengths of about 50 m (Fig. 3a). Next, the whole glacier was evenly extruded into 15 ice layers (Fig. 3b) with the surface and bedrock DEMs as the upper and lower boundaries, respectively. Thereafter, steady-state simulation was carried out to obtain glacial physical parameters (e.g., velocity field, stress field and pressure field). The observation point C2 (Fig. 2a), located on the glacier centreline, was the point of maximum velocity from early GPS observations. In this study, the measured velocity at C2 was selected as a control quantity to estimate the basal friction parameter and the Glen enhancement factor, that is, when β took a certain value, E was assigned an appropriate value; the combination of β and E was used to make simulated ice velocities match well with GPS-measured velocities at point C2. Based on the above premise, by comparing simulated ice velocities of each combination of β and E with measured velocities along the glacier centreline, the best combination with the smallest overall deviation was selected as the optimal solution for model parameters.

Fig. 3  (a) Finite element mesh, and (b) vertical stratification of the glacier AL.

The main purpose of steady-state simulation is to obtain the values of β and E. Although there is meltwater at the base of polythermal glaciers, the amount of meltwater varies with time and directly affects the basal conditions (Werder et al. 2013). Since these are unknown for AL, in this study, three basal assumptions are considered: no sliding, full sliding (the entire glacier base is sliding with a certain velocity) and partial sliding.

Figure 4 and Table 2 show that the simulation results for full sliding at the base are closer to the measured values than those for no sliding at the base. This finding indicates that sliding phenomena are expected under certain areas of AL. However, the simulation results match well with the measured values at all stakes except A1, when there is full sliding at the glacier base (Fig. 4). The measured ice-flow velocity is far less than the simulated velocity at stake A1, where the ice thickness is ca. 20 m, suggesting that the glacier base is frozen around the margins of the bed. Therefore, the steady-state simulation was performed again, accounting for this effect in places where ice thickness is less than 20 m.

Fig. 4  Simulation results for the states of no sliding and full sliding at the glacier base. The black line denotes the measured velocities. The grey line denotes the simulation results for no sliding, and the other colours denote the results for full sliding.

Table 2 Mean value and standard deviation of the difference between measured and simulated velocities (m a^–1) for the states of no sliding and full sliding at the glacier base.

	Full sliding							No sliding
βa	0.04	0.05	0.06	0.07	0.08	0.09	0.10	∞b
Ec	0.667	0.814	0.914	0.987	1.042	1.085	1.120	1.430
Mean value of difference	0.286	0.160	0.074	0.013	–0.034	–0.071	–0.099	–0.375
Standard deviation of difference	0.267	0.216	0.197	0.195	0.199	0.206	0.213	0.334
aThe basal friction parameter (unit: MPa a m-1 ). bThe glacier has no basal sliding. cThe Glen enhancement factor.

In the absence of sliding, the simulated velocities match well with the measured values at stake A1 (Fig. 5), suggesting that the glacier has a sufficiently cold ice-bed interface to prevent basal sliding when the ice thickness is less than 20 m. Standard deviation of the difference decreases at first and then increases when β is between 0.01 and 10 000, and the absolute value of the mean difference has a similar trend (except β = 0.01) (Table 3). When β is 0.05 and E is 0.865, the mean value is 0.010 m a^–1, the standard deviation is 0.175 m a^–1, and their absolute values reach minima as shown in Table 3; thus, β of 0.05 and E of 0.865 are the best parameter estimates for the model.

Fig. 5  Simulation results for partial sliding at the glacier base.

Table 3 Mean values and standard deviation of the difference between measured and simulated velocities (m a^–1) for partial sliding at the glacier base.

βa	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08
Eb	0.037	0.213	0.530	0.732	0.865	0.956	1.023	1.073
Mean value of difference	0.246	0.311	0.216	0.092	0.010	–0.050	–0.092	–0.125
Standard deviation of difference	0.814	0.469	0.267	0.195	0.175	0.176	0.185	0.196
βa	0.09	0.10	0.50	1.00	10	100	1000	10000
Eb	1.112	1.143	1.374	1.402	1.427	1.430	1.430	1.430
Mean value of difference	–0.152	–0.174	–0.332	–0.353	–0.373	–0.374	–0.375	–0.375
Standard deviation of difference	0.207	0.216	0.308	0.321	0.332	0.334	0.334	0.334
aThe basal friction parameter (unit: MPa a m-1 ). bThe Glen enhancement factor.

Despite the fact that there certainly are spatial and temporal variations in the sliding coefficient, they are difficult to constrain with the given data. Therefore, these variations were ignored in this study, and an overall optimal basal friction coefficient was selected and assumed to be invariant with time, which simplified the simulation process and improved the calculation efficiency. In this study, the model parameters of the Elmer/Ice model are as follows: E = 0.865 and β = ∞ (when the ice thickness is less than 20 m) or β = 0.05 (when the ice thickness is greater than or equal to 20 m). Figure 6 shows the error graph using these model parameters.

Fig. 6  Error graph using the best estimates of model parameters. Black points denote horizontal velocities at the stakes depicted in Fig. 2, and blue vertical lines represent the errors between the simulated and measured velocities.

Model uncertainties. Owing to assumptions and simplifications in our model, many factors have impacts on the modelling results, such as the constant SMB gradient, inter-annual mass balance variations, debris cover distribution, insolation of glacier surface, superimposed ice formation and internal accumulation. Estimating the impact of each and every factor is extremely difficult, so we used a backtracking simulation to verify modelling results with measured values.

Using the historical meteorological records (Førland et al. 2011), together with the data of past glacier terrain (Marlin et al. 2017), a backtracking simulation from 1962 to 2009 was performed to verify the feasibility of the model. In this simulation, the model parameters (β = 0.05, E = 0.865) and the bedrock terrain were assumed to be unchanged in the past 50 years. We performed a steady-state simulation under the climatic condition of 1962, and then used this steady-state simulation result as initial condition to perform a transient simulation from 1962 to 2009. The biases of glacier features (area, volume and ice thickness) between the measurements and simulation in 2009 were calculated. These biases can be yearly averaged as the unified uncertainty rates in percentage for modelling, which can be applied to the future simulation results.

With the backtracking simulation from 1962 to 2009, the biases of area, volume and ice thickness in 2009 were 3, 14 and 11.5%, respectively. Therefore, the uncertainties for the future simulations were set to 0.06% a^–1 for area, 0.3% a^–1 for volume and 0.25% a^–1 for ice thickness.

Using the optimized model parameters from the steady-state simulation, transient simulations of AL were carried out to project the glacier’s future evolution from 2010. The entire workflow of the glacier modelling is shown in Fig. 7.

Fig. 7  Flow chart for glacier modelling.

Results and discussion

Using the optimized model parameters from the steady-state simulation and taking the steady-state simulated results as the initial values, the future evolution of AL was simulated under the pessimistic, high-probability and optimistic scenarios from 2010 to glacier disappearance. The simulation results are as follows.

Glacier ELA and mass change

Although the glacier ELA varies from year to year, predicting its future variation is difficult. Therefore, a simplified linear trend was applied for the future ELA evolution according to Eqn. 3 such that glacier ELA gradually increases as in Fig. 8a. With increasing ELA, the glacial accumulation area ratio is quickly reduced to zero, as shown in Fig. 8a, which means that the glacier is under full melting state from the terminus to the head area.

Fig. 8  (a) Glacier ELA and accumulation area ratio (AAR) evolution; (b) stake elevation change under the high-probability scenario.

Regarding the mass changes in specific locations, the stakes along the centreline (A2, B2, C2, D3, E2 and F as shown in Fig. 2a) were chosen to evaluate their time series surface changes. Each stake has a decay curve of elevation, as shown in Fig. 8b. Generally, the higher the elevation of the stake, the longer the ice at stake survives. However, the ice thickness at the stake also plays an important role in the vanishing of ice at a stake position. For example, under the high-probability scenario, stake E2 is at a higher elevation than stake D3 but will disappear about four years earlier than D3 because of its thinner ice cover.

Glacier volume and area

Under all scenarios, the glacier is in a continuously decreasing state until disappearance (Fig. 9). Glacier volume and area decline slowly at first, then rapidly and finally slowly again at the end. Under the pessimistic, high-probability, and optimistic scenarios, the disappearance periods of AL are 84 years (by 2093), 111 years (by 2120), and 161 years (by 2170), respectively. During its decay (from 2009 to its demise, the same below), the mean volume loss rates are 0.0040, 0.0030 and 0.0021 km³ a⁻¹, respectively; and the mean area reduction rates are 0.055, 0.042 and 0.029 km² a⁻¹, respectively. Note that uncertainty in glacier volume is very large in long-term projections.

Fig. 9  Simulated glacier (a) volume, and (b) area as functions of time in the pessimistic, high-probability and optimistic scenarios. The dotted line in (a) denotes half the glacier volume in 2009, and the dotted line in (b) denotes half the glacier area in 2009.

Compared with the retreat of AL from 1962 to 2013 (mean volume loss rate: 0.0025 km³ a⁻¹, mean area reduction rate: 0.026 km² a⁻¹) (Marlin et al. 2017), the glacier will melt quickly in the future. The mean area reduction rate of AL is also similar to those of Midtre Lovénbreen (0.019 km² a^–1 over 1966–1990) (Moreau et al. 2008) and Austre Brøggerbreen (0.025 km² a⁻¹ over 1936–2010) (Tsuji et al. 2016). Outside the Arctic archipelago of Svalbard, Grosser Aletschgletscher, a glacier in the Swiss Alps, is expected to lose 90% of its present volume by the end of the 21st century (Jouvet et al. 2011); Gurenhekou glacier in the southern Tibetan Plateau would lose 35% of its volume by 2057 under continued steady warming (Zhao et al. 2014); the volume loss of Glaciar Zongo, a tropical glacier in the Andes, will be 69% by 2100 under the intermediate scenario RCP6.0 (Réveillet et al. 2015). These studies reveal rapid melting trends for many mountain glaciers worldwide in the 21st century, and many small mountain glaciers will probably disappear in one or two centuries.

The evolution of ice thickness under the three scenarios is illustrated in Fig. 10. The retreat rate of the glacier terminus is ca. 7.9 m a⁻¹ over 2010–2060 under the high-probability scenario. This retreat rate is of the same order of magnitude as that of the nearby glacier Midtre Lovénbreen (ca. 10 ma⁻¹ from 1977 to 2030) (Zwinger & Moore 2009). However, the retreat rate of AL is lower than that of Midtre Lovénbreen. This difference is probably because the snout of AL is more constrained and smaller than that of Midtre Lovénbreen for their respective modelling periods. To compare with a glacier in the southern Tibetan Plateau, the terminus retreat rate of Gurenhekou glacier will be ca. 9.1 m a⁻¹ during 2008–2057 (Zhao et al. 2014), which is comparable with those of small glaciers in Svalbard, for example, AL and Midtre Lovénbreen. When the breaking between mainstream and east tributary occurs, the mainstream will also be separated into two parts under the high-probability and pessimistic scenarios, whereas the mainstream will be still a whole under the optimistic scenario. By the late decay stage, only a little glacier ice remains; the remaining ice covers three regions (the middle part, the head of the mainstream and the east tributary) under high-probability scenario, whereas it covers only two regions (the head of the mainstream and the east tributary) under the optimistic scenario, and another two regions (the middle and head parts) under the pessimistic scenario.

Fig. 10  Simulated glacier area and thickness in different years: (i) the 50th year (2060), (ii) disappearance of west tributary, (iii) breaking between mainstream and east tributary and (iv) the late decay under the (a) optimistic, (b) high-probability and (c) pessimistic scenarios.

Glacier runoff

The annual precipitation, glacial meltwater and glacial runoff are estimated as shown in Fig. 11. Under the high-probability scenario, glacier runoff peaks in 2070 approximately (10.34 × 10⁶ m³ a⁻¹), whereas the runoff peaks under the pessimistic and optimistic scenarios occur in 2057 and 2096, respectively (Fig. 11c). Under the high-probability scenario, the yearly runoff increases to as much as two times the level of current year (i.e., 2009). Under the pessimistic scenario, the peak runoff will be even larger. Since the uncertainty in the future precipitation trend is unknown, the uncertainty in glacial runoff is difficult to estimate, so we do not present the uncertainty scope in Fig. 11.

Fig. 11  (a) Predicted precipitation in the glacier basin as a function of time. (b) Simulated glacial meltwater as a function of time in the pessimistic, high-probability and optimistic scenarios. (c) Annual glacial runoff as a function of time in the pessimistic, high-probability and optimistic scenarios. The black dashed line denotes the annual precipitation in the glacier basin from 2010 to 2170.

The latest research indicates that the total annual runoff from six major glaciers near Ny-Ålesund is about 0.6 × 10⁸ m³ a⁻¹ in 2010 (Pramanik et al. 2018). AL is among the six glaciers which in total cover a glaciated area of about 42 km². Note that the non-glaciated area (ca. 75 km²) is also included in this total runoff estimation. The area averaged runoff for the AL basin is ca. 5 × 10⁶ m³ a⁻¹, which is close to our result for AL in 2010 (5.5 × 10⁶ m³ a⁻¹).

The increase in runoff may have strong impacts on the landforms below the glacier. For example, the summer runoff lowers the level of the riverbed in the tundra area, and larger quantities of sediments are potentially transported into the bay, which will probably produce a new delta on the shore. The increased runoff also triggers potential floods on glacier foreland and raises the risks involved in field expeditions.

Conclusions

In this study, the Elmer/Ice ice-flow model was used to simulate the future evolution of AL during its decay. The modelling process consisted of steady-state simulation and transient simulations. Combining the measured velocities from GPS with steady-state simulation, the relevant parameters for the ice-flow model were obtained. Using these model parameters, from 2010, transient simulations of AL were performed to predict its future evolution under three climatic scenarios of pessimism, high-probability and optimism.

In all three scenarios, AL eventually disappears. Its disappearance time is likely to be 111 years (in the high-probability scenario by 2120); the fastest simulated disappearance time is 84 years (in the pessimistic scenario by 2093), and the slowest is 161 years (in the optimistic scenario by 2170). Under the three scenarios, glacier volume and area will decline slowly at first, then rapidly, and finally slowly again at the end. Under the high-probability scenario, the mean volume and area loss rates will be 0.0030 km³ a⁻¹ and 0.042 km² a⁻¹, respectively. Under the pessimistic, high-probability and optimistic scenarios, glacial runoff will peak in 2057, 2070 and 2096, respectively. The peak runoff will be double the current value under the high-probability scenario.

Most glaciers in Svalbard are polythermal glaciers, and AL is the representative of this type of glacier. Based on the modelling results for AL under three climatic scenarios, other Arctic valley glaciers with similar sizes, elevations, and the current and projected future conditions might be expected to gradually retreat to disappearance. This process is an inevitable glacier evolutionary trend under current climatic conditions. Notably, although this study presents the evolution of AL for three scenarios, many parameters are simplified in the ice-flow model. In fact, the response of glaciers to climate warming is a complex process. Therefore, continual corrections by more in situ data are needed to make glacier predictions more accurate in the future.

Acknowledgements

We greatly appreciate field expedition support from the Chinese Arctic and Antarctic Administration, State Oceanic Administration, People’s Republic of China.

Author contributions

Z.W. and S.A. conceived and designed the experiments and collected field data. G.L. processed the data and wrote the article; all authors contributed to the finalization of the article.

References

Ai S., E D., Yan M. & Ren J. 2006. Arctic glacier movement monitoring with GPS method on 2005. Chinese Journal of Polar Science 17, 61–68.

Ai S., Wang Z., E D., Holmen K., Tan Z., Zhou C. & Sun W. 2014. Topography, ice thickness and ice volume of the glacier Pedersenbreen in Svalbard, using GPR and GPS. Polar Research 33, article no. 18533, doi: 10.3402/polar.v33.18533.

Ai S., Wang Z., E D., Pang X., Zhou C., Yan M., Sun J. & Liu H. 2012. Topographic survey on the surface of glacier Austre Lovénbreen and Pedersenbreen in Svalbard based on GPS method. Chinese Journal of Polar Research 24, 53–59.

Ai S., Wang Z., E D. & Yan M. 2012. Surface movement research of Arctic glaciers using GPS method. Geomatics & Information Science of Wuhan University 37, 1337–1340.

Ai S., Wang Z., Tan Z., E D. & Yan M. 2013. Mass change study on Arctic glacier Pedersenbreen, during 1936-1990-2009. Chinese Science Bulletin 58, 3148–3154.

Barnard A., Wellner J.S. & Anderson J.B. 2014. Late Holocene climate change recorded in proxy records from a Bransfield Basin sediment core, Antarctic Peninsula. Polar Research 33, article no. 17236, doi: 10.3402/polar.v33.17236.

Bernard E., Friedt J.M., Tolle F., Griselin M., Laffly D. & Marlin C. 2011. Using ground based high resolution photography for seasonal snow and ice dynamics (Austre Lovénbreen, Svalbard, 79°N). Paper presented at the 15th Alpine Glaciology Meeting, 25 February, Munich, Germany.

Bernard É., Friedt J.M., Tolle F., Griselin M., Martin G., Laffly D. & Marlin C. 2013. Monitoring seasonal snow dynamics using ground based high resolution photography (AustreLovénbreen, Svalbard, 79°N). ISPRS Journal of Photogrammetry and Remote Sensing 75, 92–100, doi: 10.1016/j.isprsjprs.2012.11.001.

Collins M., Knutti R., Arblaster J., Dufresne J.L., Fichefet T., Friedlingstein P., Gao X., Gutowski W.J., Johns T., Krinner G., Shongwe M., Tebaldi C., Weaver A.J. & Wehner M. 2013. Long-term climate change: projections, commitments and irreversibility. In T.F. Stocker et al. (eds.): Climate change 2013. The physical science basis. Contribution of Working Group I to the fifth assessment report of the Intergovernmental Panel on Climate Change. Pp. 1055–1056. Cambridge: Cambridge University Press.

Fleming K.M., Dowdeswell J.A. & Oerlemans J. 1997. Modelling the mass balance of northwest Spitsbergen glaciers and responses to climate change. Annals of Glaciology 24, 203–210, doi: 10.3189/S0260305500012180.

Førland E.J., Benestad R., Hanssen-Bauer I., Haugen J.E. & Skaugen T.E. 2011. Temperature and precipitation development at Svalbard 1900–2100. Advances in Meteorology 17, 13–30, doi: 10.1155/2011/893790.

Gagliardini O., Zwinger T., Gillet-Chaulet F., Durand G., Favier L., Fleurian B.D., Greve R., Malinen M., Martín C., Råback P. & Ruokolainen J. 2013. Capabilities and performance of Elmer/Ice, a new generation ice-sheet model. Geoscientific Model Development 6, 1299–1318, doi: 10.5194/gmd-6-1299-2013.

Gardner A.S., Moholdt G., Wouters B., Wolken G.J., Burgess D.O., Sharp M.J., Cogley J.G., Braun C. & Labine C. 2011. Sharply increased mass loss from glaciers and ice caps in the Canadian Arctic Archipelago. Nature 473, 357, doi: 10.1038/nature10089.

Geuzaine C. & Remacle J. 2009. Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering 79, 1309-1331, doi: 10.1002/nme.2579.

Greve R. & Blatter H. 2009. Dynamics of ice sheets and glaciers. Berlin: Springer.

Huss M. & Hock R. 2018. Global-scale hydrological response to future glacier mass loss. Nature Climate Change 8, 135–140, doi: 10.1038/s41558-017-0049-x.

Huss M., Hock R., Bauder A. & Funk M. 2010. 100-year mass changes in the Swiss Alps linked to the Atlantic Multidecadal Oscillation. Geophysical Research Letters 37, L10501, doi: 10.1029/2010GL042616.

James T.D., Murray T., Barrand N.E., Sykes H.J., Fox A.J. & King M.A. 2012. Observations of enhanced thinning in the upper reaches of Svalbard glaciers. The Cryosphere 6, 1369–1381, doi: 10.5194/tc-6-1369-2012.

Jania J., Mochnacki D. & Gadek B. 1996. The thermal structure of Hansbreen, a tidewater glacier in southern Spitsbergen, Svalbard. Polar Research 15, 53–63, doi: 10.3402/polar.v15i1.6636.

Jouvet G., Huss M., Funk M. & Blatter H. 2011. Modelling the retreat of Grosser Aletschgletscher, Switzerland, in a changing climate. Journal of Glaciology 57, 1033–1045, doi: 10.3189/002214311798843359.

Kohler J., James T.D., Murray T., Nuth C., Brandt O., Barrand N.E., Aas H.F. & Luckman A. 2007. Acceleration in thinning rate on western Svalbard glaciers. Geophysical Research Letters 34, L18502, doi: 10.1029/2007GL030681.

Kohler J., Nordli Ø., Brandt O., Isaksson E., Pohjola V., Martma T. & Aas H.F. 2002. Svalbard temperature and precipitation, late 19th century to the present. Final report on ACIA-funded project. Oslo: Norwegian Polar Institute.

Lang C., Fettweis X. & Erpicum M. 2016. Stable climate and surface mass balance in Svalbard over 1979–2013 despite the Arctic warming. The Cryosphere 9, 83–101, doi: 10.5194/tc-9-83-2015.

Małecki J. 2013. Elevation and volume changes of seven Dickson Land glaciers, Svalbard, 1960–1990–2009. Polar Research 32, article no. 18400, doi: 10.3402/polar.v32i0.18400.

Małecki J. 2016. Accelerating retreat and high-elevation thinning of glaciers in central Spitsbergen. The Cryosphere 10, 1317–1329, doi: 10.5194/tc-10-1317-2016.

Marlin C., Tolle F., Griselin M., Bernard E., Saintenoy A., Quenet M. & Friedt J.M. 2017. Change in geometry of a High Arctic glacier from 1948 to 2013 (Austre Lovénbreen, Svalbard). Geografiska Annaler Series A 99, 115–138, doi: 10.1080/04353676.2017.1285203.

Matul A., Spielhagen R., Kazarina G., Kruglikova S., Dmitrenko O. & Mohan R. 2018. Warm-water events in the eastern Fram Strait during the last 2000 years as revealed by different microfossil group. Polar Research 37, article no. 1540243, doi: 10.1080/17518369.2018.1540243.

Möller M. & Kohler J. 2018. Differing climatic mass balance evolution across Svalbard glacier regions over 1900-2010. Frontiers in Earth Science 6, UNSP 128, doi: 10.3389/feart.2018.00128.

Moreau M., Mercier D., Laffly D. & Roussel E. 2008. Impacts of recent paraglacial dynamics on plant colonization: a case study on Midtre Lovénbreen foreland, Spitsbergen (79°N). Geomorphology 95, 48–60, doi: 10.1016/j.geomorph.2006.07.031.

Nuth C., Kohler J., König M., Deschwanden A.V., Hagen J.O.M., Kääb A., Moholdt G. & Pettersson R. 2013. Decadal changes from a multi-temporal glacier inventory of Svalbard. The Cryosphere 7, 1603–1621, doi: 10.5194/tc-7-1603-2013.

Oerlemans J. & Fortuin J.P. 1992. Sensitivity of glaciers and small ice caps to greenhouse warming. Science 258, 115–117, doi: 10.1126/science.258.5079.115.

Pelto M.S. 2018. How unusual was 2015 in the 1984–2015 period of the North Cascade Glacier annual mass balance? Water 10(5), article no. 543, doi: 10.3390/w10050543.

Pérez T., Mattar C. & Fuster R. 2018. Decrease in snow cover over the Aysén River catchment in Patagonia, Chile. Water 10(5), article no. 619, doi: 10.3390/w10050619.

Pramanik A., Van Pelt W., Kohler J. & Schuler T.V. 2018. Simulating climatic mass balance, seasonal snow development and associated freshwater runoff in the Kongsfjord basin, Svalbard (1980–2016). Journal of Glaciology 64, 943–956, doi: 10.1017/jog.2018.80.

Réveillet M., Rabatel A., Gillet-Chaulet F. & Soruco A. 2015. Simulations of changes to Glaciar Zongo, Bolivia (16°S), over the 21st century using a 3-D full-Stokes model and CMIP5 climate projections. Annals of Glaciology 56(70), 89–97, doi: 10.3189/2015AoG70A113.

Saintenoy A., Friedt J.M., Booth A.D., Tolle F., Bernard E., Laffly D., Marlin C. & Griselin M. 2013. Deriving ice thickness, glacier volume and bedrock morphology of AustreLovénbreen (Svalbard) using GPR. Near Surface Geophysics 11, 253–261, doi: 10.3997/1873-0604.2012040.

Sun W., Yan M., Ai S., Zhu G., Wang Z., Liu L., Xu Y. & Ren J. 2016. Ice temperature characteristics of the Austre Lovénbreen glacier in Ny-Ålesund, Arctic region. Geomatics and Information Science of Wuhan University 41, 79–85, doi: 10.13203/j.whugis20150302.

Tsuji M., Uetake J. & Tanabe Y. 2016. Changes in the fungal community of Austre Brøggerbreen deglaciation area, Ny-Ålesund, Svalbard, High Arctic. Mycoscience 57, 448–451, doi: 10.1016/j.myc.2016.07.006.

Välisuo I., Zwinger T. & Kohler J. 2017. Inverse solution of surface mass balance of Midtre Lovénbreen, Svalbard. Journal of Glaciology 63, 593–602, doi: 10.1017/jog.2017.26.

Verbunt M., Gurtz J., Jasper K., Lang H., Warmerdam P. & Zappa M. 2003. The hydrological role of snow and glaciers in alpine river basins and their distributed modeling. Journal of Hydrology 282, 36–55, doi: 10.1016/S0022-1694(03)00251-8.

Werder M.A., Hewitt I.J., Schoof C.G. & Flowers G.E. 2013. Modeling channelized and distributed subglacial drainage in two dimensions. Journal of Geophysical Research—Earth Surface 118, 2140–2158, doi: 10.1002/jgrf.20146.

World Glacier Monitoring Service (WGMS) 2017. Fluctuations of glaciers database. Zurich: World Glacier Monitoring Service. Accessed on the internet at https://wgms.ch/data_databaseversions/in 2018.

Xu M., Han H. & Kang S. 2017. Modeling glacier mass balance and runoff in the Koxkar River Basin on the South Slope of the Tianshan Mountains, China, from 1959 to 2009. Water 9(2), article no. 100, doi: 10.3390/w9020100.

Xu M., Yan M., Kang J. & Ren J. 2007. Progress in studies on mass balance of glaciers, Svalbard, Arctic. Journal of Glaciology and Geocryology 29, 730–737.

Xu M., Yan M., Ren J., Ai S., Kang J. & E D. 2010. Surface mass balance and ice flow of the glaciers Austre Lovénbreen and Pedersenbreen, Svalbard, Arctic. Chinese Journal of Polar Science 21, 147–159.

Zhang G., Li Z., Wang W. & Wang W. 2014. Rapid decrease of observed mass balance in the Urumqi Glacier No. 1, Tianshan Mountains, central Asia. Quaternary International 349, 135–141, doi: 10.1016/j.quaint.2013.08.035.

Zhang Y. 2011. The glacier mass balance and its relationship with climate change of Austre Lovénbreen and Pedersenbreen, Svalbard, Arctic. MS thesis, Shandong Normal University, Jinan, China.

Zhao L., Tian L., Zwinger T., Ding R., Zong J., Ye Q. & Moore J.C. 2014. Numerical simulations of Gurenhekou Glacier on the Tibetan Plateau using a full-Stokes ice dynamical model. Journal of Glaciology 60, 71–82, doi: 10.5194/tcd-7-145-2013.

Zwinger T. & Moore J.C. 2009. Diagnostic and prognostic simulations with a full Stokes model accounting for superimposed ice of Midtre Lovénbreen, Svalbard. The Cryosphere 3, 217–229, doi: 10.5194/tc-3-217-2009.