Monte Carlo modelling projects the loss of most land-terminating glaciers on Svalbard in the 21st century under RCP 8.5 forcing

The high Arctic archipelagos around the globe are among the most strongly glacierized landscapes on Earth apart from the Greenland and Antarctic ice sheets. Over the past decades, the mass losses from land ice in the high Arctic regions have contributed substantially to global sea level rise. Among these regions, the archipelago of Svalbard showed the smallest mass losses. However, this could change in the coming decades, as Svalbard is expected to be exposed to strong climate warming over the 21st century. Here we present extensive Monte Carlo simulations of the future ice-mass evolution of 29 individual land-terminating glaciers on the Svalbard archipelago under an RCP 8.5 climate forcing. An extrapolation of the 29 sample glaciers to all land-terminating glaciers of the archipelago suggests an almost complete deglaciation of the region by 2100. Under RCP 8.5, 98% of the land-terminating glaciers will have declined to less than one tenth of their initial size, resulting in a loss of 7392 ± 2481 km2 of ice coverage.


Introduction
The Arctic land-ice masses have contributed substantially to global sea level rise over the last decades [1,2] and their relative importance, with respect to other contributors to sea level rise, is expected to increase in the near future [3][4][5]. The coincidence of extensive glacierized areas and large ice volumes that is evident across the Arctic archipelagos is unprecedented outside the Greenland and Antarctic ice sheets [6]. This implies an immense leverage effect that makes Arctic glaciers and ice caps especially sensitive to climate change [7], resulting in the potential for considerable ice-mass decline and the production of large amounts of meltwater output. The surface-albedo feedback related to this future loss of land ice further fosters the Arctic amplification [8] and an Arctic summer-warming trend [9]. Indeed, the persistence of land-snow and sea-ice covers plays, in general, the most important role in this feedback mechanism [10], but sustained albedo decreases due to land-cover changes from glacierized to ice-free landscapes can be considered as an additional key control during the summer months. When modelling these changes, most scenario-based future projections currently rely on forcings that follow four representative concentration pathways (RCP), describing future climate variability from 2006 onwards [11]. While these scenarios were originally designed to represent a wide spread of emission scenarios, recent years have shown that we are currently heading to a future that even lies beyond the expectations of the highest emission scenario (RCP 8.5) [12]. In this context, Svalbard is among the Arctic regions for which the strongest air temperature increases over the 21st century are predicted [13]. Modelling studies of future glacier mass balance (the net result of mass gains minus mass losses) evolution under an RCP 8.5 forcing suggest increasingly negative balances with the most intensive drop being predicted for the southern half of the archipelago [14]. In line with this, global-scale studies on future glacier evolution and resulting sea level-rise contributions suggest that Svalbard will experience an enormous, and for the Arctic unprecedented, glacier recession over the 21st century [3][4][5].
At present, the ice cover of Svalbard extents over 33 775 km 2 and is formed by 1668 individual glaciers and ice-cap subbasins [15]. About 88% of these individual units are land-terminating glaciers (figure 1), although they solely account for 31% of the overall glacierized area. In recent years, both the relative (per unit area) and the total mass losses from Svalbard have been the lowest of all Arctic regions. The estimate of mass losses from Svalbard over 2003-2009, derived combining ICESat [16] and GRACE [2] data, is 130±60 kg m −2 a −1 [1]. For comparison, the losses during the same period for other Arctic regions were estimated to 310±40 kg m −2 a −1 for Arctic Canada North, 660±110 kg m −2 a −1 for Arctic Canada South, 210±80 kg m −2 a −1 for the Russian Arctic, and 420±70 kg m −2 a −1 for the glaciers peripheral to the Greenland Ice Sheet [1,17,18]. Regarding the total mass losses, expressed in terms of sea-level equivalent (SLE), the losses have been 0.014± 0.006 mm SLE a −1 for Svalbard, as compared with 0.091± 0.011 mm SLE a −1 for Arctic Canada North, 0.074±0.011 mm SLE a −1 for Arctic Canada South, 0.030±0.011 mm SLE a −1 for the Russian Arctic, and 0.105±0.019 mm SLE a −1 for the Greenland periphery [1,17,18]. However, the mass losses from Svalbard are predicted to increase at a sustained and faster rate than most Arctic regions [5]. This will lead to relative volume losses, over 2010-2100 and under an RCP 8.5 scenario, of about 82±18%, compared with values between 30±12% and 70±19% for the rest of Arctic regions [4].
In this paper we present area-change projections of all land-terminating glaciers on Svalbard until the end of the 21st century under an RCP 8.5 forcing, which is the RCP scenario that has been proven to best represent the observed climate warming in recent years [12]. We also estimate the associated sea levelrise contribution according to different scaling laws [5,19,20]. Moreover, we spatially resolve for different subregions of Svalbard in order to understand the variability of ice-cover changes across the archipelago. This works adds to previous studies (in particular [14]) in three main aspects: (1) we use detailed bedrock topography and ice volume from recent ice radar surveys [19] for 29 reference glaciers, which enables us to apply a high-resolution modelling framework to precisely calculate future ice retreat of these reference glaciers. (2) We include a calculation of changes in glacier topography and extent, which allows us to estimate the We apply comprehensive Monte Carlo simulations to account for potential inaccuracies and to facilitate a rigorous uncertainty assessment.

Data and methods
Our modelling approach can be split into two main parts. First, a climatic mass balance and glacier change model, which is forced by statistically downscaled climate data from ten global climate models (GCM), is applied to a set of 29 reference land-terminating glaciers for which detailed ice-thickness information is available. From each of the resulting area-change time series of the 29 reference glaciers we calculate the time elapsed until an area reduction to below 10% of the respective initial extent is reached. We consider an area reduction by >90% as a de facto glacier loss, because the glacier remnants become isolated dead-ice bodies. Second, the periods elapsed until glacier loss for all 29 reference glaciers are extrapolated to the entire set of 1471 land-terminating glacier basins on Svalbard. In what follows, we describe the main elements of the modelling, together with the associated individual uncertainties U1, U2 and U3. Details on the calculation of these uncertainties, and of the combined final uncertainty ranges for the periods elapsed until glacier loss, can be found in the Methods section in the supplementary material, along with more in-depth information on all the modelling methods applied.

Climatic mass balance and glacier change model
The basic data for the 29 reference glaciers, which spread across five of the eight subregions considered (figure 1), are given in table S1. For each of these reference glaciers, ice thickness and bedrock topography are known from recent ice-radar surveys [19]. This enables us to use a high-resolution modelling framework to precisely calculate their future retreat. In our modelling we use a spatial resolution of 30 m except for the three largest glaciers, for which a resolution of 60 m was used, and a daily temporal resolution over the period 2006-2100. For each glacier, we apply the same elevation-dependent climatic mass balance model, allowing for parameterised surface topography feedbacks to simulate ice surface lowering and thus glacier retreat. We focus on land-terminating glaciers only as the implemented feedback parameterisation is not designed to reproduce ice loss at the calving fronts of tidewater glaciers. Ablation is calculated using a temperature index model with fixed parameters for snow and ice surfaces. Accumulation is assumed to be equal to precipitation that falls at air temperatures below a certain threshold temperature. Refreezing is calculated using the P max approach [21]. Surface topography feedback is modelled using a parameterisation suggested by Huss et al [22], which is based on mass conservation, whereby surface geometry and extent are updated annually according to the calculated mass-balance distribution (see more details in SI methods). An overview of all model parameters is given in table S2.

Climate-data input and model runs
Forcing for the model runs over the period 2006-2100 is provided by statistically downscaled air temperature and precipitation data, following the RCP 8.5 scenario represented in ten different earth system or climate models (table S3) that are part of the Coupled Model Intercomparison Project Phase 5 (CMIP5). Statistical downscaling is done separately for each of the 29 reference glaciers (SI Methods). In situ reference for each glacier during downscaling is taken from 10 km resolution regional climate model (RCM) output provided by the observationally constrained European Arctic Reanalysis (EAR) product [23]. The associated climate grids were produced using the Polar Weather Research and Forecast model (version 3.1.1) and have already been successfully applied in climatic mass balance modelling on Svalbard [24]. From these grids, the respective closest grid point to the glacier is selected as in situ reference. By using RCM output as in situ reference, we indirectly also account for effects of orographic precipitation across the archipelago, an advantage that adds to the purely statistical downscaling of the GCM data.
For statistical downscaling of air temperature we apply a combination of multiple linear regression and variance inflation techniques [25,26]. For statistical downscaling of precipitation the local scaling method in combination with multiple linear regression [27,28] is used. Elevational distribution of both variables across the glacier is performed by applying fixed lapse rates (table S2).
The modelling framework consisting of downscaling and mass-balance modelling is tuned to massbalance observations at Austre Lovénbreen and Middre Lovénbreen from the period 2006-2010 by fitting a temperature offset for adjustment of the EAR air temperature grids (figure S2). These two glaciers were selected because they are the only ones of the 29 reference glaciers with observed mass balances available. All other parameters of the modelling framework are preselected according to appropriate values given in the literature (table S2). No spatial variability of the parameters is introduced because of a complete lack of in situ references that could have been used for calibration.

Uncertainty assessment: Monte Carlo simulations
As in situ observations of mass balance or climate data across our 29 glaciers are very scarce, we use the same model-parameter sets for all reference glaciers and apply comprehensive Monte Carlo simulations to account for potential inaccuracies and facilitate a rigorous uncertainty assessment. It consists of 900 individual runs for each of the ten different climate forcings. This number represents the minimum number of runs necessary to achieve stable results (figure S3). During each run, all model parameters were varied within their individual probability densities, representing both systematic and random uncertainties (table S4). Systematic uncertainties form constant deviations that vary from run to run, while random uncertainties are formed by deviations that vary continuously from day to day along the runs. The underlying probability density functions were predefined by appropriate standard deviations. In addition to model parameters, the air temperature and precipitation time series employed as model forcing were also altered according to systematic and random uncertainties (table S4). As direct result of the Monte Carlo simulation for each glacier, a mean time series of area change as well as its one-sigma uncertainty range is produced ( figure S4). These are based on data from 9000 different runs, i.e. ten different forcings with 900 Monte Carlo runs each. This one-sigma range is the first and main type of uncertainty (U1). It represents the uncertainty introduced by the spatially invariant model parameters and by the statistical downscaling of the climate data. As an example, figure S1 presents area changes with time according to the mean of all 9000 runs for one arbitrary reference glacier.

Extrapolation to the entire set of Svalbard landterminating glaciers
We extrapolate the modelling results of the reference glaciers to the whole set of 1471 land-terminating glaciers on Svalbard. To achieve more reliable results during the extrapolation, we concentrate on readily accessible predictors related only to the surface of the glaciers rather than relying on uncertain, modelled volume data. We relate the periods elapsed until glacier loss to five surface-related topographic variables by using multiple regression analysis. The five topographic variables are initial glacier area, median, minimum and maximum elevation of the glacier and the standard deviation of elevations across the glacier area. We apply a leave-one-out cross validation to calculate the regression, which allows us to quantify both the uncertainties related to the extrapolation and to the regression itself which, in turn, facilitates a full propagation of the Monte Carlo uncertainties. Leaveone-out cross validation means that we calculate the multiple regression on the basis of 28 reference glaciers and derive its uncertainty in terms of the root mean square error of the period until glacier loss by comparison with the 29th glacier. This procedure is repeated 29 times with iteratively leaving out one glacier. The final parameters of the multiple regression including their one sigma uncertainty ranges are then calculated from the respective 29 individual values (table S5). This second type of uncertainty (U2) is the uncertainty introduced by the extrapolation of the results from the 29 reference glaciers to all other landterminating glaciers. The third and last type of uncertainty (U3) is the uncertainty introduced by the representation of the period elapsed until glacier loss by the topographic variables. It is given by the mean root mean square error derived from the crossvalidation procedure.

Results and discussion
The area-change projections of all land-terminating glaciers on Svalbard by the end of the 21st century under an RCP 8.5 forcing are given in table 1, and their associated sea level-rise contributions, estimated according to different scaling laws [5,19,20], are presented in table 2. Our projections show that Svalbard will lose 1435±59 of its 1471 land-terminating glaciers by the end of the 21st century, which represents a relative reduction in terms of number of 97.6±4.0% (table 1). The fact that the uncertainty range extends beyond the total number indicates that even a complete loss of all land-terminating glaciers is potentially possible. The relative glacier losses reach 69.9±23.8% when estimated in terms of glacierized area (table 1), or 55%-58% in terms of volume, according to the scaling laws used for converting to sea level equivalent (table 2).
When resolving for the eight subregions (SR) considered the picture becomes a bit more diverse. Three subregions are predicted to lose all their land-terminating glaciers already well within the 21st century (table 1, figure 1). While Barentsøya and Edgeøya (SR8) are projected to already show a complete loss by 2080, Nordenskiöldland (SR4) and Vestfonna (SR6) follow closely after (figure 2). The two subregions whose land-terminating glaciers show the strongest potential for surviving the 21st century are Austfonna (SR7) and Northeast Spitsbergen (SR3), which are the subregions with either the largest and/or the highest individual glacier basins (figures 1 and 3). The higher a glacier is situated or the larger it is, the later it will be at risk of vanishing completely. The elevation range of the glacier determines the timing of frequent exposure to melting conditions, while the size of the glacier simply controls the amount of ice that needs to be melted until a loss of the ice body occurs. This linkage is also supported by the glacier evolution in SR4 and SR8. In both subregions, low-lying land-terminating glaciers dominate the setting, but the respective size distributions of the glaciers are completely different. While in Nordenskiöldland the smallest mean glacier size of all subregions (2.8 km 2 ) is found, on Barentsøya and Edgeøya the glaciers are considerably larger (18.6 km 2 ) but nevertheless vanish faster than in any other subregion ( figure 2). Interestingly, the glaciers of and around Vestfonna (SR6), which at first glance should react similarly to those of Austfonna, seem to experience a different fate. The large land-terminating basins of Vestfonna already start to vanish from the 2060s onwards. This results in a rather early complete disappearance before the end of the 21st century even if most glaciers of this size show distinctly stronger persistency in other subregions ( figure 3). Hence, saying that size matters is not a priori justified or at least not generally valid. Elevation is the key control for the timing of glacier loss over the low-lying share of landterminating glaciers. This inference is supported by a Table 1. Overview on glaciers and glacier losses on Svalbard divided into eight subregions. Glacier losses are given by the end of the 21st century and under an RCP 8.5 forcing. The subregion numbers given in column one correspond to the numbering shown in figure 1. The numbers of reference glaciers that are explicitly modelled in the Monte Carlo simulations are shown in column two and the mean area per individual glacier in column three. Numbers (upper value) and total extent in km 2 (lower value) of all glaciers regardless of type are shown in column four. The shares of land-terminating glaciers (LTG) are given as both absolute (column five) and relative (column six) numbers and extents. The shares of LTG lost by 2100 are given accordingly (columns seven and eight) but together with the overall uncertainty estimates. These estimates combine the uncertainties induced by the Monte Carlo simulations, by the extrapolation of the individual results to the entire archipelago and by the multiple regression that relates the period until glacier loss with different topographic parameters.   recent study that predicts a loss of Vestfonna's accumulation area and thus the start of a final downwasting of the ice cap already for the 2040s if an RCP 8.5 forcing is assumed [29]. Over recent years, the ice-mass losses across the eastern parts of the archipelago were found to be rather limited [16]. However, when taking our results together it becomes obvious that the rate of these mass losses is expected to increase substantially over the oncoming decades ( figure 1), a fact that is in line with findings by Førland et al [30], who predicted the highest future air temperature increases to occur over the eastern parts of Svalbard. By contrast, the recent positive glacier surface-elevation changes over NE Spitsbergen [16] seem to persist during the near future. In particular, the intensively glacierized ridge along the northern part of this subregion shows land-terminating glaciers that will survive until the late 22nd century (figure 1). Only 46.1±20.2% of the glacierized area in NE Spitsbergen is predicted to be lost by the end of the 21st century, which is the smallest share among all subregions (table 1).
Calculations of the sea-level rise going along with these large projected glacier losses show values ranging between 2.00±1.26 mm and 2.50±1.67 mm when integrated over the entire archipelago (table 2). NE Spitsbergen (SR3) and the area around Vestfonna (SR6) show the largest contributions to this total while Andrée & Dickson Land (SR2) and Nordenskiöldland (SR4) show the smallest contributions. These numbers of SLE are rather small compared with the currently assumed total glacier volume on Svalbard (13-24 mm SLE [19]) even if the land-terminating glaciers lost by the end of the 21st century hold a share of 21.9±7.5% of the currently glaciated area on the archipelago (table 1). This seeming inconsistency is, however, explained by the strongly nonlinear relationships between area and volume of the glaciers that are assumed for the calculations of SLEs (table 2). The tidewater glaciers on the archipelago are, on average, considerably larger than the land-terminating glaciers and thus hold a higher proportion of the total volume. However, the almost complete deglaciation of certain subregions of Svalbard before the end of the 21st century (figure 1) clearly indicates the relevance of the rather small land-terminating glaciers at least for the local-scale albedo feedback. These glaciers show a scattered distribution thus increasing the integrated albedo over the entire area, as discussed by Lang et al [14]. Hence, their disappearance will considerably change the local-scale radiation balance and thus likely the entire local-scale climate. The induced increases of local air temperatures will then trigger a feedback loop leading to increased glacier melt. This scenario can be considered as an early precursor of what might happen to similar Arctic glacierized regions.
When comparing our projections with those stemming from regionalized global mass-balance projections [3][4][5], we have to keep in mind that our results are specific for land-terminating glaciers, while the mentioned studies consider the entire set of Svalbard glaciers, i.e. including tidewater glaciers. Second, we note that the relative volume losses that they present refer to differing estimates of the total volumes of Svalbard glaciers, so they are not directly comparable among each other. This is why we have recalculated the relative losses estimated by those authors, referring all to two common total volume estimates of Svalbard glaciers. We selected for these a low-range estimate of 16.6 mm SLE, based on volume-area scaling [19], and a high-range one of 19.9 mm SLE, based on inversion from surface geometry and mass balance data [6] but recalculated for the newest Randolph Glacier Inventory 5.0 as done in [4].
Marzeion et al [3] projected, under RCP 8.5 forcing, 14.0 mm SLE total mass losses from Svalbard glaciers during 2006-2100, which is a 70%-84% reduction from the current ice volume of Svalbard glaciers. Radić et al [5] estimated larger losses of 15.8 mm SLE, for the same period and scenario, i.e. a 79%-95% volume reduction. Although both of these estimates include tidewater glaciers, their estimated losses only consider surface mass losses (as we do), and neglect frontal ablation (the combination of iceberg calving and submarine melting at the marine-terminating glacier fronts). Frontal ablation is included in the estimate of 16.4 mm SLE (i.e. 82%-99% volume reduction) by Huss and Hock [4]. However, the latter losses translate into an effective increase in sea level of 13.9 mm when accounting for ice located below sealevel presently displacing ocean water [4]. Compared with these results, our estimated losses of 2.0-2.5 mm SLE (according to the scaling laws used; see table 2) could seem insignificant in absolute terms, but not when considered as volume losses relative to the total volume of Svalbard land-terminating glaciers, which amount to 55%-58%. On the other hand, two other facts reinforce our estimates. First, smaller glaciers, as is the case of land-terminating glaciers in Svalbard, are expected to experience comparatively larger mass losses [4]. Second, accelerated melt rates have been identified as the main component driving changes in net mass loss of Svalbard glaciers [5].
Regarding shortcomings of our approach, we note that about two thirds of the 29 reference glaciers are located in one single subregion, which introduces further but unquantifiable uncertainty into the final modelling results, beyond those already considered in section 2 (U1, U2 and U3). During the extrapolation of the modelling results of the individual reference glaciers, the local pattern of future air temperature and precipitation trends in Nordenskiöldland (SR 4) is transferred to the entire archipelago. This induces uncertainties in two different ways. First, the substantial spatial variability across the archipelago, which exists regarding future climate trends [30], is partly masked. Second, the set of reference glaciers consists of valley and cirque glaciers which are small when compared to land-terminating glacier units in SR 3, SR 6 and SR 7 (figure 1). This implies neglecting certain mass balance-relevant processes. Snow redistribution by wind drift is especially important across the individual land-terminating glacier units of the large ice fields and ice caps located in SR 3, SR 6 and SR 7 [e.g. 31] and is able to significantly alter the mass balance of these units by mass transfers to or from neighbouring glacierized areas. Moreover, these glacier units are characterised by large areas that are situated above the regional equilibrium line [24], and thus anyway show completely different accumulation characteristics [e.g. 32], which can hardly be captured by extrapolation of modelling results from the much smaller reference glaciers.
We also note that several studies have suggested that temperature-index models might be oversensitive to temperature change and thus overestimate future glacier melt [33,34]. In a recent study [4], it has been estimated that using a simplified energy-balance model instead of a temperature index model in predicting glacier evolutions over the 21st century under an RCP 8.5 scenario results, in the case of Svalbard, in a reduction of the projected mass loss by 24%. Consequently, our results could be overestimating future glacier reductions by similar amounts. Nevertheless, global glacier-mass projections available for Svalbard [3][4][5] use temperature-index models instead of full energy balance models.
Our approach stands out from the rest by using detailed bedrock topographies and accurate initial glacier volumes. As discussed in Huss and Hock [4], one of the main shortcomings of global glacier-mass projections is the uncertainty of the initial glacier volume. In particular, they find that Svalbard glaciers present the largest sensitivity to the initial volume among all the studied regions, i.e. a ±30% deviation in the initial volume results in a −24% to +29% change in the SLE contribution.

Conclusions
We modelled the evolution of 29 individual Svalbard land-terminating glaciers over the 21st century by applying a parameterising mass-balance model embedded in a Monte Carlo simulation environment. Climate forcing for these projections was provided by downscaled data from ten different GCMs, all representing the RCP 8.5 future climate scenario. This modelling extends beyond already published studies by including glacier-retreat calculations that account for explicit knowledge on glacier bedrock topography and ice thickness. The results for these 29 glaciers were then extrapolated to all land-terminating glaciers on Svalbard using multiple regression, with glacier topography-related parameters as predictors. This allowed us to determine the periods elapsed until glacier loss for all land-terminating glaciers on Svalbard under an assumed RCP 8.5 scenario.
Our projections suggest that by the end of the 21st century 1435±59 of the 1471 land-terminating glaciers on Svalbard will be lost. The associated area reduction amounts to 7392±2521 km 2 , i.e. a loss of ∼70±24% with respect to the initial area of 10 575 km 2 . Due to the nonlinear relationship between glacierized areas and glacier-ice volumes, the projected overall volume loss ranges only between 55±36% and 58±37%, depending on the volumearea scaling law considered. The associated sea levelrise contribution was estimated to range between 2.00±1.26 and 2.50±1.67 mm SLE.
We found that the low-lying glaciers of Nordenskiöldland (SR4) and Barentsøya and Edgeøya (SR8) will be first to disappear. By contrast, the large, and high-lying land-terminating glaciers along the main ridge of NE Spitsbergen (SR3) are projected to partly last until the middle of the 22nd century.