Accelerating future mass loss of Svalbard glaciers from a multi-model ensemble

Abstract

Projected climate warming and wettening will have a major impact on the state of glaciers and seasonal snow in High Arctic regions. Following up on a historical simulation (1957–2018) for Svalbard, we make future projections of glacier climatic mass balance (CMB), snow conditions on glaciers and land, and runoff, under Representative Concentration Pathways (RCP) 4.5 and 8.5 emission scenarios for 2019–60. We find that the average CMB for Svalbard glaciers, which was weakly positive during 1957–2018, becomes negative at an accelerating rate during 2019–60 for both RCP scenarios. Modelled mass loss is most pronounced in southern Svalbard, where the equilibrium line altitude is predicted to rise well above the hypsometry peak, leading to the first occurrences of zero accumulation-area ratio already by the 2030s. In parallel with firn line retreat, the total pore volume in snow and firn drops by as much as 70–80% in 2060, compared to 2018. Total refreezing remains largely unchanged, despite a marked change in the seasonal pattern towards increased refreezing in winter. Finally, we find pronounced shortening of the snow season, while combined runoff from glaciers and land more than doubles from 1957–2018 to 2019–60, for both scenarios.

Introduction

Due to sea-ice related feedbacks, the climate on the archipelago of Svalbard, situated at the southwestern boundary of wintertime Arctic sea-ice, experiences amplified climate change compared to the global mean (AMAP, 2017; IPCC, 2019). During 1971–2017, temperatures in Svalbard have risen by 3–5^°C (Hanssen-Bauer and others, 2019), with fastest changes in northern and eastern Svalbard. At the same time, precipitation trends at weather stations in recent decades have not been significant or show moderate positive significant trends (Hanssen-Bauer and others, 2019; Førland and others, 2020). Ensemble future projections from statistical downscaling and regional climate models (RCMs) indicate a 6–7^°C warming and a precipitation increase of on average 45% from 1971–2000 to 2071–2100 for a medium RCP 4.5 emission scenario. With strongest winter and spring warming, a main characteristic of the projected precipitation increase is a substantial rise in frequency and intensity of winter-time rainfall events (Hansen and others, 2014; Vikhamar-Schuler and others, 2016; Bintanja and Andry, 2017).

Air temperature and precipitation trends are the main drivers of glacier mass change. Glacier mass change resulting from atmosphere–surface–subsurface interactions is best described by the climatic mass balance (CMB), determined by precipitation, moisture fluxes and runoff (Cogley and others, 2001). In contrast to surface mass balance (SMB), commonly observed using the glaciological method, CMB accounts for internal accumulation, i.e. refreezing below the previous summer surface. A number of recent modelling studies have quantified the past SMB/CMB of all glaciers in Svalbard (Lang and others, 2015b; Aas and others, 2016; Østby and others, 2017; Möller and Kohler, 2018; Van Pelt and others, 2019; Noël and others, 2020), which cover an approximate area of 33 775 km² or 57% of the total land area (Fig. 1). These studies consistently show a close to zero or weakly negative SMB/CMB in recent decades, but with a clear negative trend (Schuler and others, 2020). Schuler and others (2020) combined observational records from various regions in Svalbard to estimate a total CMB of −7 ± 4 Gt a⁻¹ during 2000–2019, and a total glacier mass balance, i.e. including frontal ablation, of −8 ± 6 Gt a⁻¹.

Fig. 1. Topographic map of Svalbard. Different colour scales are used for glacier- and non-glacier areas. The archipelago is split into three regions, with borders given by the brown lines: northwest (NW), northeast (NE) and south (S). UTM easting and northing coordinates are in zone 33X. Elevations are extracted from the DEM S0 Terrengmodel Svalbard, provided by the Norwegian Polar Institute (Norwegian Polar Institute, 2014). Glacier outlines are extracted from the Global Land Ice Measurements from Space (GLIMS) database (König and others, 2014). The orange circle marks the location of the time series shown in Figure 3.

Although many modelling or observational studies present past and current CMB and SMB of all Svalbard glaciers, future projections are relatively scarce. Day and others (2012) used a statistical approach to find dramatic increases in precipitation offsetting most of the increase in glacier melting during the 21st century. Lang and others (2015a) used a RCM to simulate SMB evolution up to 2100 under a high-emission RCP 8.5 scenario to find that the entire glacierized area of Svalbard would be prone to a negative SMB from 2085 onwards. Most recently, RCM output for the time-slices 1970–2000 and 2070–2100, included in the Climate in Svalbard 2100 report (Hanssen-Bauer and others, 2019), showed a drop in Svalbard CMB from −14 to −87 Gt a⁻¹ and an equilibrium line altitude (ELA) increase from 653 to 1067 m a.s.l. between the two periods for the RCP 8.5 scenario.

Similar to glacier mass balance, seasonal snow conditions on glacierized and non-glacierized terrain in Svalbard depend heavily on trends in temperature and precipitation, which control precipitation type and amount, snowmelt and the length of the snow and snow-free season. Modelling studies covering all of Svalbard show that since the 1950s snow arrives later in autumn by 1–2 d per decade, while no significant trend is found for snow disappearance in spring/summer (Van Pelt and others, 2016a, 2019). Future consequences of projected climate change on snow conditions in Svalbard include, among others, continued snow season shortening, an earlier peak of maximum snow depth in spring, increased avalanche risk and an increased frequency and intensity of winter-time rain-on-snow events (Hanssen-Bauer and others, 2019).

This study is a follow up on a recent study by Van Pelt and others (2019), presenting a dataset of CMB, snow conditions and runoff in Svalbard for the period 1957–2018. Here we force our simulations by using output from the Arctic CORDEX RCMs (Hanssen-Bauer and others, 2019), which downscale several different general circulation model (GCM) simulations. We extract climate trends and combine them with historical downscaled reanalysis data to construct a future climate forcing for two emission scenarios (RCP 4.5 and RCP 8.5). With these two forcings a coupled energy-balance – snow/firn model (Van Pelt and others, 2012, 2019) simulates glacier CMB, seasonal snow conditions and runoff for the period 2018–60. We employ the same spatial resolution (1 × 1 km), fixed geometry, time step (3-hourly) and parameter setup as in Van Pelt and others (2019). The latter has extensively been calibrated against a range of observational data across Svalbard (Van Pelt and others, 2019). The 1-km spatial resolution used provides unprecedented spatial detail of future glacier mass change and snow conditions in Svalbard; for comparison, Day and others (2012), Lang and others (2015a) and Hanssen-Bauer and others (2019) used spatial resolutions of 25, 10 and 2.5 km respectively.

Methods

Model and setup

A surface energy-balance model coupled with a multilayer snow/firn model (EBFM; Van Pelt and others, 2012, 2019) is used to simulate the future CMB, seasonal snow and firn conditions, and runoff for the entire land-area of Svalbard. Detailed descriptions of the model physics and numerical implementation can be found in Van Pelt and others (2019) and references therein.

The model uses meteorological input of air temperature, relative humidity, cloud cover, precipitation and air pressure, to solve the surface energy balance and estimate surface temperature and melt. The energy-balance routine determines net solar radiation, net thermal radiation, turbulent heat exchange and the subsurface heat flux (Van Pelt and others, 2012). A recently updated albedo scheme (Van Pelt and others, 2019), based on Bougamont and others (2005), introduces the effect of temperature and water content dependent ageing of snow on surface albedo. The fraction of precipitation falling as snow and rain changes linearly within a ± 1^°C band around a calibrated threshold temperature (0.6^°C; Van Pelt and others, 2019). The surface model is two-way coupled to a subsurface routine, computing the multi-layer evolution of density, temperature and water content. The subsurface model simulates densification, due to refreezing and compaction, as well as temperature changes, due to heat conduction and refreezing. The model accounts for irreducible water storage, i.e. water held in place by capillary forces, while percolating water is transported vertically using a combined bucket-type and deep-percolation scheme (Marchenko and others, 2017), and becomes runoff if it reaches the base of the vertical domain or if the underlying layers are impermeable ice or soil. Seasonal snow modelling in non-glacier terrain became possible with the inclusion of a soil thermal model (Westermann and others, 2011), as described in Pramanik and others (2018).

After earlier applications of the model to single or small sets of glacier basins in central Svalbard (e.g. Van Pelt and others, 2012, 2018; Marchenko and others, 2017) and western Svalbard (e.g. Van Pelt and Kohler, 2015; How and others, 2017; Pramanik and others, 2018; Winsvold and others, 2018; Köhler and others, 2019), the model has recently been applied to the entire glacierized and non-glacierized area of Svalbard (Van Pelt and others, 2019). In the latter, the model and the downscaling of meteorological fields were extensively evaluated against in situ stake mass-balance data, weather station records and firn core profiles across Svalbard, which helped to reduce model biases and errors.

Here, we use the same calibrated set of parameter values as in Van Pelt and others (2019) for future simulations. Using the 3-hourly meteorological time series described in the next subsection, we will perform simulations for the emission scenarios RCP 4.5 and RCP 8.5 covering the period from 1 September 2018 to 31 August 2060. Subsurface conditions at the end of the historical simulation (31 August 2018) in Van Pelt and others (2019) are used as initial conditions in the current experiments. As in Van Pelt and others (2019), the model is run with a 3-hourly temporal and 1 × 1-km spatial resolution, and the subsurface model consists of 50 vertical layers on a Lagrangian grid, extending to maximum depths of 5 and 20 m in non-glacier and glacier terrain, respectively.

As in Van Pelt and others (2019), glacier outlines, extracted from the Global Land Ice Measurements from Space (GLIMS; König and others, 2014) database, and grid elevations, extracted from the S0 Terrengmodel Svalbard (Norwegian Polar Institute), are held fixed throughout the simulations. This arguably has an impact on e.g. simulated CMB and runoff, even though the two main effects (glacier thinning and retreat) have impacts of opposite sign. To test this, we perform an additional sensitivity experiment where the model simulates only large (>20 km²) land-terminating glaciers and converts annual CMB values into elevation changes using a parameterization given in Huss and others (2010) and Huss and Hock (2015). Ice thickness is taken from the Svalbard glacier thickness database presented in Fürst and others (2018). The simulation is performed for the RCP 4.5 scenario and covers the period 2018–60; elevations and outlines in 2018 are given by the original (previously fixed) datasets. We do not attempt to simulate tidewater glacier retreat, to avoid additional uncertainties related to the description of frontal ablation, and small (<20 km²) land-terminating glaciers are excluded due to a lack of 1 km² grid cells per glacier. Glacier thinning/thickening is applied at the end of every negative/positive CMB year on 31 August while grid cells that reach zero ice thickness are no longer considered part of the glacier and the mask is updated accordingly.

Climate input

Our model experiments require spatially distributed and time-dependent meteorological input, which includes the meteorological parameters of air temperature, relative humidity, cloud cover, precipitation and air pressure. These time series are constructed using a combination of historical data, describing short-term fluctuations, and future RCM projections, from which long-term trends are extracted. This involves the following steps.

We extract long-term trends of air temperature and precipitation from output of five RCMs, contained in the Arctic CORDEX database (http://climate-cryosphere.org/activities/polar-cordex/arctic). Arctic CORDEX is an international initiative that provides atmospheric data for the Arctic region from an ensemble of RCM results. Five models included Svalbard in their domain and generated output for both RCP 4.5 and RCP 8.5 emission scenarios for 1950–2100. The model realizations differ in the use of RCMs and forcing from GCMs. The same set of RCM model runs is used here as in the Climate in Svalbard 2100 report (Hanssen-Bauer and others, 2019), except that we exclude output from three models that were only run for the RCP 8.5 but not for the RCP 4.5 scenario. For more details about the individual RCM simulations, we refer to Table A1.1 in Hanssen-Bauer and others (2019).

For every RCM model run, we use linear regression to calculate the temperature and precipitation trends at each 50 × 50 km grid cell, for each month. The resulting trends are downscaled to our 1-km grid by means of natural neighbour interpolation. The spread of the models for temperature and precipitation is shown in Appendix Figures 14 and 15 respectively, showing annual mean trend distributions. All considered models predict positive temperature and precipitation trends for all of Svalbard; four model realizations show trends of comparable magnitude, with one model realization showing stronger trends for both temperature and precipitation. Trends from the five RCM runs are averaged to generate single maps of trend distributions per month. The above steps are applied for both RCP scenarios, generating single monthly maps of temperature and precipitation trends, shown in Figures 16 and 17 respectively for the RCP 4.5 scenario. Strongest warming is projected for November, December and January, while the summer months (June–August) exhibit weakest trends in both temperature and precipitation. Annually averaged trend distributions for temperature and precipitation for RCP 4.5 and RCP 8.5 scenarios (Fig. 2) reveal a southwest-northeast gradient in warming, with strongest warming (up to 1.2^°C decade⁻¹ for RCP 8.5) and precipitation increase (up to 8.3% per decade) in northeast Svalbard. Spatially averaged warming increases modestly from 0.9 to 1.1^°C decade⁻¹ from RCP 4.5 to RCP 8.5, while the spatial mean precipitation trend increases from 4.7 to 5.9% decade⁻¹.

Fig. 2. Model-averaged annual temperature (left) and precipitation (right) trends for RCP 4.5 (a-b) and RCP 8.5 (c-d) for 1988–2060.

Next, the monthly temperature and precipitation trend maps for the two emission scenarios are superimposed on the 3-hour meteorological time series. For this we select the final 30 years (1988–2018) of meteorological input that was previously used to force the model in Van Pelt and others (2019). The meteorological data are based on a simulation with the High Resolution Limited Area Model (HIRLAM; Reistad and others, 2011) forced with European Centre for Medium Range Weather Forecasts (ECMWF) reanalysis data (Uppala and others, 2005). The HIRLAM data are downscaled from the original ~10 km resolution to the 1-km resolution using interpolation and elevation corrections as in Van Pelt and others (2019). The 1988–2018 time series of temperature and precipitation are first detrended by subtracting linear trends on a monthly basis (e.g. for January the linear trend for all January data between 1988 and 2018 is subtracted). Then, the 30 year detrended time series are repeated during 2018–48 and the first 12 years of the record are repeated again during 2048–60.

Finally, the monthly linear trend maps are superimposed on the detrended 3-hourly time series of temperature and precipitation. Trends are applied from 1 September 2018 onwards, and calculated relative to the mid-point of the reference period in 2003. For the other meteorological parameters (relative humidity, cloud cover and air pressure), the 1988–2018 time series are directly repeated during 2018–48 and 2048–60, and no detrending or long-term trends have been applied. Figure 3 shows an example of 1-year temperature time series with daily resolution during 1999, 2029 and 2059 at a site at 520 m a.s.l. on Nordenskiöldbreen in central Svalbard (location marked in Fig. 1), revealing the seasonal inhomogeneity in warming.

Fig. 3. Time series of air temperature in 1999, 2029 and 2059 at 520 m a.s.l. on Nordenskiöldbreen in central Svalbard (location marked in Fig. 1).

This approach, combining short-term historical variability with long-term trends from ensemble RCM output to generate future time series, has been used previously by Van Pelt and others (2016b). The technique used here is also similar to the Statistical Analogue Resampling Scheme (STARS; Orlowsky and others, 2007; Lutz and Gerstengarbe, 2014). Both methods superimpose linear trends on historical variability, but differ in how short-term weather variations are repeated. In general, the main advantages of the employed method are that only one run is required with the energy balance–snow model per RCP scenario (rather than running the model for multiple individual future RCM realizations), and that potential erroneous short-term meteorological variability in individual future RCM realizations is not affecting our results since we only extract long-term trends from the RCM output. A potential drawback is a lack of trends in cloud cover, relative humidity and air pressure. Even though there are indications that future cloud cover and relative humidity may increase (e.g. Cao and others, 2017), these processes depend heavily on local circulation patterns and orography and as a result no reliable trends could be extracted from the coarse-resolution Arctic CORDEX RCMs. The coarse resolution of the Arctic CORDEX RCMs currently also limits the representation of feedbacks from the surface to the atmosphere over land, e.g. regarding local atmospheric warming due to albedo reductions. Higher-resolution RCM output would facilitate generation of more detailed air temperature trend maps, exhibiting more local variability than in the current maps shown in Appendix Figures 14 and 16. Finally, by directly repeating historical variability in the future, we implicitly assume that the frequency of weather events remains consistent in time. This could potentially introduce biases in our model results, e.g. when the frequency of warm spells changes. However, a recent study by Isaksen and others (2016) did not detect a significant change in the frequency of ‘warm’ atmospheric circulation types, and ascribed the recent warming rather to increased sea surface temperature and a general background warming trend.

Results and discussion

In this section, we present, analyse and discuss the results from two future simulations for the scenarios RCP 4.5 and RCP 8.5, and compare results with existing output for the period 1957–2018 (Van Pelt and others, 2019). Most emphasis will be given to the results for the medium-emission RCP 4.5 scenario. Finally, we discuss uncertainties and the sensitivity of the results to the use of a fixed glacier geometry.

Climatic mass balance

Van Pelt and others (2019) found that the mean CMB for Svalbard was +0.09 m w.e. a⁻¹ for 1957–2018, with a negative trend of −0.06 m w.e. a⁻¹ decade⁻¹. Here, we find that for 2019–60 the spatial mean CMB becomes strongly negative with values of −0.82 and −1.06 m w.e. a⁻¹ for RCP 4.5 and 8.5 respectively (Fig. 4). This corresponds to a total mass change of −31 and −40 Gt a⁻¹ to which mass loss due to frontal ablation needs to be added (assuming negligible basal melting). Blaszczyk and others (2009) estimated frontal ablation at −0.20 m w.e. a⁻¹ (− 6.8 Gt a⁻¹) for the period 2003–09. An updated frontal ablation estimate for the more recent period 2013–18 is considerably higher, at ~15 Gt a⁻¹ (G. Moholdt and others, manuscript in preparation). Differences between 1957–2018 and 2019–60 CMB for RCP 4.5 (Fig. 4c) show strongest future changes at low elevations in the south and east of Svalbard (down to −2 m w.e. a⁻¹); this spatial distribution is in line with the pattern Van Pelt and others (2019) found for 1957–2018. Time series of CMB for 1957–2060 (Fig. 5a) reveal an acceleration of mass loss despite linear temperature and precipitation trends, illustrating the impact of a feedback resulting from ablation zone expansion on the CMB (Van Pelt and others, 2012; Van Pelt and Kohler, 2015).

Fig. 4. Spatial distribution of CMB, averaged over the periods 1957–2018 (a) and 2019–60 (b) for the RCP 4.5 emission scenario. The difference between the two periods (Δ) is shown in panel (c).

Fig. 5. Annual time series of total CMB (a) and porespace (b) on glaciers, together with total runoff (c) and refreezing (d) for glacier and non-glacier areas. Time series are shown for the historical run (1957–2018) and for the RCP 4.5 and 8.5 future scenarios (2018–60).

The difference in future CMB between the two emission scenarios is only 21% of the total change between the historical run and the mean of the RCP 8.5 scenario, which is comparable to the relative changes in both temperature and precipitation. It is noteworthy that temperature and precipitation discrepancies between the two RCP scenarios are small for the period considered here compared to projections for 2071–2100, as further discussed in Hanssen-Bauer and others (2019).

Increasingly negative CMB implies a rise in the ELA. On average for Svalbard, we find an increase in ELA from 372 m a.s.l. (1957–2018) to 581 m a.s.l. for the RCP 4.5 scenario (2019–60). When splitting Svalbard into three regions (Fig. 1), we find the lowest ELA in northeast Svalbard, which experiences the smallest increase, from 324 to 513 m a.s.l., between the historical and future periods (Fig. 6). Although the highest ELAs are in northwest Svalbard, the strongest increase is in southern Svalbard (from 394 to 644 m a.s.l.). The mean 2019–60 ELA for southern Svalbard lies way above the peak in hypsometry (Fig. 6e), while for NW and NE Svalbard the ELA is at or just above the hypsometry peak (Figs 6a, c). As argued by Noël and others (2020), the hypsometry of Svalbard glaciers, which lies predominantly at relatively low elevations, increases their sensitivity to future climate warming. With increasing ELA the accumulation-area ratio (AAR), i.e. the ratio of the accumulation zone area to the total glacier area, declines. This effect is most noticeable in southern Svalbard, and is amplified by the comparatively sharp peak in hypsometry. Our results suggest that the AAR is likely to reach zero with increasing frequency from the 2030s onwards, leading to negative CMB for the entire glacier area. The higher elevations of NW and NE Svalbard mean that zero AAR values occur later, in the 2050s. The higher climate sensitivity of CMB in low-elevation ablation zones compared to high-elevation accumulation areas (Van Pelt and others, 2012) explains the general steepening of the mass-balance gradient (Fig. 6a, c, e), which is most pronounced in NE Svalbard.

Fig. 6. Left: Elevation profiles of glacier hypsometry and CMB for northwest (a), northeast (c) and south Svalbard (e) for the RCP 4.5 scenario. Right: Time series of ELA and AAR for northwest (b), northeast (d) and south Svalbard (f) for the RCP 4.5 scenario. The regions (NW, NE and S) are indicated in Figure 1.

Glacier snow and firn conditions

Rising ELA implies that the firn line retreats to higher elevations, while the remaining firn densifies as an indirect response to higher melt rates and increased significance of refreezing. Both firn line retreat and snow/firn densification affect the amount of pore volume, i.e. the amount of air, in snow and firn across Svalbard, which we use here as a parameter for further analysis. Here, pore volume is defined between the surface and a depth of 12 m below the surface. Annual snapshots of pore volume, expressed in m³ m⁻², in 2018 and 2060 (Fig. 7a, b), show a complete loss of pore volume in southern Svalbard by 2060 for the RCP 4.5 scenario, and dramatic pore volume reductions on Holtedahlfonna and Austfonna Ice Caps in northwest and northeast Svalbard, respectively; only high-elevation Lomonosovfonna Ice Cap in central Svalbard maintains a substantial firn area by the end of the simulation. Time series of pore volume (Fig. 5b) show a rapid reduction in pore volume, dropping from 68 km³ in 2018 to 21 km³ (RCP 4.5) or 14 km³ (RCP 8.5) in 2060.

Fig. 7. Spatial distribution of snow/firn pore volume on glaciers in 2018 (a) and 2060 (b) for the RCP 4.5 scenario. The difference between the two distributions (Δ) is shown in panel c.

Refreezing has been shown to be most prominent in the accumulation zones on Svalbard (e.g. Østby and others, 2017; Van Pelt and others, 2019). Despite predicted reductions in the accumulation area and pore volume, refreezing does not show a significant upward or downward trend (Fig. 5d) in either of the future scenarios. Total refreezing for all of Svalbard changes from 12.5 Gt a⁻¹ for 1957–2018 to 12.2 Gt a⁻¹ for both RCP scenarios, and refreezing on glaciers alone increases weakly from 8.6 to 9.0 Gt a⁻¹. On glaciers, a combination of factors explains the absence of a clear trend. On the one hand, there is reduced potential for refreezing in snow and firn due to (1) less potential for water storage in (denser) snow and firn, and (2) less pronounced cooling of snow and firn during winter reducing the potential for refreezing at the start of the following melt season. On the other hand, an increased frequency of winter melt and rainfall events across Svalbard (Fig. 8a) increases the potential for winter-time refreezing. The above factors appear to compensate each other in a future climate. One consequence is that the seasonal distribution of refreezing, as shown in Figure 8b for all of Svalbard, changes markedly, revealing reduced refreezing during the melt season and increased refreezing during the cold season. Comparatively high winter-time refreezing rates in early years of the future simulation are explained by the fact that after removing the linear trend the winters from 1990 to 1994 in the reference period remain above-average warm, with a relatively high frequency of winter warm spells; repeating these time series with future climate trends applied to them yield relatively high melt, rainfall and refreezing rates during 2020–24 compared to the surrounding years.

Fig. 8. Diagrams showing monthly melt and rain (a), refreezing (b) and runoff (c) for the RCP 4.5 scenario and averaged for the entire land area of Svalbard.

With a reduction in firn areal extent, the potential for deep water storage in perennial firn aquifers (PFAs) is reduced. PFAs have been observed in western Svalbard, on Holtedahlfonna Ice Cap (Christianson and others, 2015) and on Lomonosovfonna Ice Cap in central Svalbard (Uppsala University; unpublished data), and develop in locations where continuously temperate and deep firn prevails (Kuipers Munneke and others, 2014), and where runoff is slow. Although temperate firn conditions at 12-m depth are a common feature in the present (Fig. 9a), a considerable reduction in the extent of temperate firn zones is projected (Fig. 9b). At locations with PFAs, water tables are likely to reach closer to the surface in the near future, possibly becoming surface lakes in local depressions. The latter effect is the combined result of increased melt water supply and reduced pore volume.

Fig. 9. Spatial distribution of subsurface temperature at 12-m depth on glaciers in 2018 (a) and 2060 (b) for the RCP 4.5 scenario. The difference between the two distributions (Δ) is shown in panel (c).

Runoff and snow season length

Increased glacier melt and rainfall, and a negligible change in refreezing, lead to a strong positive trend in runoff from glacierized areas in Svalbard (Fig. 10). The average total runoff during 2019–60, including contributions from glacier and non-glacier terrain, amounts to 97 (RCP 4.5) or 109 (RCP 8.5) Gt a⁻¹, which is more than two times the amount of 46 Gt a⁻¹ during 1957–2018 (Fig. 5c). The relative contribution of runoff from non-glacierized areas drops from 31% during 1957–2018 to 20% (RCP 4.5) or 19% (RCP 8.5) during 2019–60. As for CMB, the largest changes are expected at the lower elevations on glaciers in southern Svalbard, which leads to runoff increases between the historical and future periods that in some locations exceed 2 m w.e. a⁻¹ (Fig. 10c). The seasonal signature of runoff (Fig. 8c) shows an increased duration and intensity of the runoff peak in summer, caused by the combined effects of increased melt and rainfall (Fig. 8a) and reduced summer-time refreezing (Fig. 8b). Furthermore, although winter and spring runoff is negligible in the first decades of the historical run, it becomes increasingly significant, particularly towards the end of the future simulation, with major spikes that mostly represent large winter rainfall events.

Fig. 10. Spatial distribution of runoff, averaged over the periods 1957–2018 (a) and 2019–2060 (b) for the RCP 4.5 scenario. The difference between the two periods (Δ) is shown in panel (c).

Changes in snowfall and runoff affect the duration of the snow and snow-free season in the glacier ablation zones and in non-glacier terrain. The spatial mean snow-free season length is found to increase from 37 d in 1957–2018 to 87 d (RCP 4.5) or 102 d (RCP 8.5) in 2019–60. The spatial distribution of the snow-free season length for the two periods for the RCP 4.5 scenario reveals that largest changes are projected for southeastern Svalbard with local changes exceeding 150 d, i.e. more than 5 months, while the remaining firn zones in 2060 experience a zero change (Fig. 11c). The seasonal distribution of snow cover fraction (Fig. 12), defined as the fraction of the total land area in Svalbard covered with >0.01 m w.e. of snow, which shows that earlier snow disappearance in spring and later snow onset in autumn both contribute to the increased snow-free season duration. This is different to the results of Van Pelt and others (2019), who found that snow onset came later in autumn, but that there was a negligible change in snow disappearance date. Figure 12 further shows that the minimum snow cover fraction in August nearly halves from 58 to 33% (RCP 4.5) or 31% (RCP 8.5) from 1957–2018 to 2019–60, while occasional snow-free conditions may occur locally during November to March in a future climate.

Fig. 11. Spatial distribution of the snow-free-season duration, averaged over the periods 1957–2018 (a) and 2019–60 (b) for the RCP 4.5 scenario. Differences between the two periods (Δ) are shown in panel (c).

Fig. 12. Monthly snow cover fraction for the RCP 4.5 scenario during 1957–2018 and 2019–60. Here, the snow cover fraction is defined as the fraction of the total land-area covered with >1 cm w.e. of snow. Monthly values are determined as the mean of daily snow cover fraction values.

Uncertainties

A discussion of uncertainties affecting the results for the historical run (1957–2018) was given in Van Pelt and others (2019), and included a description of uncertainty induced by the use of a fixed glacier geometry, the uncertainty of the climate forcing, and modelling errors. In the ‘Methods’ section we have described the choice for the approach to generate future climate forcing, in which advantages and drawbacks are mentioned. For more discussion on this the reader is referred to Van Pelt and others (2016b), where a similar method to construct a future meteorological forcing was first used. In the remainder of this section, we discuss the impact of using a fixed geometry for the glaciers in our simulations.

The use of a constant glacier geometry implies that potential changes in glacier area and elevation are ignored. As discussed in Van Pelt and others (2019), this means that our model output is calculated for a so-called ‘reference surface’ (Elsberg and others, 2001). A drawback of this is that, in case the actual surface and reference surface start to deviate much, that the calculated values no longer represent nature, but rather become virtual quantities. An advantage is that the effect of climate trends on glacier mass change can be isolated, i.e. that the results do not show the mixed signature of climate-induced trends and trends induced by geometric changes.

To quantify the impact of ignoring changes in glacier geometry, we performed a sensitivity experiment as described in the ‘Methods’ section. In the experiment, volume changes of all large (>20 km²) land-terminating glaciers in Svalbard are modelled over the period 2019–60 for the RCP 4.5 scenario. Glacier geometries evolve by applying annual CMB values to update surface heights at the end of every mass-balance year (31 August) according to a relation by Huss and others (2010). A comparison of total runoff and CMB with and without evolving geometry is shown in Figure 13. We find that the compensating effects of thinning-enhanced surface melting and the melt reduction following retreat of glacier fronts have a nearly balancing effect on both runoff and CMB. By the end of the simulation in 2060, when largest discrepancies prevail, runoff increased from 16.9 Gt a⁻¹ for a static surface to 17.6 Gt a⁻¹ for a dynamic surface, equivalent to a 4% increase. For comparison, with a dynamic geometry the total volume of the considered glaciers drops by as much as 24% in 2060 between the two experiments (Fig. 13b), while the glacier area shrinks by 14%. It should be noted that the sensitivity experiment above only applies to a limited sample of glaciers in Svalbard, comprising 18% of the total present-day glacier volume in Svalbard, and that potentially different sensitivities may apply for tidewater glaciers and small land-terminating glaciers. This could not be quantified here, and ideally requires coupling of the current model with an ice flow model, which is left for future study.

Fig. 13. Time series of total runoff (a) and CMB (b) for simulations of large land-terminating glaciers with a fixed and dynamic surface for the RCP 4.5 scenario. Additionally, total glacier volume is shown on the right axis in panel (b).

Conclusions

We presented projections of future CMB, snow conditions and runoff for Svalbard. With a coupled energy balance–subsurface model simulations for two emission scenarios (RCP 4.5 and RCP 8.5) were performed and output for 2019–60 is compared against previously published results in Van Pelt and others (2019) for the period 1957–2018. Meteorological input for the future simulations is constructed from a combination of long-term trends, extracted from a set of future RCM simulations, and historical short-term variability.

We find that the average CMB for Svalbard glaciers, which was weakly positive during 1957–2018, becomes negative at an accelerating rate during 2019–60 for both scenarios. A major factor contributing to this mass loss acceleration is the ELA to lie at or even above the peak elevation band in the glacier hypsometry, which makes glacier CMB more sensitive to climate trends. This effect is most strongly seen in southern Svalbard, where the model predicts the strongest ELA rise and where the hypsometry peak is most narrow. As a result of rising ELA, the AAR is expected to reach zero for the first time in the 2050s in northern Svalbard, while in southern Svalbard, zero AAR is predicted to occur already in the 2030s. This has a major impact on the state of snow and firn; by 2060 the total snow and firn pore space in the upper 12 m of the model is found to be reduced by 70% (RCP 4.5) or 80% (RCP 8.5) relative to 2018, with the only significant remaining firn area to be found on Lomonosovfonna Ice Cap in central Svalbard. We further find that although annual refreezing in seasonal snow and firn remains largely unchanged, its seasonal distribution changes markedly, with higher refreezing rates in winter and spring in connection with increasingly frequent rainfall events, and lower refreezing in summer due to a reduction in the cold content. With shrinking temperate firn zones and the aforementioned loss of firn pore volume in a future climate, the extent and volume of PFAs is expected to decline. For seasonal snow on glaciers and in non-glacier terrain, we project a more than doubling of the snow-free season length from 1957–2018 to 2019–60, with comparable contributions from earlier snow disappearance in spring and later snow arrival in autumn. Total runoff is found to increase markedly from 46 to 97 (RCP 4.5) or 109 (RCP 8.5) Gt a⁻¹, with an increasing relative contribution from glacier-covered areas. Finally, a sensitivity experiment with dynamic glacier geometry for large land-terminating glaciers was performed, which suggests that assuming constant glacier geometry only has a limited impact on the simulated runoff and CMB.

Data

The data shown in Figures 2–13 and Appendix Figures 14–17 are published in the online data repository (Van Pelt and others, 2021). A complete list of available model output variables at daily temporal and 1-km spatial resolution can be found at http://www.wardvanpelt.com/model_output.txt. These data are provided on request by contacting the main author.