Ice thickness measurements were collected by towing a Blue Systems ice-penetrating radar (Blue Systems Integration Ltd., 2013) with a low-frequency (2.5 MHz) antenna and an on-board GPS behind a snowmobile (about 20 km h⁻¹). The transmitter (Tx) and receiver (Rx) were in two different sledges in line with each other. The antenna arms were protected by flat PVC tubing. The spacing between Tx and Rx was 60 m, and each of the two antenna arm lengths was 16.8 m. Combined, the total drag after the snowmobile was about 120 m. There is a schematic figure of the radar setup in Supplementary Appendix S1.

The antenna spacing was tested in the field, and 60 m proved to be the optimal distance for a clear reading. The signal’s vertical range was set to 50 mV, and the sampling rate was set to 250 MHz. The trigger selection edge regulated the level (mV) at which the measurement should be triggered (trigger level). The trigger level was set to 40 mV, but it had to be adjusted a few times in the field to trigger the measurement. Both 5 and 10 MHz antennas were tested, but they proved insufficient for penetrating to the bedrock in large regions where the ice thickness exceeded about 100 m, probably due to the high water content in the ice (Rignot et al., 2013). Profiles measuring the transition from the ice edge to thicker ice (about 50 m) were collected using a Malå ground penetrating radar (GPR) system with a 25 MHz Rough Terrain Antenna (RTA) series antenna with a 12 channel GPS (± 5 m accuracy). Data from Sørfonna were collected for 5 days in March of each year from 2015 to 2018. The data were collected in a grid, with longitudinal spacing of 500 m (NE-SW) and transverse spacing of 1,000 m (NW-SE). Data from Nordfonna were collected by Førre (2012). Post-processing and analysis of the radar profiles were conducted using Reflexw software (Sandmeier, 2016) with standard radar processing tools (e.g., interpolation of GPS measurements to all traces, trace interpolation and removal of standstill traces, static correction, application of an AGC gain filter, equidistance traces, Stoltz migration, and time-to-depth conversion).

The precision of the ice thickness dataset was investigated by comparing ice thickness measurements at all locations where radar profiles crossed (i.e., crossover analysis). The radar GPS elevations were compared to a surface DEM from 2013 (Kartverket, 2020). The subglacial terrain map was interpolated and constructed according to the protocol described in Supplementary Appendix S2. The interpolation method ‘Topo to Raster’ was used to build DEMs (ESRI, 2016a) suitable to perform a hydrological analysis on the subglacial topography map. The interpolated subglacial terrain map was integrated with a DEM (2013; Kartverket, 2020) of the present landscape to visualize the Folgefonna peninsula as ice-free. Both the interpolated subglacial terrain map and the surface DEM have a resolution of 10 m. The hydrological toolbox in ArcMap 10.6 was used to interpolate and to perform hydrological analyses (e.g., mapping of drainage catchments and mapping of future lakes; ESRI, 2016b). The drainage catchments are delineated from the topography which includes the surface of ice cap, both at present and, in future scenarios when the ice cap still exists. Subglacial drainage and water path alterations made by the hydropower companies are not accounted for.

The ice thickness of the ice cap was mapped in two ways: first, by subtracting the subglacial terrain map from the glacial surface DEM, and second, by interpolating from the ice thickness measurements directly (with the glacial boundary set to zero). Mean ice thickness and ice volume were calculated (in Python) by integrating the ice thickness raster.

Issues related to ice thickness differences in crossing profiles were manually examined after the first interpolation of the subglacial topography map (step seven in Supplementary Appendix S2) to remove obvious errors, such as non-natural features (e.g., spirals or overlapping contours) caused both by general and specific uncertainties in radar measurements. For example, the radar system is sensitive to situations in which the antenna distance is compromised, such as turns or steep slopes, so precautions were taken when driving the snowmobile during the survey. The radar itself does not gather data at one point directly beneath the system, but within an elliptical area, called the Fresnel zone (Robin et al., 1969). Therefore, the roughness of the bed and undulating topography can contribute to the theoretical error of the ice-penetrating radar measurements. A migration filter was applied while processing the radar profiles to reduce the effect of the Fresnel zone and shift the data points to the correct coordinates (Lindbäck et al., 2014). In steep terrain, applying the migration filter works to correct the slope of dipping bedrock topography. This may lead to differences in interpreted ice thickness between crossing radar profiles, as after migration only the profiles running perpendicular to the dip direction of the subglacial topography will display the correct ice thickness.

The horizontal accuracy of mapping is limited by the GPS precision. The Blue Systems ice-penetrating radar uses two accuracies: GPS Standard Position Service, with an accuracy of less than 15 m (95%), and Wide Area Augmentation System (WAAS), with an accuracy of less than 3 m (95%). The system automatically chooses what is available. The WAAS was used for most but not all traces, and therefore the dataset’s accuracy was less than 15 m. The number of satellites used by the GPS to determine its position was between eight and eleven. Overall, the uncertainties related to the horizontal accuracy of the radar measurements and the applied measurement grid suggest that the data are good enough to capture landscape features ranging from tens to a few hundreds of meters. Therefore, we discuss landscape formations rather than individual point features.

The distance to the glacial bed was calculated using the two-way travel time and standard propagation velocity for electromagnetic waves through ice. The velocity in ice is dependent on temperature and density; therefore, it is slightly different for each glacier (Gruber and Ludwig, 1996). Because the exact velocity of electromagnetic waves travelling through the ice at Folgefonna is unknown, and because we did not have the right set of circumstances in which to measure it, we assumed homogeneous ice cover and ice density. Dry ice is prone to higher velocities (1.68–1.75 × 10⁸ ms⁻¹), and wetter ice is prone to lower velocities (1.55–1.65 × 10⁸ ms⁻¹; Magnússon et al., 2016). However, studies conducted in northwestern Iceland successfully used velocities of 1.70 × 10⁸ ms⁻¹ or even up to 1.75 × 10⁸ ms⁻¹, depending on the condition of the ice (Magnússon et al., 2016). A two-way travel time of 1.68 × 10⁸ ms⁻¹ was chosen for Sørfonna based on previous ice-penetrating radar studies conducted in Greenland and Canada (Wilson et al., 2013; Lindbäck et al., 2014; Yde et al., 2014) and the previous study at Nordfonna (Førre, 2012).

The crossover analysis was performed to check the consistency of the ice thickness measurements and to estimate vertical uncertainties in ice thickness. The crossover analysis identified 275 crossover points, as seen in Figure 2, and an average ice thickness difference of 20 m, with values ranging between 0 and 152 m (see distribution in Figure 2). Minor differences are to be expected due to different measurement years and horizontal GPS uncertainty, and most differences are below 50 m (Figure 3). Only one point, located at the deepest part of the glacier in steep terrain, has a difference of more than 100 m. Here the profile along the valley measures 404 m in ice thickness, while the radar profile perpendicular to the valley measures a 556 m ice thickness (152 m difference). This large difference is most likely a result of different profile orientations, as only dipping reflectors in the profile perpendicular to the subglacial topography will be corrected by the migration routine. The interpolation is based on the full dataset, but step seven in the interpolation method (Supplementary Appendix S2) was used to clear these kinds of differences. Nevertheless, manual work is not always fully reliable, even though precautions were taken to be as thorough as possible.

As the elevation accuracy of the radar GPS is difficult to quantify, we compared these measurements to the DEM from 2013 (Kartverket, 2020). The GPS elevation measurements for each radar profile were subtracted from the DEM (Figure 3). The maximum difference is ±30 m, but the difference between most values is within 10 m, and the mean is 5.5 m. There are unavoidable differences between the GPS elevation measurements and the DEM, as the DEM is from 2013 (Kartverket, 2020), and the ice cap has been melting over the past few years (Nesje et al., 2008; Andreassen et al., 2012). As discussed earlier in this section, a 10 m accuracy is within the uncertainty of radar measurements (see the crossover analysis in Figure 2); therefore, both the DEM and GPS data will provide a similar result when used for interpolation and calculation of ice thickness. At Sørfonna, an accuracy of ±10 m is important at its edges, where the ice cap is thinner; in the middle part, where the ice thickness is between 300 and 600 m, 10 m will not have too much effect on the overall shape of the glacier and the glacial bed. For this reason, we analyse large-scale bedrock features and avoid discussing exact numbers for ice thickness when they are smaller than tens of meters. We also cannot discuss individual outlet glaciers’ specific ice thicknesses, as our dataset does not cover most of the outlet glaciers.

The error sources described in this section may have affected the accuracy of the estimated mean ice thickness and volume. To test this, we compared the ice thickness map made in this study to measurements obtained by NVE and Statkraft using radar between 1964 and 2004 on the northwestern part of Sørfonna (Brekke and Nord, 2008). The comparison shows a similar ice thickness pattern in this area (see Supplementary Appendix S3). Of course, the exact numbers differ somewhat because the ice cap and outlet glaciers have both advanced and retreated in past decades, and on the previous map there are more measurement points on the outlet glaciers but fewer on the actual main body of the glacier.

The model is based on the continuation equation for mass conservation and is a combination of the two-dimensional shallow ice flow model developed by Meur and Vincent (2003) and the one-dimensional model described by Oerlemans (2001). The model requires five input parameters: glacial bed elevation, glacier surface elevation, mass balance data, quantification of internal deformation, and subglacial sliding. In this study, the mapped subglacial bed, a glacier surface (Kartverket, 2020), and mass balance data (Andreassen et al., 2012), which are described further in the next section, were used as input data. Internal deformation and subglacial sliding were calculated based on empirical parameters following the method described by Laumann and Nesje (2014). The model has previously been used on the plateau glacier Spørteggbreen in western Norway (Laumann and Nesje, 2014) and on Hardangerjøkulen (Laumann and Nesje, 2017). The model is derived from the model developed by Giesen (2009) which has for example been used on Hardangerjøkulen by Giesen and Oerlemans (2010).

To calculate how climate will affect the ice cap, temperature and precipitation must be transformed into mass balance. A simple degree-day model was used for this purpose, with daily temperature and precipitation in Bergen (65 km west of Folgefonna; Smith-Meyer and Tvede, 1996; Braithwaite, 2011; Figure 1) used as input parameters. Bergen was chosen since meteorological stations located in the immediate proximity to Sørfonna would not be applicable for the whole ice cap due to significant difference in local climate close to the ice cap. The model calculates snow accumulation and melting at different elevations on the ice cap. To calibrate the model, the temperature and precipitation gradients (grt, grP), the melting factors for snow and ice (cs, ci), and a correction factor (corr.) for the distance between Bergen and Sørfonna were used to tune the data. This was done by integrating actual mass balance measurements taken during field surveys (Kjøllmoen, 2014; Kjøllmoen, 2018; Norwegian Water Resources and Energy Directorate, 2019) of the outlet glacier Svelgjabreen (southern part of Sørfonna), for which there are 11 years of continuous measurements (2007–2017), and the outlet glacier Gråfjellsbrea (northern part of Sørfonna), for which there are 10 years of continuous measurements (2003–2012).

Tuning each of the above-mentioned parameters is also partly dependent on the tuning of each parameter separately. However, using mass balance data, it is possible to independently calculate the grt for four of the 11 years. The annual measured mass balance for the outlet glacier Svelgjabreen was positive for the three highest elevations. This means that in those years, only snow (not ice) melted at these locations. On average, the ratio between measured melting at the lowest elevation (1,300 masl.) and at the highest elevation (1,600 masl.) is 1.55 m w.e. Assuming a constant melting factor for snow, cs, this will also be the ratio between degree-days at the same elevations. To find this number, grt is the only unknown, and it is calculated to be –0.62°C (100 m) ⁻¹. With known grt, grP is calculated using measured accumulations over 11 years. The two remaining parameters (cs and ci) are found by comparing the average measured and calculated summer balances. The parameters used for Svelgjabreen in the model are as follows: grt = –0.62 ^oC per 100 m, grP = 8% per 100 m, corr. = 1.00, cs = 0.0041 m w.e. ^oC⁻¹d⁻¹, and ci = 0.0059 m w.e. ^oC⁻¹d⁻¹. The parameters used for Gråfjellsbrea are estimated in the same way as for Svegjabreen, and they are as follows: grt = –0.6°C per 100 m, grP = 9% per 100 m, corr. = 0.65, cs = 0.0044 m w.e. ^oC⁻¹d⁻¹, and ci = 0.0060 m w.e. ^oC⁻¹d⁻¹.

Mass balance measurements from NVE show that the mass balance of Sørfonna has a large gradient from north to south (Østrem et al., 1988). The mass balance at every grid point on Sørfonna is thus calculated as a function of elevation and have a latitudinal “lapse rate”.

The model can be tested by modelling a documented time period that has already passed. The glacial extent of Sørfonna was mapped both in 1959 and 2013 (Smith-Meyer and Tvede, 1996). Therefore, the model was tested using the glacial extent in 1959 as input for the model and running it annually until 2013 using the existing calculated and measured mass balance. A comparison was then made between the modelled output and the measured glacial extent in 2013 to verify the model’s robustness (Figure 4).

The radar profiles (Figure 5) used in the interpolation process of Sørfonna cover the ice cap in a grid with line spacing of approximately 500 × 1,000 m, as seen in Figure 5A. The middle part, where we began each day in the field (from Fonnabu; Figure 1), features denser line spacing. Some lines are discontinuous in the southern part, making this area slightly less surveyed. As the example profile illustrates (Figures 5B,C), the subglacial topography is far from flat and fluctuates from valleys to mountain peaks. The measurements show that the ice cap is thickest in the south-central part, reaching about 200–600 m, and thins towards the edges as the outlet glaciers flow down into valleys (Figure 6). The northern part of the ice cap is generally thinner, similar to Nordfonna, but it still reaches an ice thickness of more than 300 m in the central parts (Figure 6).

The results of the two ice thickness maps (interpolated versus subtracted from DEM) were similar (mean difference = –3 m, standard deviation = 22 m; see figure in Supplementary Appendix S4). The ice thickness map presented in Figure 6 is derived using the DEM surface (see Section 3.1). The maximum measured ice thickness is 568 ± 20 m, the mean ice thickness of the ice cap is 188 ± 20 m, and the volume of the ice cap is 28 km³. Uncertainties are reviewed further in the discussion section.

The interpolated bed topography of both Nordfonna and Sørfonna shows that the smooth dome-shaped glacier surface covers a landscape with deep valleys, hanging valleys, cirques, and glacial over-deepenings, similar to the landscape surrounding the ice cap (Figure 1). The subglacial landscape of Sørfonna has two main valleys draining towards the southwest, with valley floors at 700–800 masl. North, east, and west of these valleys are subglacial mountain peak areas with peaks between 1,200 and 1,500 masl. (Figure 7), resulting in ice thicknesses of less than 50 m in the central parts of the southern dome (Figure 6). There are several high-elevation areas (> 1,400 masl.) in the northern part of the subglacial landscape that are incised by valleys (800–900 masl.) oriented towards the southeast and the northwest (Figure 7). Some valleys and over-deepenings in the landscape are naturally surrounded by water divides and thus may be the locations of lakes in the future. When the ice cap has completely melted and the landscape is ice free, the hydrological analysis suggests that up to thirty new lakes of different sizes will be formed. The largest lake will have an area of 4.7 km² and a volume of 0.28 km³ (Figure 7; Table 2). A future ice-free landscape, including future lakes, of both Sørfonna and Nordfonna is presented in Figure 7. The subglacial landscape of Nordfonna is based on the dataset from Førre (2012) but interpolated using the approach described in Supplementary Appendix S2.

Figure 8 shows the modelled Sørfonna ice cap for the years 2120 and 2170, as well as drainage catchments and lakes. Under climate scenario 1, a large part of the 164 km² ice cap is still present in 2120 (101 km²) and in 2170 (76 km²). Under scenario 2, the ice cap is very small (31 km²) in 2120 and likely to be completely melted by 2170 (Table 3).

Using the 1959 extent as input for the model (Figure 4), we successfully reproduced the glacial extent in 2013. Some of the outlet glaciers are not at their full size, and some are larger compared to the DEM measured in 2013 (Kartverket, 2020), but the ice body and the elevation contours of the ice cap concur with the measured and calculated glacial size from 2013. Running the model for the two climate scenarios showed that the southern part of Sørfonna is relatively stable, while the northern part of Sørfonna melts away more quickly (Figure 8).

The present drainage system of the landscape where Sørfonna lies is divided between three major drainage catchments and seven smaller ones (Figure 1). These shift in size when the ice cap melts and therefore altering drainage divides defined by the ice cap. In both climate scenarios, the southwestern drainage catchment, which includes the two major subglacial valleys in the southern part of the glacier, will increase in size, gaining area from the neighbouring catchments (Figure 9). In climate scenario 1, the southwestern drainage catchment will increase by 8 km² by 2120 and 30 km² by 2170. In climate scenario 2, the southwestern catchment will increase by 36 km² by 2120 and 59 km² by 2170. In the northern part of Sørfonna, the drainage configuration is more stable and will remain more or less unchanged as the ice cap shrinks. This change in catchment size indicates a shift towards a more southerly drainage in this landscape compared to present conditions. For a visualization of changes in the drainage catchments and glacial extent for every tenth year, see Supplementary Appendix S5.

The ice thickness map of Sørfonna (Figure 6) shows a maximum ice thickness of above 550 m (the measured maximum depth is ∼570 m), which significantly exceeds the previously measured maximum of 375 m and the estimated maximum of 400 m (Brekke and Nord, 2008; Norwegian Water Resources and Energy Directorate, 2019). Based on our dataset the mean ice thickness is close to 190 m, in contrast with the previous estimate of 140 m (Kennett and Sætrang, 1987; Andreassen et al., 2015). A mean ice thickness of ∼190 m is similar to other ice caps in Norway. For example, the smaller ice cap Hardangerjøkulen (71 km² in area), ∼70 km northeast of Sørfonna, has a mean ice thickness of 150 m and a maximum depth of 380 m (Andreassen et al., 2012). The larger ice cap Jostedalsbreen (474 km² in area), ∼160 km north of Folgefonna, has a mean ice thickness of 160 m and a maximum depth approaching 600 m (Andreassen et al., 2015, and references therein). The previous volume estimate for Sørfonna was 7 km³ (Andreassen et al., 2015), which is lower than Hardangerjøkulen (11 km³) and Jostedalsbreen (49 km³). The calculated volume of close to 30 km³ for Sørfonna presented here is a reasonable volume for an ice cap with the observed area and ice thickness.

This study, combined with the results from Nordfonna (Førre, 2012), shows that the Folgefonna ice cap covers a complex landscape, with deep valleys, hanging valleys, cirques, mountain ridges, and over-deepenings. Qualitative analyses of the subglacial topography below Sørfonna reveal a general consistency with Nordfonna subglacial topography and the landscape surrounding the ice caps. The main valleys of the Sørfonna subglacial topography follow a southwest-northeast direction, lining up with structural weakness zones in the Hardangerfjorden area (Fossen and Hurich, 2005).

The extent and thickness of ice caps controls the routing of meltwater into surrounding valleys and proglacial and distal glacier-fed lakes. It is therefore possible to constrain bedrock thresholds identified on the subglacial topography map based on numerous existing records of lake sediments surrounding Sørfonna (Figure 10). On the west side of the ice cap, in the valley Bondhusdalen (below Bondhusbreen; Figures 1, 10), radiocarbon ages were reported for the distal glacier-fed Lake Bondhusvatnet (Figure 10) based on a 2.5 m long sediment core (Simonsen, 1999). The uppermost 2 m of the core contains inorganic, partly laminated sediments with low organic content, indicating that the sedimentation is dominated by silt and clay eroded by the ice cap. The lowermost 50 cm of the core contains a high level of organic sediments with high loss-on-ignition values (> 70%). This stratigraphy indicates that there was no input of glacier meltwater into the lake prior to 2240 ± 40 cal yr BP, and that this marks the time when Lake Bondhusvatnet first became a lake in response to the increased input of meltwater from the ice cap (Simonsen, 1999). This is later than the first regrowth of Sørfonna after the Holocene Thermal Maximum (HTM), at ca. 6,900 cal yr BP, which was recorded in Lake Buervatnet (Figure 10) on the eastern side of Sørfonna (Røthe et al., 2019). The subglacial topography map of Sørfonna confirms a threshold below Bondhusbreen, which potentially will redirect its meltwater from the high mountains, from north to south, when the glacier has retreated behind the threshold. Exactly when this will happen can be determined from denser measurements of ice thickness along the threshold, followed by more accurate ice flow modelling in this area, combined with down-scaled climate predictions based on earth system modelling.

On the east side of Sørfonna, two outlet glaciers (Øvre and Nedre Buerbreen; Figures 1, 10) drain towards the distal glacier-fed Lake Sandvinvatnet (Figure 10). The mapped subglacial topography shows a bedrock threshold in the upper part of Buerbreen with an over-deepening to the southwest, which will likely be filled with water when the glacier retreats (Figure 10). Analysis of the sedimentology in Sandvinvatnet indicates that the first glacially derived sediments entered the lake around 4,800 cal yr BP (Ekblom Johansson et al., 2020). Before 4,800 cal yr BP, only organic sediments were present, indicating no glacier in this sub-catchment. Three short recurring periods (2200–1900, 4,600–4,400, and 4,800–4,100 cal yr BP) also show only organic material, which strengthens the theory of a threshold below the outlet glaciers. This is confirmed by observations in the field, where exposed bedrock is visible in the icefall (see image in Supplementary Appendix S6). The glacial return at 4,800 cal yr BP occurred later than for many other glaciers in Norway and differs both from the returns of similar dates at Nordfonna (Bakke et al., 2005) and from glacial reconstruction from the nearby Lake Buervatnet (Røthe et al., 2019). Naturally, glacier growth in the region currently covered by Folgefonna would have begun as isolated patches of ice in high-elevation areas that in time grew to become larger and lower glaciers that eventually merged into the ice cap we see today. Therefore, different glacial return dates after the HTM are possible around the boundary of the present ice cap.

In the southwestern corner of Sørfonna, there is an outlet glacier (Møsevassbreen; Figures 1, 10) with a calving terminus leading into Lake Møsevatn (Figure 10). Reconstruction of former ELA variations based on a study of the sediments deposited in this lake shows a remarkable drop in ELA from 1446 to 1300 cal yr BP, which led to considerable glacier growth (Bjønnes, 2006). After the readvancement of Møsevassbreen, the glacier has been continuously present in the lake catchment. The first part of the record, which covers the past 7,000 years, indicates that no glacier meltwater drained into the lake. This confirms our interpretation of the subglacial topography map; we argue that there is a threshold in the upper part of Møsevassbreen that will direct meltwater from the mountains in an easterly direction when the glacier is smaller than it is today. Inorganic sediment layers from the lake could indicate smaller glacier events from 6,612 to 5,791 cal yr BP, 4,722–4,181 cal yr BP, and 3,000–2900 cal yr BP (Bjønnes, 2006).

A similar study (Wittmeier, 2010) was conducted down-valley from the outlet glacier Svelgjabreen (also called western Blomsterskardsbreen in some studies; Figures 1, 10) based on sediment cores from Lake Midtbotnvatnet (Figure 10). By studying the sediments deposited in this lake, the authors concluded that glacial clay and silt had been continuously deposited from the earliest date (6,295 ± 115 cal yr BP) until the present (Wittmeier, 2010). Their results confirm that the valley below Svelgjabreen drains meltwater from the thickest part of the Sørfonna ice cap. However, it is not possible to conclude whether the glacier was melted completely during the HTM. This contrasts with the work of Nesje et al. (2008) and Bakke et al. (2005), who conclude that all studied glaciers in Norway melted away during the early- and mid-Holocene and regrew during the mid- and late-Holocene.

By simulating the evolution of Hardangerjøkulen (the closest other ice cap to Folgefonna) from the mid-Holocene to the present, Åkesson et al. (2017) found that Hardangerjøkulen is highly sensitive to climate change. With a projected future warming of 3°C, the model indicates that the Hardangerjøkulen ice cap will disappear almost entirely before the end of the 21^st century, and that not even an improbable 100% increase in precipitation can compensate for the effect of the warming (Giesen and Oerlemans, 2010).

In southern Norway, other glacier model simulations have been carried out using one-dimensional models of Spørteggbreen (Laumann and Nesje, 2014); Nigardsbreen, an eastern outlet glacier of the Jostedalsbreen ice cap (Oerlemans, 1997); Briksdalsbreen, a western outlet glacier of the Jostedalsbreen ice cap (Laumann and Nesje, 2009a; Laumann and Nesje, 2009b); and Storbreen glacier in Jotunheimen (Andreassen and Oerlemans, 2009). Similar to the results from Hardangerjøkulen, sensitivity experiments at Spørteggbreen for a model run from 2011 to 2100 show that an increase in winter precipitation of about 70% is needed to compensate for a warming of 2.3°C (Laumann and Nesje, 2014). The results from Sørfonna presented here show that, at least in the central and southern parts, the ice cap is more resistant to climate change than Nordfonna and other glaciers in western Norway. Runs of the scenario 2 model in our study indicate some ice remnants in 2120, and in scenario 1, there is still ice at Sørfonna in 2170 (Figure 8). Given a present ice volume of 28 km³ (Table 3), Sørfonna has an estimated loss of 43% (17 km³) and 92% (2.5 km³) in climate scenarios 1 and 2, respectively.

Simulations conducted with the Global Glacier Evolution Model (GloGEM)—forced by RCP 2.6 (similar to our climate scenario 1), RCP 4.5, and RCP 8.5 (similar to our climate scenario 2) emission scenarios–yield a global glacier volume loss of between 25 and 48% from 2010 to 2100 (Huss and Hock, 2015). The RCP 4.5 scenario indicates glacier volume losses of 20% in the Arctic, Antarctic, and Subantarctic regions, and 90% in low-latitude regions such as Central Europe (Huss and Hock, 2015). When the evolution of twenty-nine glaciers on Svalbard is modelled under the RCP 8.5 scenario and extrapolated to all land-terminating glaciers on Svalbard, Möller et al. (2016) projections suggest that 1435 ± 59 of the 1471 land-terminating glaciers on Svalbard will be lost by 2100, indicating that it is possible that all the glaciers could disappear. The projected overall volume loss, however, only ranges between 55 ± 36% and 58 ± 37%, due to the nonlinear relation between glaciated areas and ice volumes (Möller et al., 2016). Other studies indicate larger volume losses of the Svalbard glaciers by 2100. For example, Marzeion et al. (2012) project a reduction of 70–84% using RCP 8.5 forcing, and Radić et al. (2014) estimate a reduction of 79–95% following an RCP 8.5 scenario to 2100.

In sum, these results indicate that glaciers with a hypsometry containing small areas above the ELA, such as Spørteggbreen, and thinner ice caps, such as Hardangerjøkulen, are more sensitive to future climate change than Sørfonna. Both the thick ice and the strong precipitation gradient, with more winter precipitation in the south (about 5,000 mm) compared to north (about 3,350 mm), may explain the robustness of the southern part of Sørfonna compared to its thinner and more climate-sensitive northern part. The high amount of winter precipitation may also explain why Sørfonna is likely to evolve more similarly to glaciers at higher latitudes, such as those on Svalbard, and to show more resilience to future warming than other glaciers in southern Norway and Central Europe.

This study presents the results of an ice flow model based on two scenarios, both of which involve linear forcing and anticipate a gradual rise in temperatures. Despite their simplicity, the model scenarios provide a decent foundation for strategies to adapt to a future with no ice present in the landscape, which are needed by both hydropower companies and society, dependent on the infrastructure around the ice cap. However, every model is a simplification of reality, and the results must be discussed with care. The actual retreat rate depends on the climate outcome which will likely be between the two climate scenarios. Nonetheless, the model provides a strong indication of how and when the ice cap retreats over the area.

We suggest that the resilience of the Sørfonna ice cap is due to the increase in precipitation predicted in the near future, which initially compensates for the increase in temperature. Using climate scenario 2, the ice cap is modelled to be even larger in 2070 compared to the same year using climate scenario 1 (Figure 8). When the ice cap shrinks below its ELA, melting rapidly increases because the accumulation area decreases. When the glacier is still present, drainage catchments will keep draining towards both the north and the south. However, in model scenario 2, the ice cap is completely melted by 2170 (Figure 8D), and the drainage is likely to flow southward. It is possible that the ice cap will become dynamically inactive towards the end. Already in 2120, in climate scenario 2, the ice cap might be too thin (20–100 m; Figure 8C) to be considered an active glacier. In this case, the remaining ice could still act as a water divide.

Although the coverage of the ice-penetrating radar dataset is good for the central parts of the glacier, the outlet glaciers are not as well covered. One example is the outlet glacier Bondhusbreen, where Statkraft AS gathers meltwater directly beneath the glacier at 900 masl. (Figure 1). Our subglacial topography map suggests a ridge underneath Bondhusbreen that acts like a water divide and a boundary for two future lakes. The ice thickness dataset from 1964 to 2004 (Brekke and Nord, 2008), as well as the lake sediment data from Lake Bondhusvatnet (Simonsen 1999), also indicate that there is a bedrock threshold in the upper part of Bondhusbreen. However, to be certain about the elevation of the threshold, more data are needed, with a minimum of one profile oriented perpendicularly over the ridge. The height of this ridge is important, as it determines whether the drainage of the large future lake (no. 1 in Figure 7) just south of the ridge will flow north or south, a distinction of major interest to the hydropower company. If there is ice covering the ridge, the ice flow will be oriented towards the north, allowing Statkraft to continue to intake water directly from Bondhusbreen. When the ice cap retreats beyond the threshold, according to our subglacial topography map of Sørfonna, the drainage will flow to the south. Future adaptation strategies might therefore include more tunnels–if the national park allows their construction. For Sunnhordaland Kraftlag AS, which generates hydropower in the southern sector of the glacier, the water supply seems to be less affected, especially because the valley beneath Svelgjabreen drains a huge area of Sørfonna (Figure 1). However, based on the interpreted threshold north of Møsevassbreen, it is likely that the meltwater drainage will stop in this catchment when the glacier retreats beyond the threshold (Figure 10).

For the tourist industry the area around Folgefonna will change as the glacier slowly retreats. Most of the outlet glaciers will disappear within the next century, and only the highest-lying parts of ice will survive. Based on the above discussion of thresholds under the ice, it is likely that Øvre Buerbreen, Nedre Buerbreen, Bondhusbreen, and Møsevassbreen will be the first outlet glaciers to disappear, along with the northern part of the ice cap. Buerbreen and Bondhusbreen are popular tourist attractions, particularly important to the local economy (e.g., food, accommodations, and guide services).

When Folgefonna retreats, it will expose the area around the glacier to new geohazard threats, such as floods and rock falls. The unstable period of landscape readjustment from glacial to nonglacial conditions is known as the paraglacial period and may last from decades to millennia (Ballantyne, 2002). Even in a small, rapidly responding cirque glacier system, Bakke et al. (2009) showed that the sedimentary input into proglacial Lake Kråkenes was still detectable after 90 years. In larger systems such as Folgefonna, the paraglacial period is deemed to last for millennia, with an exponential shift from hazardous GLOFs and unstable slopes to more subtle changes in glacially conditioned sediment availability (Ballantyne, 2002).

Action needs to be taken to monitor geohazards in areas with rapidly changing glaciers. Known areas with ice-dammed lakes near Sørfonna are mapped by NVE and shown in Figure 11. In the future, when a warmer climate leads to a rapid retreat of the outlet glaciers at Sørfonna, GLOFs will represent a potential danger to infrastructure and society, and the potential locations of GLOFs will shift with the glacier as it retreats. The glacial extent in 2070, shown in Figure 11, would not allow for GLOFs at Buerbreen, as the glacier is beyond the threshold. Just west of Bondhusbreen, however, where historical observations confirm that GLOFs occurred previously (Figure 11), the retreat of the ice cap will form small glacier-dammed lakes in both scenarios 1 and 2, which may cause potential hazards for the heavily touristed areas down-valley in Bondhusdalen. The retreat of Møsevassbreen (Figure 11) might also cause GLOFs when the glacier retreats up the mountain and crosses two new lakes, potentially creating unstable glacier damming. The retreat of Svelgjabreen may create GLOFs in three ways: glacier damming by the main or lateral valley due to asynchronous retreat, ice fall from calving fronts in some of the new lakes, and ice damming of lateral valleys caused by avalanches from steep slopes in the main local valley. The model experiment conducted in this study does not allow for detailed analyses of every outlet glacier, but based on our subglacial topography map, some areas need more attention and, possibly, regular monitoring in the future. See Figure 10 for the subglacial terrain of Sørfonna together with Sørfonna’s current glacial extent (2013) and the outlet glaciers discussed in this section.