Spatio-temporal variability in geometry and geodetic mass balance of Jostedalsbreen ice cap, Norway

Abstract

The Jostedalsbreen ice cap is mainland Europe's largest ice cap and accommodates 20% (458 km² in 2019) of the total glacier area of mainland Norway. Jostedalsbreen and its meltwater contribute to global sea-level rise and to local water management, hydropower and tourism economies and livelihoods. In this study, we construct a digital terrain model (DTM) of the ice cap from 1966 aerial photographs, which by comparing to an airborne LiDAR DTM from 2020, we compute changes in surface elevation and geodetic mass balances. The area mapped in both surveys cover about 3/4 of the ice cap area and 49 of 82 glaciers. The measured glacier area has decreased from 363.4 km² in 1966 to 332.9 km² in 2019, i.e. a change of −30 km² or −8.4% (−0.16% a⁻¹), which is in line with the percentage reduction in area for Jostedalsbreen as a whole. The mean geodetic mass balance over the 49 glaciers was −0.15 ± 0.01 m w.e. a⁻¹, however, large variability is evident between glaciers, e.g. Nigardsbreen (−0.05 m w.e. a⁻¹), Austdalsbreen (−0.28 m w.e. a⁻¹) and Tunsbergdalsbreen (−0.36 m w.e. a⁻¹) confirming differences also found by the glaciological records for Nigardsbreen and Austdalsbreen.

Introduction

Glaciers and ice caps worldwide have in general been retreating and losing mass (Zemp and others, 2019; Hanna and others, 2020), but the changes are far from uniform and there have been occasional periods when glaciers have advanced due to changes in atmospheric circulation or volcanic activity (e.g. Solomina and others, 2016; Mackintosh and others, 2017). Mainland Norway is one region where several periods of glacier growth have taken place around 1910, around 1930, in the 1970s and in the 1990s (Andreassen and others, 2005). Advances of glaciers in Norway have been attributed to increased winter precipitation or periods with colder temperatures (Andreassen and others, 2005; Nesje and others, 2008). The latest advance culminated around the year 2000 and since then, glaciers in Norway have had a pronounced retreat of glacier termini (Andreassen and others, 2020; Kjøllmoen and others, 2022).

Understanding spatio-temporal variability in glacier changes requires multi-temporal and spatially distributed datasets, which unfortunately are exceptionally resource-consuming to acquire on the ground. The glaciological method measures the mass balance at point locations and data are interpolated over the entire glacier surface to obtain glacier-wide averages (e.g. Cogley and others, 2011). Point-based glaciological mass-balance records are only available worldwide for selected glaciers (WGMS, 2021). The geodetic method measures the cumulative mass balance for a period by differencing digital terrain models (DTMs) and converting the volume to mass change by assuming a density. The glaciological method measures the surface mass balance, whereas the geodetic method measures the sum of surface, internal and basal mass balances. Several studies have shown that neglecting the internal and basal balance can lead to significant biases for temperate glaciers (Oerlemans, 2013, Andreassen and others, 2016, Jóhannesson and others, 2020).

While traditional photogrammetric processing requires significant manual operator input in aligning and georeferencing imagery, advances in automated photogrammetric pipelines such as structure from motion have made it possible to analyse historical and contemporary satellite and aerial datasets at extensive temporal and spatial scales (Mölg and Bolch, 2017). This has increased the availability of geodetic mass-balance records allowing glacier changes to be assessed in data-scarce regions, as well as analyses of glacier changes at continental to global scales (Hugonnet and others, 2021; Thompson and others, 2021; Berthier and others, 2023). Norway has extensive archives of historical aerial photography, which date back to the 1950s and 1960s in many regions, but most of that archive has not been utilised for assessing glacier changes. A combination of multiple data sources and methods can be used to estimate mass balances dating back to the Little Ice Age (LIA; e.g. Carrivick and others, 2019, 2020; Aðalgeirsdóttir and others, 2020; Lee and others, 2021).

The aim of this work is to present geometric and geodetic changes of the Jostedalsbreen ice cap between 1966 and 2020 based on large-scale photogrammetric processing of historical aerial photographs and airborne LiDAR surveys. To complement the geodetic estimate, we present seasonal corrections from the mapping dates until the end of the melt season for individual glaciers. We also present estimates of internal and basal ablation from dissipation of energy due to flow of ice and water over the ice cap from 1966 to 2020, including variations between outlet glaciers. We put elevation changes into a longer-term context by co-analysing glacier outline datasets from LIA, 1966, 2006 and 2019 to compute area and length changes.

Study area

Jostedalsbreen is the largest ice cap in mainland Europe and currently (2019) covers 458 km², thereby comprising about 20% of the total glacier area of mainland Norway (Andreassen and others, 2022). Jostedalsbreen has been divided into more than 80 glaciers, many of them with individual names (Andreassen and others, 2012). The ice cap had its maximum LIA extent between 1740 and 1860 and a median date of 1755 (Bickerton and Matthews, 1993; Nesje and Dahl, 1993; 2003; Nussbaumer and others, 2011; Winkler, 2021; Carrivick and others, 2022) (Fig. 1). At this time, it had an area of 568 km² and an estimated ice volume of between 61 and 91 km³ (Carrivick and others, 2022). Between the LIA and 2006, the major outlet glaciers had in combination lost at least 93 km² or 16% of their area and 14 km³ or 18% of LIA volume (Carrivick and others, 2022). The glacier area has been further reduced by 3% from 2006 to 2019 (Andreassen and others, 2022).

Fig. 1. Location map of Jostedalsbreen. The inset map shows the location of the ice cap in southern Norway. The glacier extents in 2019, 2006, 1966 and Little Ice Age (LIA) are shown for Jostedalsbreen. For surrounding glaciers, only the 2019 extent is shown. LIA new are new outlines in this study. Background mountain shadow is from the 100 m national DTM. Glacier ID from Andreassen and Winsvold (2012). Coordinate system geographical coordinates on inset and UTM 33N, datum ETRS_1989 on main map.

Continuous records of field-measured glaciological mass balance exist on two outlet glaciers, Nigardsbreen and Austdalsbreen, since 1962 and 1988, respectively, and continue today. Austdalsbreen calves into a hydro-power reservoir, and consequently the glacier is influenced by the lake-level regulations (Laumann and Wold, 1993; Andreassen and Elvehøy, 2021). Calving is accounted for in the glaciological records (e.g. Andreassen and others, 2016; Kjøllmoen and others, 2022). Additionally, glaciological mass balance has been measured at Supphellebreen, Tunsbergdalsbreen and Vesledalsbreen for shorter periods (Andreassen and others, 2005; Kjøllmoen and others, 2021, Table 1).

Table 1. Overview of periods and number of years (n) of glaciological mass balance and front variation measurements at outlet glaciers from Jostedalsbreen up to and including 2020

Glacier ID refers to Andreassen and others (2012). Data from Kjøllmoen and others (2022). See Figure 1 for location of glaciers.

Jostedalsbreen has a maritime climate with high precipitation rates and mild winters. The glacier mass-balance records show high mass turnover with on average 8 m w.e. ablation at the glacier tongue of Nigardsbreen (600 m a.s.l.) and snow accumulations of 7 m at the plateau. The glaciological mass-balance series have been reanalysed and categorised as ‘original’ (as published in Glasiologiske undersøkelser i Norge/Glaciological investigations in Norway (e.g. Kjøllmoen and others, 2022)), ‘homogenised’ (recalculated using homogenised methods and datasets) or ‘calibrated’ (periods are calibrated with geodetic observations) according to their reanalysis status (Andreassen and others, 2016). For Nigardsbreen, the series were calibrated for the periods 1984–2013 and 2013–2020 due to significant deviation between glaciological and geodetic mass balance (Andreassen and others, 2016; Kjøllmoen, 2016, 2022) following Zemp and others (2013). Part of the deviation was attributed to internal ablation (Andreassen and others, 2016).

Glacier front variation (length change) is presently measured at seven outlet glaciers: Austerdalsbreen, Brenndalsbreen, Fåbergstølsbreen, Nigardsbreen, Stigaholtbreen, Tuftebreen and Vetle Supphellebreen, and has been previously measured at nine other glaciers (Table 1). Several of the front variation series have been terminated in recent years because glacier termini have retreated such that accurate or meaningful measurements are no longer possible (Andreassen and others, 2020; Andreassen and Elvehøy, 2021). The series of Briksdalsbreen, Fåbergstølsbreen and Stigaholtbreen have been revised due to differences between length changes measured on maps and the field observations (Kjøllmoen and others, 2007, 2019, 2020; Andreassen and Elvehøy, 2021). Advances of Jostedalsbreen outlet glaciers were recorded around 1910, around 1930, around 1950, in the second half of the 1970s and in the 1990s (Andreassen and others, 2005; 2020; Nesje and others, 2008; Winkler and others, 2009).

Ice thickness measurements have been carried out in several field campaigns using hot steam drilling at selected points and ground penetrating radar (e.g. Østrem and others, 1976; Sætrang and Wold, 1986) revealing ice thicknesses of up to ~600 m (Andreassen and others, 2015).

Many of the Jostedalsbreen outlet glaciers showed advances that culminated around 2000. Since then, glaciers have reduced in both length and mass (Andreassen and others, 2020). Geodetic mass-balance results for three outlet glaciers, which together cover 20% of the total Jostedalsbreen area, have revealed large differences in surface elevation change and geodetic mass balance. Whereas Nigardsbreen had a small mean surface elevation change of −2.2 m from 1964 to 2013, Tunsbergdalsbreen had a mean surface lowering of 40 m in the same period (Andreassen and others, 2020). The surface elevation change of Austdalsbreen was −17.4 m for the shorter period 1966–2009 (Andreassen and others, 2020). Therefore, to improve the spatial and temporal coverage and spatial resolution of glacier surface elevation changes and so to analyse controls on glacier mass loss, this study presents the results of geodetic mass balance for a larger part of the ice cap from 1966 to 2020.

Data and methods

Glacier outlines

Little Ice Age (LIA) outline

LIA glacier outlines for the most prominent outlet glaciers of Jostedalsbreen were obtained from Carrivick and others (2022), who mapped the former glacier extent based on geomorphological evidence such as moraine ridges and trimlines. For glaciers without such evidence, the 2006 outlines from the then most recent inventory were kept. The mapping was mindful of sites where moraines had been dated by lichenometry and using some historical records and used 1 m resolution LiDAR-derived topography data and sub-metre resolution optical satellite images. Uncertainty in the LIA outline positions was assessed by Carrivick and others (2022) to be <10% for the total area of Jostedalsbreen, and substantially better than this in valleys with dated moraines, owing to the high resolution and precision of these datasets, the presence of dated moraines in many valleys and the (mostly very high) clarity of the geomorphological evidence.

The LIA outlines by Carrivick and others (2022) do not have a mapped LIA maximum extent for Austdalsbreen and Sygneskarsbreen due to the hydropower lakes obscuring the geomorphological evidence (Fig. 2). In this study, we were able to identify LIA outlines for Austdalsbreen, Sygneskarsbreen and glacier 2485 with recourse to the 1966 datasets (Figs 1, 2). Specifically, we digitised the LIA outline for the three glaciers based on trimlines identified in orthophotos from 1966 (see next section) and other images available on norgeibilder.no and as observed in the field.

Fig. 2. LIA outlines were digitised for Sygneskarsbreen, Austdalsbreen and ID 2485 based on 1966 orthophotos (left figure). The newest glacier outlines were based on Sentinel-2 images from 27 August 2019 (right figure). Lake Kupvatnet and lake Austdalsvatnet (part of the reservoir Styggevatnet) were regulated in 1988. Source: /norgeibilder.no/Copernicus Sentinel data 2019/.

1966 outline

Glacier outlines for 1966 were obtained from the Norwegian Water Resources and Energy Directorate (NVE) (Winsvold and others, 2014). The 1966 dataset was digitised based on the first edition 1: 50 000 topographical maps in the N50 series from the Norwegian Mapping Authority which were based on aerial photographs and previously used in Paul and others (2011). For Jostedalsbreen, the maps were based on aerial photographs from July 1966 and the outlines for five map sheets (1317-I, 1318-II, 1418-I, 1418-III, 1418-IV) covering Jostedalsbreen were digitised (Paul and others, 2011). Uncertainties in the 1966 glacier outline dataset are georeferencing and digitising accuracy as well as manual interpretations of snow and ice extent by the cartographers that produced the maps. Uncertainties within the interpretation of the maps of the digitiser are assumed small (Winsvold and others, 2014). There were some adverse snow conditions in the 1966 maps, in particular the northern map sheet (1418-IV) that may result in overestimation of the 1966 area (Paul and others, 2011). All images of the photo series from July 1966 (Contact no. WF1833) have later been orthorectified and are available at https://norgeibilder.no (since November 2022) (Hexagon, 2022).

We compared the 1966 orthophotos and the 1966 outlines and found some differences in the position of the glacier outlines, which can be attributed to geolocation errors in both the outline and in the orthophoto (Fig. 3). It could be possible to revisit the orthophotos and manually update the outlines, but we decided not to do so in our study because these inaccuracies are small compared to the uncertainties in the vicinity of the snowfields and as discrepancies also can be due to geolocation errors in the orthophotos now made available. Neither Paul and others (2011) nor Winsvold and others (2014) give uncertainty estimates for the 1966 outlines of Jostedalsbreen. We estimate an uncertainty of ~4% in the total area of the ice cap by applying a buffer of ±25 m around the perimeter. Individual glaciers may have larger uncertainties, depending on size and snow conditions along the margin. It is likely that the area is rather overestimated than underestimated due to the snow conditions in the 1966 aerial photographs.

Fig. 3. Subset of the 1966 orthoimages showing Lodalsbreen and Stigaholtbreen, with inset of Stigaholtbreen. The 1966 outlines are digitised from the N50 topographical maps. Coordinate system UTM 33N, datum ETRS_1989. Source: /norgeibilder.no/ Widerøe's flyveselskap AS/.

2006 outline

The 2006 outline was taken from the inventory that originated from a Landsat image from 16 September 2006 using a semi-automatic mapping method (Paul and others, 2011). This inventory was updated and included in the national glacier inventory by NVE with Jostedalsbreen divided into 82 glacier units (Andreassen and others, 2012). We refer to the final inventory dataset by Andreassen and others (2012) and not the dataset described by Paul and others (2011) due to some differences in what parts that were defined as Jostedalsbreen (see discussion). The uncertainty in the outlines for all of Norway was estimated to be within ±3%. For Jostedalsbreen, we estimate the uncertainty to be within one pixel (30 m). Applying a buffer of ±15 m gives an uncertainty in total area of ±2.5%.

2019 outline

The 2019 outline was obtained from the national glacier inventory derived from Sentinel-2 images taken 27 August 2019 using a standard semi-automatic mapping method selecting optimum threshold values of the red band divided by the shortwave infrared band and threshold in the blue band 2 (Andreassen and others, 2022). The dataset was checked using orthophotos and Sentinel composites. Manual edits were made to correct for parts in shadow, debris and lake outlines. Here, Jostedalsbreen was divided into 81 glacier units, two glaciers from 2006 (ID2367 and 2369) were disconnected from the main glacier since 2006 (see Fig. 4a) and one glacier unit was added (ID6762) due to shrinking of ID2399 (Fig. 4b). The uncertainty in the outlines is estimated to be within one pixel (10 m). Applying a buffer of ±5 m gives an uncertainty in total area of ±1.1% (Andreassen and others, 2022).

Fig. 4. Comparison of glacier inventory 2006 and 2019 detached parts and ice divides for (a) section around Bøyabreen where IDs 2367 and 2369 were included in Jostedalsbreen in 2006 and detached from it in 2019, and (b) section around Nigardsbreen where ID2299 was split into 2299 and 6762 in the new 2019 inventory and ice divides were updated. The dark grey outline shows what is included as Jostedalsbreen in 2019. Coordinate system UTM 33N, datum ETRS_1989.

Ice divides and glacier drainage basin homogenisation

For the 2019 glacier inventory, the ice divides from 2006 were used as much as possible, but divides were updated for some glaciers to separate units more naturally. Ice divides were also updated for Nigardsbreen (Fig. 4b) and Austdalsbreen and some neighbouring glaciers to fit with divides used in the glaciological mass-balance calculations (Andreassen and others, 2022). In this study, the ice divides for glacier drainage basins from 2019 were used to update the 1966 basins to make it consistent. For our study, this resulted in minor modifications for glacier basins around Nigardsbreen (2280, 2289, 2296, 2297) and Austdalsbreen (2478, 2480 and 2485). It should be noted that the mass-balance glacier basin of Nigardsbreen consists of IDs 2297 + 2299 + 2311 (Fig. 4b).

Outlet glacier centrelines

For the 1966 and 2006 inventories, glacier centrelines were calculated by Winsvold and others (2014) using a three-step costgrid-least-cost route approach which requires glacier outlines and a DTM as input (Kienholz and others, 2014). Glaciers smaller than 1 km² from the 1966 dataset were excluded in their study to reduce the noise from seasonal snow cover. In the most recent inventory from 2019 over Jostedalsbreen, centrelines were calculated using a geometrical routing algorithm described by Maussion and others (2019) based on applying the same three-step costgrid-least-cost route approach (Kienholz and others, 2014). The 2019 centrelines are available in the digital NVE Atlas (Andreassen, 2022).

A challenge with centreline methods, whether automatic or manual, is that position changes along a glacier front can be uneven and length calculations will therefore vary depending on the position of the centrelines. Furthermore, nunataks and shifting ice divides and hence changing drainage basins can cause problems when adjusting a centreline to multiple datasets. To ensure consistency in our calculations, we adjusted the centrelines derived from 2019 by manually lengthening them to fit the 2006, 1966 and LIA outlines (Fig. 5). We used the 2019 centrelines as they are based on the updated 2019 basins. In some cases, the centrelines were modified in the lower part to best fit the different glacier inventories. Some glaciers that were one glacier during the LIA maximum such as Bergsetbreen and Baklibreen were split to have one centreline for each present glacier (Fig. 5).

Fig. 5. Centrelines in LIA, 1966 and 2006 were digitised or adjusted using centrelines from 2019. Background orthophoto from 1966. Note the advance of Tuftebreen where the 2006 outline is outside the 1966 outline. Source: /norgeibilder.no / Coordinate system UTM 33N, datum ETRS_1989.

1966 DTM

The 1966 DTM was constructed from aerial photos acquired on 19 July and 21 July 1966 taken by Widerøe's flyveselskap AS at scale 1: 38 000 (Contract NoWF1833). The images had a forward overlap of 60% and a lateral overlap of 20%. The camera used was a Wild RC5/RC8 and the date of calibration was 17 January 1966. The images were scanned by the Norwegian Mapping Authority at a resolution of 12.5 μm resolution giving a pixel resolution of ~48 cm. The general image quality was variable, some of the images were blurry and had artefacts that probably originated from the development of the film. The image contrast also varied among the stripes, maybe due to varying scanning parameters of the different flight strips. Some steep parts outside glaciers were too dark to contain usable information, whereas in other photos parts of glaciers were too bright. There was no fresh snow in the images.

The images were further processed by Terratec As in 2021 (Terratec, 2021a, 2021b). First, a smaller part of the area around Austerdalsbreen was tested, before the remaining parts of the study area of broader Jostedalsbreen were processed. The 1966 DTM covers about 2/3 of Jostedalsbreen (Fig. 6).

Fig. 6. Overview of Jostedalsbreen dataset with the outlines from 1966 and 2019, the extent of the DTM from 1966 and 2020 and LiDAR flight lines over Jostedalsbreen on 9, 10 and 15 August 2020. Note the missing lowermost part of Tunsbergdalsbreen. Tun-Tunsbergdalsbreen, Nig-Nigardsbreen, Aus-Austdalsbreen. The inset shows the mapped part of Jostedalsbreen.

The aerotriangulation was performed in two steps. First, the internal orientation was solved using the camera calibration report, removing the radial lens distortion in the images. Thereafter, the external orientation was solved using 24 ground control points selected from the stable terrain visible in both the 1966 photos and in aerial orthophotos from 2017 (TT-14233, resolution 25 cm) that had been orthorectified earlier with accuracies of 0.17 and 0.14 m in horizontal and vertical directions, respectively. This resulted in a final root mean square error (RMSE) of 0.55 m in xy and 0.60 m in z for the 1966 DTM.

The result was manually edited to remove blunders in areas with poor contrast, resulting in some data voids holes in the 1966 DTM (Fig. 6).

2020 DTM

Airborne laser scanning (LiDAR) data were collected by Terratec over 3 d in August 2020 covering Jostedalsbreen and surrounding area with a minimum point density of 2 p m^–2. The central parts (blue) were scanned on 9 August, the north-western parts on 10 August (green) and the south-eastern parts on 15 August (Fig. 6). The LiDAR data were processed by Terratec and made available from https://hoydedata.no as pointcloud (las files) and a gridded raster (Terractec, 2020). We opted to use the 10 m resolution DTM which is now part of the national digital terrain model of Norway (DTM 10) by the Norwegian Mapping Authority. The 2020 DTM dataset covers all of Jostedalsbreen, except for the lower tongue of Tunsbergdalsbreen (Fig. 6). A quality control of the dataset showed that 99.6% of all 10 × 10 m grid cells had a point density of 2 p m^–2 or better. As the LiDAR survey was performed over several days, snow and ice melting cause some differences in surface elevation. Terratec adjusted for differences between flight strips and there was no systematic difference in elevation between different flight lines. The absolute georeferencing accuracy was determined by comparing the LiDAR data to control points taken from vehicle tracks measured with differential GNSS (Global Navigation Satellite System) along roads (Terratec, 2020). The comparison revealed a median difference of −0.16 m, which following adjustment was reduced to 0.00 m with a std dev. of 0.052 m and a dh of ±0.12 m. Overall, the absolute georeferencing accuracy of the final point cloud is assumed to be within ±0.1 m.

Other DTMs

NVE has airborne LiDAR surveys for two other dates as part of the mass-balance programme on Nigardsbreen and Austdalsbreen. For the Nigardsbreen LiDAR DTMs, the survey dates were 17 October 2009 and 10 September 2013, and for the Austdalsbreen LiDAR DTMs, the survey dates were 17 October 2009 and 27 August 2019.

The 2009 LiDAR survey, covering both Nigardsbreen and Austdalsbreen, had an average point density of 0.32 p m^–2. The LiDAR data were processed by Blom Geomatics AS (Blom Geomatics, 2009). The absolute georeferencing accuracy was determined by comparing the LiDAR data to seven control points measured by GNSS on the Nigardsbreen glacier surface. The comparison revealed an average difference of +0.16 m (Kjøllmoen, 2016).

The 2013 LiDAR survey, covering Nigardsbreen, Tuftebreen, Baklibreen, Bergsetbreen and Tunsbergdalsbreen, had an average point density of 1.0 p m^–2. The LiDAR data were processed by Terratec (Terratec, 2014). The absolute georeferencing accuracy was estimated to be 0.1 and 0.25 m in vertical and horizontal directions based on theoretical considerations (Terratec, 2014). The absolute georeferencing accuracy was also determined by comparing the LiDAR data to 13 control points measured by GNSS on the Nigardsbreen glacier surface. The comparison revealed an average difference of −0.10 m (Kjøllmoen, 2016).

The 2019 LiDAR survey, covering Austdalsbreen, had an average point density of 1.5 p m^–2. The LiDAR data were processed by Terratec. The absolute georeferencing accuracy was estimated to be 0.10 m (z) and 0.25 m (x,y) based on theoretical considerations (Terratec, 2019). The LiDAR dataset was compared with the LiDAR dataset from 2009 in stable areas. The comparison revealed a systematic difference of 0.05 m and accordingly the 2019 dataset is lowered 0.05 m (Terratec, 2019).

DTM differencing and calculation of geodetic mass balance

To ensure that the 1966 and 2020 DTMs were optimally aligned, we performed a co-registration routine following the method set out by Nuth and Kääb (2011) which seeks to minimise the RMSE slope normalised elevation biases over stable terrain. The 2020 DTM was considered the reference. The process was iterated until the std dev. (SD) over stable terrain changed by <2%. We implemented this co-registration through the demcoreg Python module (https://github.com/dshean/demcoreg) which resulted in minor shifts of −3.39, −3.47 and −0.93 m in the x, y and z directions of the 1966 DTM, respectively. Upon examining the elevation biases over stable terrain, no non-linear biases were observed.

Voids were filled using a global hypsometric polynomial method (after McNabb and others, 2019) where a third-order polynomial was fitted between 50 m elevation bins and the change in surface elevation. We opted to perform the void filling on Jostedalsbreen as one unit as opposed to running it on each glacier catchment as the latter would introduce erroneous jumps in pixel value at ice divides. The void-filled elevation change raster was then cleaned following the method set out by Gardelle and others (2013) whereby pixels were removed that were >3 SDs of the elevation change over stable terrain per 50 m elevation bins. The pixels removed by this process were primarily in areas of shadow or steep terrain.

We calculated the geodetic mass balance, B _geod:

(1)$$B_{{\rm geod}}\,{\rm} = {\rm \Delta }V\,\times \,f_{{\rm \Delta V}}{\rm /}A$$

where ΔV is the volume change and A is the average glacier area of the two surveys assuming a linear change in time. In the absence of glacier outlines from 2020, the average glacier area was calculated based on the outlines from 1966 and 2019. We used a conversion factor, f _ΔV, of 850 ± 60 kg m⁻³ that is considered appropriate to use for converting volume change to mass change for a wide range of conditions (Huss, 2013).

We calculated the uncertainty of the volume change (EΔv_i in m³) by summing up the standard error (EΔh_i in m) per 50 m elevation bins, multiplied by the area of each elevation bin (A_i) to account for the hypsometry:

(2)$$E\Delta v_i = \mathop \sum \limits_i^n E\Delta h_i \times A_i$$

where EΔh _i is derived from the std dev. over stable ground (σ _stable), divided by the effective number of observations (N) (Bolch and others, 2011)

(3)$$E\Delta h_i = \displaystyle{{\sigma _{{\rm stable}}} \over {\sqrt N }}$$

N is calculated using the number of pixels (N _tot) in the DTM differencing, the pixel size (PS) and the distance of spatial autocorrection, d, which is commonly assumed to be equal to 20 pixels (Rounce and others, 2018; King and others, 2019).

(4)$$N = \displaystyle{{N_{{\rm tot}} \times {\rm PS}} \over {2d}}$$

The total uncertainty of the geodetic mass balance (E _Bgeod) also considered uncertainties relating to the volume to mass conversion (E_ρ), which following Huss (2013) was taken to be ±60 kg m⁻³. We also assumed an error in delineating the glacier outline (E _a), the ±25 and ±5 m for the 1966 and 2019 outlines translate to ±3.8 and ±1.1% uncertainties in glacier area over the entire ice cap, resulting in an area uncertainty of ±4.0% when the two DTMs are compared. Therefore, we have:

(5)$$E_{{\rm Bgeod}} = \sqrt {E\Delta v_i^2 + E_\rho ^2 + E_a^2 } $$

Estimation of seasonality correction and dissipative melting

To facilitate the comparison between geodetic and glaciological mass-balance records, we applied a temperature-index model to estimate, for each of the outlet glaciers, internal and basal ablation from dissipation of potential energy due to ice and water flow over the period 1966–2020 and seasonality correction from the mapping dates to the end of the 1966 and 2020 melt seasons. The seNorge_2018 dataset of gridded daily total precipitation and daily mean temperature from the Meteorological Institute of Norway (Lussana and others, 2019) was used as input to the model. The model was run on the 1 km DTM of the seNorge_2018 model with a time series of glacier masks created by intersecting the 1966, 2006 and 2019 glacier outlines with the DTM. Each mask was used for half of the period before and after the mapping year.

Daily melt, M_i, in a gridcell was calculated as the sum of temperature and radiation-based components when the seNorge temperature in the cell, T_{sn ,i}, was above 0°C (Hock, 1999):

(6)$$M_{{\rm snow/ice}, i} = \left\{ \matrix{ ( {\rm MF} + C_{\rm snow/ice}\,I_i ) T_{sn, i}, &T_{sn, i} > \,0^\circ {\rm C}\cr 0, &T_{sn, i} \le \,0^\circ {\rm C}\cr }\right. , $$

where MF is the melt factor and C _snow/ice are radiation coefficients relating the amount of melt of snow and ice in a gridcell to the potential direct incoming solar radiation, I_i, in the cell. Firn melt in a gridcell was computed as the average of the potential snow and ice melt. Daily accumulation in each gridcell was calculated from the daily total precipitation, assuming that the fraction of precipitation that falls as snow is a linear function of the daily mean temperature in a 2°C interval around 1°C. In addition to MF and C _snow/ice, we calibrated global precipitation and winter temperature correction (P _corr and T _corr) to the seNorge_2018 daily precipitation and winter temperature as these variables are considered uncertain over remote, mountainous areas (Lussana and others, 2019).

Potential energy released by downward motion of meltwater and ice is likely the largest heat source for internal ablation (Cogley and others, 2011). For a temperate glacier, we assume that all of the energy dissipated due to the flow of ice and water contributes to meltwater production within or at the base of the glacier. In this study, we estimate dissipative melt associated with loss in potential energy following Oerlemans (2013). Hence, internal and basal ablation due to dissipative melting was calculated from the difference in the potential energy of precipitation deposited over the glacier surface and the potential energy the precipitation has as it leaves the glacier snout in the form of meltwater runoff, assuming a steady-state glacier and that the water that enters and leaves the glacier is at the freezing point. Internal and basal ablation from dissipation of potential energy over a glacier, M _int, was calculated as the sum of potential energy loss over each glacier gridcell as:

(7)$$M_{{\rm int}} = \displaystyle{{\mathop \sum \nolimits_i^{} g\;P_iA_i( {h_i-h_{{\rm min}}} ) } \over {A_{{\rm tot}}L_m}}\;$$

where P_i is the corrected seNorge_2018 precipitation in cell i, A_i and h_i are the glacierised area and elevation of the cell, respectively, h _min is the minimum elevation of the glacier, A _tot is the total area of the glacier, g is the acceleration of gravity and L_m is the latent heat of fusion.

The mass-balance model was calibrated following the Monte Carlo approach of Engelhardt and others (2014) with 25 000 model runs. In situ seasonal surface mass-balance measurements are available for five outlet glaciers of Jostedalsbreen: Nigardsbreen (1962–2020; 59 years), Austdalsbreen (1988–2020; 33 years), Supphellebreen (1964–1967; 4 years), Tunsbergdalsbreen (1966–72; 7 years) and Vesledalsbreen (1967–1972; 6 years). To find the best global parameter sets (P _corr, T _corr, MF, C _snow and C _ice) for the entire ice cap, every second year of the mass-balance records was used for model calibration or validation. The resulting coefficient of variation between the modelled and observed glacier-wide winter (summer) balance was 0.20 (0.14–0.16) for the calibration years and 0.22–0.23 (0.17–0.19) for the validation years for the 100 parameter sets with the best ranked performance (Table S1). Dissipative melting and glacier-wide seasonality corrections for 1966 and 2020 to the end of the melt season were computed as the mean of the model results from the 100 runs with the best performing parameter sets (see Tables S1–S2). The end of the melt season in 1966 and 2020 was defined as the date at which the modelled net accumulation over each outlet glacier was at a minimum.

Results

Length changes

Length changes along central flowlines for the 18 outlet glaciers with mass balance or front variation measurements reveal an overall reduction in length of 2.6 km, or 28%, from their mean LIA length of 9.6 km (Fig. 7). For the subperiods LIA–1966, 1966–2006 and 2006–2019, the total change (change per year) is 2047 m (−10 m a⁻¹), 253 m (−6 m a⁻¹) and 343 m (−26 m a⁻¹), respectively. Although the time spans are very different in the subperiods (211, 40 and 13 years, respectively), the data show that for all the glaciers the fastest retreat among these subperiods is for the last period 2006–2019.

Fig. 7. Length change along flowlines for 18 outlet glaciers of Jostedalsbreen ice cap with mass balance or front variation measurements (see Table 1).

The dataset captures the advances of Bødalsbreen, Tuftebreen, Briksdalsbreen and Melkevollbreen between 1966 and 2006, and suggests advances of Bergsetbreen and Supphellebreen between 2006 and 2019. The advances of Bødalsbreen, Briksdalsbreen and Supphellebreen are also visible in the front variation records (NVE, 2023), whereas for the other glaciers the advances are not recorded in the front variation measurements. For Supphellebreen, the centreline reveals a net advance of 32 m from 2006 to 2019, whereas the front variation measurements show a net retreat of 37 m between 2006 and 2014 (4 years with advance, 4 years with retreat) whereafter the observations were terminated. Although the time period is not the same and the positions of front variation points and centreline are not identical, the positive change for this glacier is likely an error following the coarser satellite-based outlines of 30 and 10 m resolution, respectively.

Annual front variation observations show that Nigardsbreen retreated 3.0 km since measurements started in 1899. Between 1966 and 1975, the glacier retreated 461 m, and then slowed down and retreated only 45 m until 1988. Nigardsbreen advanced 259 m between 1988 and 2000. The glacier retreated slowly until 2006 (−8 m a⁻¹) and then faster up to 2019 (−41 m a⁻¹). Several other glaciers experienced the same pattern. Only the longest glaciers, Tunsbergdalsbreen and Lodalsbreen, did not advance in the 1990s. On the other hand, several glaciers, i.e. Briksdalsbreen, Bergsetbreen and Brenndalsbreen advanced both in the 1960s and 1970s and in the 1990s.

Area changes

Our addition of LIA outlines for three more glaciers resulted in a LIA area increase of 0.8% from 568 to 572 km² compared to the previous LIA inventory of Carrivick and others (2022). Jostedalsbreen ice cap has reduced from 572 km² at its LIA maximum extent to 458 km² in 2019 or 456 km² when including disconnected parts (Table 2). This gives a total area reduction of 116 km² or 20%, or an area change rate of −0.08% a^–1 using the median LIA year 1755. Similarly, the 1966 glacier area is smaller when only including parts within the LIA limits (Table 2). The reason for this is that the LIA inventory used the 2006 inventory as basis and not 1966 which extends beyond LIA in parts of Jostedalsbreen. The reduction in area from 1966 to 2019 within the LIA limits is 29.5 km², or 6%, showing an area reduction of −0.11% a⁻¹. The relative changes are thus larger from 1966 to 2019 (−0.11% a⁻¹) compared to the period from LIA to 1966 (−0.08% a⁻¹). However, it must be emphasised that the total area and calculated area changes will vary depending on what mapping year is used as reference and what is considered part of the continuous Jostedalsbreen. Moreover, change rates will depend on the LIA year used for mapping.

Table 2. The total area, area change and area change rates relative to LIA (~1755) of the Jostedalsbreen ice cap

n denotes the number of glacier units. Area* includes disconnected parts that are within the LIA outlines.

The total area of the 49 glacier units covered by DTM differencing represents 73% of the total area of Jostedalsbreen, and this percentage is nearly the same in 1966, 2006 and 2019 (Table 3). The sample of 49 glaciers has reduced in total area from 363.4 km² in 1966 to 332.9 km² in 2019, a reduction of 30 km² or 8.4%. The percentage reduction is the same as for the total area of Jostedalsbreen. Looking at the subperiods 1966–2006, and 2006–2019 also give comparable values in total percentage reduction.

Table 3. Glacier area (A) of the sample of 49 glaciers and all Jostedalsbreen for 1966, 2006 and 2019

Over the period of DTM comparison, 1966–2019, all glaciers reduced in size except for one (ID 2255) (Table 4, Fig. 8). The reason for this is that the 1966 outlines treat this glacier as two parts: 2255 connected to Jostedalsbreen and a disconnected part of size 0.57 km². According to the orthophotos from 1966, the disconnected part should be included in the continuous Jostedalsbreen. The Sentinel-2 image used for the 2019 outlines has heavy shading and it remains uncertain whether the parts were connected. Looking at subperiods, five of the 49 glaciers had an increase in area from 1966 to 2006 and seven had an increase in area from 2006 to 2019. The increase in area for some of the glaciers is often related to variable inclusion of parts attached to the glacier and as such length changes provide a more accurate notion of glacier advance. All glaciers reduced their area from the LIA maximum extent (Fig 1).

Fig. 8. Change in area between 1966 and 2019 for the 49 study glaciers covered by DTM in 1966 and 2020.

Table 4. Area (A) in 1966 and 2019, elevation changes (ΔH), geodetic mass balance (B _geod) and estimated error (E _Bgeod) in geodetic balance for the 49 glaciers compared between the 1966 DTM and the 2020 DTM

Values are shown both without and with void filling, void filled values are noted vf. Glacier IDs and names from the glacier inventory (Andreassen and others, 2012; 2022).

Surface elevation changes and geodetic mass balance

The surface elevation changes over the 54-year period from 1966 to 2020 show large spatial variability among the glaciers (Fig. 9, Table 4). The most negative changes are for the detached part of Brenndalsbreen (ID 2301) (−43 m), Vesledalsbreen (−30 m), Lodalsbreen (−28 m) and Tunsbergdalsbreen (−24 m) when using void-filled values. Nigardsbreen including basins 2299 and 2311 has a surface change of −4.0 m, the part of 2297 only is −3.0 m. Austdalsbreen has a surface change of −19 m. The overall change of all the 49 glaciers is −10.2 m (−0.19 m a⁻¹) with voids filled and −10.8 (−0.20 m a⁻¹) without filling.

Fig. 9. Elevation differences of Jostedalsbreen from 1966 to 2020. The insets show details of Lodalsbreen (upper right) and Briksdalsbreen and Brenndalsbreen (lower right). The glacier extents 2019 and 1966 are shown for the Jostedalsbreen ice cap. Background: mountain shadow from the 100 m national DTM.

Three of the glaciers (IDs 2485, 2488 and 2489) have overall positive balances. These are small, higher elevated glaciers in the northern part of the ice cap and might be caused by errors in the 1966 DTM due to poor contrast in the orthophotos.

Converting the elevation change values into geodetic mass balance shows that the mean area-weighted value of all the 49 glaciers is −0.15 (0.16) ± 0.02 m w.e. a⁻¹ using void-filled (non-void-filled) values. The geodetic mass balance is less negative for Nigardsbreen than the overall and more negative for Austdalsbreen and Tunsbergdalsbreen than the overall (Fig. 10).

Fig. 10. Geodetic mass balance for Nigardsbreen (Nig), Tunsbergdalsbreen (Tun), Austdalsbreen (Aus) and all 49 glaciers covered by the 1966–2020 DTM differencing (All). Bnf – not void filled. Bvf – void filled. Results are for the period between the surveys and not adjusted for seasonal correction. See Table 4.

Examining the changes by elevation bands shows that the large thinning of the Tunsbergdalsbreen tongue (Fig. 9) affects the overall hypsometric curve for Jostedalsbreen with its tongue covered by bands from 600 m a.s.l. and higher (Fig. 11a). Since data are missing for the lowermost part here due to lacking LiDAR data from 2020 (Figs 6, 9), an inclusion of this part would probably change the curve in the lower intervals due to the area reduction and large thinning here.

Fig. 11. Comparison of mean surface elevation change per elevation bin for Jostedalsbreen ice cap and Austerdalsbreen, Nigardsbreen and Tunsbergdalsbreen between (a) this study and (b) Hugonnet and others (2021). Shaded error bands represent 1 std dev. of elevation changes per elevation bin. Only glacier catchments fully covered by both datasets were included in Jostedalsbreen.

Seasonality correction for 1966 and 2020

Dates of modelled minimum mass balance of the 49 surveyed outlet glaciers on Jostedalsbreen in 1966 ranged from 6 to 22 September, except for lower parts of Briksdalsbreen (ID 2301) where the minimum date was estimated to 24 October 1966. In 2020, estimated minimum dates ranged from 4 September to 17 October (Table S2).

The date of modelled minimum mass balance for Nigardsbreen in 2020 was estimated to be 17 October, which gives an additional melt of −0.57 ± 0.04 or −0.36 ± 0.04 m w.e. from the respective mapping dates of higher (9 August) and lower (15 August) elevations in 2020 until the end of the ablation season. For 1966, the date of minimum was estimated to be between 7 (ID 2297) and 14 (ID 2299 and 2311) September, giving a total correction of −1.01 ± 0.04 m w.e. For Austdalsbreen, the minimum date in 2020 was estimated to be 17 October, which amounts to a correction of −0.71 ± 0.05 m w.e. from the mapping date 9 August. For 1966, the estimated minimum date was 14 September, which constitutes a seasonality correction of −1.22 ± 0.09 m w.e. Kjøllmoen (2022) estimated melt of −0.57 m w.e. from 9 August to 15 October 2020 (date of glaciological measurements of minimum mass balance) for Nigardsbreen (basin consisting of ID 2297, 2299, 2311 and 6762) in the reanalysis of the glacier's mass-balance record, which corresponds well with our estimated seasonality correction between 9 August and 17 October 2020.

Discussion

How large is Jostedalsbreen?

The digital glacier outlines used in this study revealed reductions in the total area of the ice cap, from 572 km² at its LIA maximum extent, to 500 km² in 1966, 474 km² in 2006 and 458 km² in 2019 (Table 2). As demonstrated, both the total areas and calculated change rates will differ depending on the reference used for comparison. The available digital sources of both the outlines and the source material used to derive the outlines (orthophotos, topographic maps and satellite images) are key to be able to compare the data properly for individual glaciers, but this is also time-consuming for the analyst.

The accuracy in area will depend on the spatial resolution of the source. Satellite imagery with a finer resolution may make it easier to differentiate between continuous and discontinuous ice bodies but may also require more editing (Paul and others, 2016). Studying glacier changes over longer time periods requires the use of different sources which can represent challenges due to different methodologies influencing glacier mapping results. To compare the changes in Jostedalsbreen over time is of value, but few sources are available. The maritime climate of Jostedalsbreen with persistent clouds and in many years heavy seasonal snow cover has led to very few suitable satellite scenes being available. As a result, only one inventory is so far available from the Landsat era (Paul and others, 2011; Winsvold and others, 2014).

Another point to mention is that the size of Jostedalsbreen is reported differently in other studies using the same sources. In the first list of numbers and areas of glaciers in Norway by Liestøl (1962), the size of Jostedalsbreen was measured from a county map of scale 1: 200 000 based on surveys from 1863–71 and revised by aerial photographs from 1945 available at that time. The resulting area of ‘The main plateau of Jostedalsbreen with outflows’ was 473.3 km², thus smaller than our 1966 area. A more detailed inventory was published in 1969 based on the same aerial photographs from 1966 and topographic map sheets. Here, Jostedalsbreen was mapped to be 486.3 km² and for the first time was divided into 45 glacier units (p157–160 in Østrem and Ziegler, 1969). An updated inventory using the same methods was published in 1988 based on aerial photographs from 1984 (Østrem and others, 1988). Here, Jostedalsbreen was divided into 63 glacier units totalling an area of 486.7 km². Both these inventories were analogue and available only as digital tabular data, not as digital outlines.

The 1966 photographs used in the 1969 inventory were used to create the main maps of Norway that were digitised by Paul and others (2011). The 1966 estimate of Jostedalsbreen is 520 km²_, whereas we report 500 km² for the continuous ice cap. The difference is due to Winsvold and others (2014) using the 2006 inventory as baseline and therefore excluded parts covering Tindefjellsbreen and Skålebreen in the north-west. Results can differ depending on what is included or excluded as the study domain and on how to define the continuous ice cap.

Glaciological versus geodetic mass balance

Nigardsbreen is the only glacier with glaciological mass-balance data covering the entire period 1966–2020. As mentioned, previous comparisons of glaciological and geodetic mass balance resulted in calibration of the data for periods after 1984. The homogenised and calibrated glaciological mass balance for the period 1966/67 to 2019/20 for Nigardsbreen shows a mass surplus of 1.4 or 0.03 m w.e. a⁻¹ during this period. Using just the homogenised values for the whole period gives a more positive surplus of 16.6 or 0.25 m w.e. a⁻¹. Our geodetic mass balance is −3.52 m w.e. when voids are not filled and −3.41 m w.e. when voids are filled. Adjusting for the additional melt in 1966 (−1.01 m w.e.) and 2020 (−0.57 m w.e.) results in geodetic mass balances of −3.08 and −2.97 m w.e. for not filled and filled voids, or −0. 057 and −0.055 m w.e. a⁻¹, respectively. Our geodetic results are thus more negative than the homogenised and calibrated values for this period and the differences are even larger using the non-calibrated values and confirm the discrepancy found in the previous studies (Andreassen and others, 2016; Kjøllmoen, 2022). It should be noted that the period from 1966 to 2020 is not directly comparable with the subperiods used in reanalyses: 1964–1984–2009–2020 (Kjøllmoen, 2016, 2022). The differences between the homogenised glaciological mass balance and geodetic mass balance are partly ascribed to internal and basal ablation not being accounted for in the glaciological method, as well as the spatial interpolation to convert individual point measurements to glacier-wide values (e.g. Zemp and others, 2013). Density conversions are an uncertainty in both the glaciological and geodetic methods. As is standard with geodetic mass-balance observations, this study assumed a fixed density assumption. Observations of thickness and density of firn are in general uncommon in Norway so it is difficult to validate our density assumption. One 47 m core from 1987 is available from Nigardsbreen which revealed densities varying between ~500 and ~900 kg m⁻³ up to a firn/ice transition at 30 m depth (Kawamura and others, 1989). This data are however a single point measurement, and in the absence of repeat data, we cannot use these data to validate our density assumption or speculate whether changes to the thickness and extent of the firn reservoirs over time could influence our geodetic results.

Internal accumulation is another source of uncertainty that is not accounted for in the glaciological method and may also play a role (e.g. Schneider and Jansson, 2004). Our estimate of internal and basal ablation from dissipation of potential energy for Nigardsbreen) amounts to −6.06 m w.e. or −0.11 m w.e. a⁻¹, which accounts for less than half of the difference between the glaciological and geodetic results. Estimates of dissipative melting for Nigardsbreen vary over the period 1966–2020, with minimum and maximum values of −0.19 m w.e. in 1990 and −0.06 m w.e. in 1977, respectively. The year 1990 corresponds to a year of higher-than-average winter mass balance in the glaciological records of Nigardsbreen, while 1977 coincides with a year of lower-than-average winter mass balance (Fig. 12).

Fig. 12. Annual balance (Ba) values for the five Jostedalsbreen outlet glaciers with glaciological mass-balance records. Upper figure shows the period 1962–1987, and lower figure shows the period 1988–2020. Ba values are colour coded.

Overall, the mean estimated internal and basal ablation from dissipative melting for the part of Jostedalsbreen covered by the 1966 DTM over the period 1966–2020 is −0.07 m w.e a⁻¹ (Table S2). Generally, dissipative melting is larger for larger outlet glaciers with significant elevation differences (Fig. 13), such as Nigardsbreen. The correlation between internal ablation, accumulation and elevation difference is as expected given the linear relationship between these variables (Eqn (7)).

Fig. 13. Estimated dissipative melting versus elevation range with glaciers displayed according to area. Aus-Austdalsbreen, Tun-Tunsbergdalsbreen and Nig-Nigardsbreen.

Over the period 1966–2020, mean annual values for individual outlet glaciers range from −0.12 m w.e. a⁻¹ (Austerdalsbreen) to −0.01 m w.e. a⁻¹ (Strupbreen). However, due to the simplicity of the method and uncertainty in precipitation input data over complex terrain, values should be considered as coarse estimates. We do not account for dissipation of energy from meltwater that may enter the glacier from surrounding ice-free areas. For smaller glaciers at high elevations, we consider this contribution to be negligible. However, for valley-terminating glaciers with large catchments and elevation differences, such as Nigardsbreen, the contribution of meltwater from ice-free areas may not be negligible. In the reanalysis of the glaciological mass-balance records, the rate of dissipative melting from flow of ice and water at Nigardsbreen over the period 1964–2013 was quantified to −0.16 m w.e. a ⁻¹ (Andreassen and others, 2016). The uncertainty was assumed to be one-third of the estimated dissipative melting, which amounts to ±0.05 m w.e. a⁻¹. For Austdalsbreen, the rate of dissipative melting was estimated to be −0.03 ± 0.01 m w.e. a⁻¹ over the shorter period 1988–2009, while our results reveal an average rate of −0.02 m w.e. a⁻¹ over the period 1966–2020. Since we employ the same principle for assessing internal and basal ablation and regard the uncertainty in our input data to be comparable, we also consider the uncertainty in our estimates comparable.

For Nigardsbreen and Austdalsbreen, our estimates of dissipative melt are similar in magnitude to those of Andreassen and others (2016). Their slightly higher estimates may be attributed to their use of precipitation data from the seNorge dataset (v. 1.1), which shows overall larger precipitation amounts than the newest version seNorge_2018 used in this study. In addition, we calibrate seNorge_2018 precipitation to available estimates of winter balance for outlet glaciers of Jostedalsbreen and calculate dissipative melting over the model grid, rather than over the elevation intervals used in the calculation of surface mass balance in glaciological records. Oerlemans (2013) estimated typical rates of 0.23 m w.e a⁻¹ for Nigardsbreen, thus twice the amount of our estimate reflecting the uncertainty in such estimates.

Nevertheless, our results show that melting due to dissipation of potential energy in the flow of ice and water is a significant contribution to mass change of glaciers with high accumulation rates and large elevation differences.

Spatio-temporal variability in glacier change of Jostedalsbreen

Our results highlight a large variability in glacier response for the 49 glaciers studied between 1966 and 2020, with geodetic mass balances ranging from −0.54 ± 0.07 m w.e a⁻¹ (ID 2301) to 0.05 ± 0.01 m w.e a⁻¹ (ID 2489). In general, thinning is more pronounced in the north-east and the smaller area in the north-west (Fig. 9). To examine the topographic drivers of this change, we calculated statistics (slope, aspect, minimum, median and maximum elevation, elevation range and solar radiation) from the 1966 and 2020 DTM and compared with the geodetic results (Table S3). From these DTMS, we found that glaciers that have the lowest median glacier elevation have more negative geodetic balance than glaciers with higher median elevation (Fig. 14a). We also calculated a hypsometric index (HI) based on the maximum (H _max), median (H _med) and minimum (H _min) elevations of each glacier following the method of Jiskoot and others (2009):

(8)$${\rm HI} = {\rm \;}\displaystyle{{H_{\rm max}-{\rm \;}H_{\rm med}} \over {H_{\rm med}-{\rm \;}H_{\rm min}}}{\rm \;and\;if\;}0 < {\rm HI} < 1{\rm \;then\;HI} = \displaystyle{{-1} \over {{\rm HI}}}$$

Fig. 14. Geodetic mass balance for the 49 glaciers plotted against (a) median glacier elevation and (b) hypsometric index category (HI). The boxes reflect 25 and 75% percentiles of the mass balance, the lines within the boxes reflect the median mass balance, and the whiskers reflect the spread of the data beyond the interquartile range. Categories: VTP, very top heavy; TH, top heavy; ED, equidimensional; BH, bottom heavy; VBH, very bottom heavy; n, number of glaciers in each category. Median glacier elevation and HI are calculated from DTM1966 (see Table S3).

This allowed us to characterise glaciers as very top heavy (HI < −1.5), top heavy (−1.5 < HI < −1.2), equidimensional (−1.2 < HI < 1.2), bottom heavy (1.2 < HI < 1.5) or very bottom heavy (HI > 1.5).

Of these topographic drivers, median elevation in 1966 and 2020 showed the strongest correlation coefficient with geodetic mass balance (0.42 and 0.40, respectively), and HI in 1966 had the third highest correlation coefficient (0.18) (Table S3). Other drivers were not that well correlated. In general, glaciers that have the lowest median glacier elevation have more negative geodetic balance than glaciers with higher median elevation (Fig. 14a). By binning glaciers by their hypsometric category, it can be seen that glaciers with very bottom heavy hypsometries had on average more negative mass balances (Fig. 14b), although there were only two glaciers within this category (Strupebreen and detached part of Brenndalsbreen (ID 2301)) (Table S3).

Temporal variability of the geodetic mass balance of Jostedalsbreen over our study period July 1966 to August 2020 is difficult to assess due to a lack of observations covering subperiods of our analysis. The global geodetic mass-balance dataset by Hugonnet and others (2021) using multiple ASTER stereo imagery offers a shorter-term context for the ice cap over the period 1 January 2000 to 1 January 2020. They provide cumulative monthly values and annual to decadal rates for every glacier of the world using outlines of the Randolph Glacier Inventory (RGI) 6.0 that for Jostedalsbreen are the same as our 2006 glacier outlines (RGI, 2017). A comparison with their results is useful to discern changes in glacier mass loss rates over the different study periods and shows that surface elevation changes (for the same 49 glaciers that we analysed) are more negative between 2000 and 2019 than between 1966 and 2020 (Fig. 11). The mean (median) elevation changes measured by both methods are −0.16 (−0.11) m w.e. a^–1 (void-filled) for 1966–2020 period and −0.63 (−0.56) m w.e. a⁻¹ for the 2000–2019 period.

Comparison of our mean surface elevation changes with those of Hugonnet and others (2021) also reveals that the latter period is more negative with enhanced thinning of the lower part of the glacier outlets (Fig. 11b). The mean curve for all 49 glaciers studied reveals thinning of all elevation bands in the 2000–2019 period.

A recent study of the lower section of Austerdalsbreen from 1966 to 2021 covering six subperiods reveal enhanced surface lowering of this outlet glacier. The rates of surface lowering increased (for each period) since the period 1986–1997, with a maximum change rate occurring in the last observation period 2020–2021 (Seier and others, 2023).

The glaciological mass-balance records of Jostedalsbreen reveal interannual variations as well as differences between the glaciers (Fig. 12). The mean annual mass balance for the joint period 1967–72 (6 years) for Nigardsbreen, Tunsbergdalsbreen and Vesledalsbreen was +0.07, +0.05 and −0.39 m w.e. a⁻¹. The small difference between Nigardsbreen and Tunsbergdalsbreen over this period contrasts with the large difference in geodetic mass balance between the two glaciers between 1966 and 2020, −0.05 and −0.36 m w.e. a⁻¹, respectively.

The glaciological mass-balance record of both Nigardsbreen and Austdalsbreen shows a transient surplus in the period 1989–1995, thereafter many years with negative annual balances, but marked mass surpluses were recorded in 2012, 2015 and 2020 (Fig. 12). The mean annual mass balance for the joint period of observations 1988–2021 (34 years) is −0.50 m w.e. a⁻¹ for Austdalsbreen and +0.01 m w.e. a⁻¹ for Nigardsbreen.

Overall, the front observation and length change records, the glaciological mass balance records, the data from Austerdalsbreen surface lowering and our comparison with the results of Hugonnet and others (2021) for 2000–2020, all evidence that glacier change rates have become more negative since 2000. Although the sensors and measurements differ, all results reveal increased glacier thinning, more negative glacier mass balance(s) and a faster reduction in length and area since 2000. However, as mentioned above, there have been years with positive mass balances, as demonstrated in the glaciological records for Nigardsbreen and Austdalsbreen. A similar tendency of transient mass surplus was found for southern Folgefonna between 2013 and 2017 (Andreassen and others, 2022).

Hugonnet and others (2021) provide cumulative monthly values for every glacier of the world and we use these data to compare with the airborne LiDAR surveys of Nigardsbreen (17 October 2009 and 10 September 2013) and Austdalsbreen (17 October 2009 and 27 August 2019) to see how the surface elevation changes compare between the global dataset and very accurate repeat airborne LiDAR. For the LiDAR DTM differencing calculations, we used the RGI basins 6.0 (RGI consortium, 2017) used by Hugonnet and others (2021) to ensure that identical basins were used. It should be noted that the dates of the ASTER images start 1st of each month without the exact dates disclosed. We averaged the value for the months prior to and after the dates of the LiDAR surveys. Nigardsbreen and adjacent glacier ID 2299 and 3311 for the period 2009–2013 shows a surface lowering of 4.00 ± 0.14, 4.03 ± 0.14 and 4.54 ± 0.14 m over the period 2009–2013 (Kjøllmoen, 2016), whereas ASTER results in the same period show a lowering of 2.0 ± 1.75, 2.2 ± 2.59 and 2.6 ± 2.97 m, respectively (Fig. 12). For Austdalsbreen, the surface lowering calculated from LiDAR DTMs (2009–2019) is 11.5 ± 0.14 m (NVE unpublished, this study), whereas the ASTER results in the same period show a lowering of 7.7 ± 2.69 m. The uncertainties in the LiDAR DTMs are <0.1 m per DTM, whereas the errors in the dh cumulative estimate for the individual months range between 1.2 and 2.2 m, which gives a higher overall error due to error propagation. The short timescale and the smaller spatial domains are affected by temporal autocorrelation and are shorter than the 5-year periods used in the Hugonnet study (Hugonnet and others, 2021). Even if we are comparing the elevation differences and thus avoid density conversion uncertainties associated with geodetic mass-balance values for short periods (Huss, 2013), the uncertainties are much larger for the repeat ASTER than the repeat LiDAR, and the results must be interpreted with caution. Nevertheless, results for both Nigardsbreen and Austdalsbreen show that the LiDAR differencing reports much more (up to twice as much) negative balance than the repeat ASTER results (Fig. 15).

Fig. 15. Comparison of elevation change in between repeat LiDAR (this study) versus ASTER comparisons from Hugonnet and others (2021). Error bars calculated from propagation of errors using errors from the elevation DTMs for each mapping year. The mapping period for Nigardsbreen (ID 2297, 2299 and 2311) is 2009–2013 and for Austdalsbreen (Aus 2478) is 2009–2019.

Conclusions

Jostedalsbreen is mainland Europe's largest glacier and accommodates 20% of all glacier ice in Norway. This study has revealed reductions in the total area of the ice cap, from 572 km² at its LIA maximum extent (median year of 1755), to 500 km² in 1966, 474 km² in 2006 and 458 km² in 2019. Total areas and change rates will differ depending on the reference used for comparison. In total, the ice cap area is reduced since the LIA by −20 or −0.08% a⁻¹. The size of Jostedalsbreen is reported differently in other studies using the same sources. Such deviations exemplify how results can differ depending on which parts to include or exclude as the continuous ice cap and highlight the difficulties in comparing glacier inventories derived from analogue sources.

Length changes for 18 outlet glaciers with front variation measurements reveal an overall reduction, a change of −2.6 km, or −28%, from their mean LIA length of 9.6 km. Although the time spans are very different in the periods studied (211, 40 and 13 years, respectively), the fastest retreat rate of the glaciers was during the period 2006–2019.

Between 1966 and 2020, Nigardsbreen had a less negative geodetic mass balance (−0.05 m a⁻¹) than most of the rest of Jostedalsbreen, whereas Austdalsbreen (−0.28 m a⁻¹) and Tunsbergdalsbreen (−0.36 m a⁻¹) had a more negative mass balance than the mean. Glaciers with lower median elevation had more negative geodetic mass balance. Calving has enhanced the mass loss of the Austdalsbreen. A comparison of repeated airborne LiDAR for shorter time periods in the 2000s shows much more negative mass balance for Nigardsbreen (2009–2013) and Austdalsbreen (2009–2019) than suggested by repeat ASTER image processing. The glaciological mass-balance records for the joint period of observations between 1988 and 2020 (33 years) also reveals the difference in mass change between the glaciers with mean values −0.46 m w.e. a⁻¹ for Austdalsbreen and +0.027 m w.e. for Nigardsbreen. After 2000, both glaciers experienced mass deficits, but marked mass surpluses were recorded in 2012, 2015 and 2020 illustrating temporal variations in mass balance.

Our geodetic results are more negative than the homogenised and calibrated values for Nigardsbreen for this period in line with previous studies. Internal and basal ablation accounts for some of the differences and were estimated to be −0.11 ± 0.04 m w.e. a⁻¹ for Nigardsbreen, which is a bit lower than estimated in previous studies and is sensitive to the precipitation dataset used. The mean estimated internal and basal ablation from dissipation of energy for the 49 outlet glaciers covered by the 1966 DTM was on average −0.07 w.e. a⁻¹, with overall higher rates for larger glaciers with the largest surface elevation differences. Our results show that melting due to dissipation of potential energy in the flow of ice and water is a significant contribution to mass change of glaciers with high accumulation rates and large elevation differences.

Overall, these datasets of surface elevation change, area and length changes, glaciological and geodetic mass balance are an important reference that can be used for modelling the past and future development of Jostedalsbreen. Repeat surveys of surface elevation changes are recommended to determine the evolution of individual glaciers and Jostedalsbreen as a whole.