Spatiotemporal mass-balance variability of Jostedalsbreen Ice Cap, Norway, revealed by a temperature-index model using Bayesian inference

Abstract

Jostedalsbreen in western Norway is the mainland Europe's largest ice cap and a complex system of more than 80 glaciers. While observational records indicate a significant sensitivity to climate fluctuations, knowledge about ice-cap wide spatiotemporal mass changes and their drivers remain sparse. Here, we quantify the surface mass balance (SMB) of Jostedalsbreen from 1960 to 2020 using a temperature-index model within a Bayesian framework. We assimilate seasonal glaciological SMB to constrain accumulation and ablation, and geodetic mass balance to adjust model parameters for each glacier individually. Overall, we find that Jostedalsbreen has experienced a small mass loss of −0.07 m w.e. a⁻¹ (−0.21 to +0.08 m w.e. a⁻¹), but with substantial spatiotemporal variability. Our results suggest that winter SMB variations were the main control on annual SMB between 1960 and 2000, while increasingly negative summer SMB is responsible for substantial mass losses after 2000. Spatial variations in SMB between glaciers or regions of the ice cap are likely associated with local topography and its effect on orographic precipitation. We advocate for models to leverage the growing availability of observational resources to improve SMB predictions. We demonstrate an approach that incorporates complementary datasets, while addressing their inherent uncertainties, to constrain models and provide robust estimates of spatiotemporal SMB and associated uncertainties.

1. Introduction

Jostedalsbreen Ice Cap is the largest glacier in mainland Europe and constitutes ~20% of the glacierized area in Norway (Andreassen and others, 2022). Situated in a sparsely populated area in western Norway, the ice cap is a major tourist attraction, stimulating local business and supporting livelihoods, in addition to providing meltwater runoff for hydropower production, agriculture and ecosystems. Jostedalsbreen is a complex glacier system divided into more than 80 units (Andreassen and others, 2022), some of which have been monitored through glaciological, geodetic or front position surveys over shorter or longer time periods during the past century (e.g. Winkler, 1996; Andreassen and others, 2020, 2023; Kjøllmoen and others, 2022).

Owing to the maritime climate in the region, with relatively mild summers and precipitation-rich winters, glaciers of Jostedalsbreen experience substantial mass-turnover and are sensitive to climate fluctuations (Oerlemans, 1992; Nesje and others, 2000; Winkler and others, 2009). The most notable example is a period of mass gain during the 1990s, documented in the long-term glaciological mass-balance records of Nigardsbreen and Austdalsbreen, and the subsequent advances of several outlet glaciers (e.g. Andreassen and others, 2005, 2020; Winkler and others, 2009; Kjøllmoen and others, 2022). Since the early 2000s, glaciological mass balance, derived from interpolation of stake measurements on the glacier surface, and front position measurements indicate significant mass loss and retreat, although years with mass surplus are still registered, for example, 2012, 2020 (Andreassen and others, 2020; Kjøllmoen and others, 2022).

Current and future climate change is expected to accelerate glacier mass loss and retreat, which in turn may alter runoff regimes of glacierized catchments in Norway (e.g. Nesje and others, 2008; Giesen and Oerlemans, 2010; Engelhardt and others, 2015; Hanssen-Bauer and others, 2017; Compagno and others, 2021; Nesje, 2023). Future retreat of major outlet glaciers of Jostedalsbreen or complete disintegration of the ice cap would have strong ecological and economic implications in the region, and could increase the risk of glacier and/or paraglacial hazards (Jackson and Ragulina, 2014; Haeberli and Whiteman, 2021). Understanding the response of glaciers and ice caps to climate change requires knowledge about mass changes in space and time. For Jostedalsbreen, mass changes from existing observational records are difficult to reconcile, as these only provide temporal and spatial snapshots. For example, long-term glaciological mass-balance records only exist for two glaciers (Nigardsbreen and Austdalsbreen; Kjøllmoen and others, 2022), and meaningful geodetic mass-balance estimates are available only at multi-year intervals (e.g. Hugonnet and others, 2021; Andreassen and others, 2023) and with incomplete spatial coverage (Andreassen and others, 2020, 2023). Modelling studies on Jostedalsbreen are also limited to selected outlet glaciers (e.g. Oerlemans, 1997; Laumann and Nesje, 2009; Engelhardt and others, 2014; Li and others, 2015; Trachsel and Nesje, 2015; Sjursen and others, 2023).

Glacier mass-balance models are valuable tools to investigate glacier mass changes and provide surface mass balance (SMB) estimates with complete temporal and spatial coverage. Common modelling approaches range from physics-based energy-balance approaches (e.g. Andreassen and Oerlemans, 2009; Giesen and Oerlemans, 2010; Zolles and others, 2019; Eidhammer and others, 2021) to relatively simple temperature-index models (e.g. Schuler and others, 2005; Engelhardt and others, 2014; Huss and Hock, 2015; Geck and others, 2021) that parameterize the relationship between temperature and melt (see e.g. Hock, 2005; Zekollari and others, 2022). While energy-balance models provide complete representation of the underlying physical processes, they often suffer from a lack of detailed in situ meteorological data and coarse-resolution climate model data (Réveillet and others, 2018). Therefore, simpler, less input-demanding temperature-index models are often preferred, as they only require temperature and precipitation as meteorological input, both of which are more readily available in many areas of the world. However, the performance of temperature-index models heavily depends on model parameter values and their calibration to site-specific mass-balance observations (Schuster and others, 2023).

In situ observations of glacier mass change, for example, through the glaciological method, are sparse and concentrated to a handful of well-monitored regions (WGMS, 2024). This challenge is only recently starting to become alleviated by increased spatial coverage of geodetic mass balance derived from satellite-sensed surface elevation changes (e.g. Dussaillant and others, 2019; Shean and others, 2020; Hugonnet and others, 2021). Consequently, satellite-borne geodetic mass balances are increasingly used to constrain temperature-index model parameters (e.g. Rounce and others, 2020a, 2023; Compagno and others, 2021). However, these observations represent multi-year signals of mass change, integrated over the glacier area, and are afflicted with relatively large uncertainties. Therefore, these multi-year geodetic data provide only coarsely-resolved spatiotemporal variability and limited constraints on model parameters (Sjursen and others, 2023), such that mass-balance models still require seasonal mass-balance signals to adequately constrain accumulation and ablation.

The goal of this work is to provide a reconstruction of the SMB of Jostedalsbreen Ice Cap from 1960 to 2020 that is coherent in space and time, and in line with several observational datasets. We aim to capture the spatiotemporal variability in SMB in detail over seasonal timescales and to assess SMB variability in response to potential climatic and topographic drivers. To achieve this, we model the SMB of Jostedalsbreen using a temperature-index model forced by high-resolution (1 km) daily mean temperature and daily total precipitation from the seNorge_2018 dataset (Lussana and others, 2019). To ensure that the modelled SMB conforms with observational datasets we employ a Bayesian framework to estimate model parameters in a two-step procedure: (1) seasonal glaciological mass-balance measurements are used to estimate a global model parameter set that constrains accumulation and ablation on the ice cap, and (2) the global parameter set is employed as prior to estimate glacier-specific parameter sets using geodetic mass-balance observations for each glacier. The procedure allows us to quantify the uncertainty in simulated SMB that arises from uncertainty in model parameters and other sources such as uncertainties that arise from limitations in the model structure and input data. We thus demonstrate an approach that leverages observational datasets with complementary mass-balance information to provide robust spatiotemporal SMB estimates across a diverse region.

2. Study area and data

2.1 Study area

Jostedalsbreen stretches along a southwest to northeast axis in western Norway (Fig. 1), covering an area of 458 km² in 2019 (Andreassen and others, 2022) and with an estimated volume of 70.6 km³ (~2020; Gillespie and others, 2024). In the latest glacier inventory (2019) the ice cap is divided into 81 glaciers ranging in area from <0.1 km² to 46.2 km² (Tunsbergdalsbreen; Andreassen and others, 2022). In previous inventories (1966: Winsvold and others, 2014, 2006: Andreassen and others, 2012) the ice cap was divided into 82 units, but two disconnected units were removed and one was added in the 2019 inventory (Andreassen and others, 2023). In this work we consider the 82 glaciers from the 1966 and 2006 inventories.

Figure 1. Overview of Jostedalsbreen Ice Cap in western Norway with glacier outlines from 1966 (Winsvold and others, 2014), 2006 (Andreassen and others, 2012) and 2019 (Andreassen and others, 2022). Hatched areas show glaciers with glaciological SMB observations. The coordinate systems are geographical coordinates in the inset and UTM33N, datum ETRS89 in main map.

Jostedalsbreen consists of three main parts that are connected by relatively narrow bands of ice (Gillespie and others, 2024). The three parts are hereafter referred to as North (northeast of Lodalsbreen, Norwegian Water Resources and Energy Directorate (NVE) glacier ID2266; 24 glaciers), South (south of Grensevarden/ID2332; 19 glaciers) and Central (39 glaciers; Fig. 1) following Gillespie and others (2024), see list of IDs in Appendix A. Surface elevation extends from below 400 m a.s.l. at the tongues of the largest outlet glaciers to above 1900 m a.s.l. on the ice-cap plateau. The northwestern side of the ice cap is characterized by mainly short and steep glaciers, while the southeastern side holds several large valley glaciers, of which Nigardsbreen is most extensively studied (e.g. Østrem and others, 1976; Oerlemans, 1997; Engelhardt and others, 2014; Li and others, 2015; Gjerde and others, 2023).

2.2 Mass-balance data and glacier inventory

Observations of mass change on Jostedalsbreen during the past 60 years are available from several sources with different spatial coverage and temporal resolution. Glaciological SMB measurements have been performed by NVE on parts of the ice cap since 1962 (Table 1; e.g. Kjøllmoen and others, 2022). The outlet glaciers Nigardsbreen and Austdalsbreen have continuous long-term SMB records since the mass-balance years 1962 and 1988, respectively, providing observations of annual and seasonal (winter and summer) SMB over 59 (1962–2020) and 33 (1988–2020) consecutive mass-balance years over the study period (Kjøllmoen and others, 2022). Note that when referring to mass-balance years, we refer to the year marking the end of a mass-balance year (e.g. 1962 refers to the mass-balance year 1961/1962). Three other glaciers have shorter SMB records: Tunsbergdalsbreen (1966–1972; 7 years), Vesledalsbreen (1967–1972; 6 years) and Supphellebreen (1964–1967, 1973–1975, 1979–1982; 11 years; Kjøllmoen, 2017). Mass loss due to lake-calving is generally negligible for Jostedalsbreen, except for Austdalsbreen where calving is estimated annually by NVE and reported to account for up to 20% of the annual ablation (e.g. Kjøllmoen and others, 2022).

Of the five original glaciological SMB records, four have been homogenized: Nigardsbreen, Austdalsbreen, Tunsbergdalsbreen and Vesledalsbreen (Andreassen and others, 2016; Kjøllmoen, 2017, 2022) and Nigardsbreen has been partly calibrated due to significant differences between geodetic and glaciological mass-balance records (Andreassen and others, 2016; Kjøllmoen, 2022). Supphellebreen has not been homogenized due to lack of data (Kjøllmoen, 2017). In this study, we consider the homogenized and calibrated records for all glaciers except for Supphellebreen, where we use the original record. For Nigardsbreen, glaciological SMB is measured for the basin consisting of Nigardsbreen (ID2297) and two smaller glaciers (ID2299 and 2311). In this study, we employ the same basin as the glaciological record when comparing modelled and glaciological SMB for Nigardsbreen.

Geodetic mass-balance estimates are available for Jostedalsbreen, or parts of the ice cap, for time periods of various length. Satellite-borne geodetic mass balance from repeat ASTER Digital Elevation Models (DEMs) is available for all glaciers for the period 2000–2019 (Hugonnet and others, 2021). In addition, geodetic mass balance is available for 49 glaciers for the period 1966–2020 (Andreassen and others, 2023). The latter estimates are based on aerial photographs from 1966 and airborne LiDAR surveys in 2020, and cover central and northern parts of the ice cap. Geodetic mass balances for Nigardsbreen, Tunsbergdalsbreen and Austdalsbreen are also available for other periods (e.g. Andreassen and others, 2016, 2020, 2023), but are not included in this study since they provide limited additional information.

Glacier outlines required for modelling SMB of Jostedalsbreen and individual glaciers are available from 1966 (Winsvold and others, 2014), 2006 (Andreassen and others, 2012) and 2019 inventories (Andreassen and others, 2022). Ice divides in the 2019 inventory are mostly aligned with 2006 outlines but updated for some glaciers (Nigardsbreen, Austdalsbreen and neighbouring glaciers) to harmonize with those used in glaciological SMB calculations (Andreassen and others, 2022; Kjøllmoen, 2022). The 1966 ice divides have been homogenized with the 2019 inventory (Andreassen and others, 2023). Additional glacier outlines used in calculation of glaciological SMB are available for Nigardsbreen (1964, 1974, 1984, 2009, 2013, 2020) and Austdalsbreen (1988, 2009) (e.g. Kjøllmoen and others, 2022).

2.3 Meteorological forcing data

As meteorological forcing for the SMB model we employ gridded 1 km resolution daily mean temperature and daily total precipitation from seNorge_2018 version 21.09 (Lussana and others, 2019; Lussana, 2021). The seNorge (https://www.senorge.no/) collection of datasets is provided by the Norwegian Meteorological Institute (MET Norway) and are based on spatial interpolation of measurements from a large network of weather stations across the Norwegian mainland, while also leveraging monthly precipitation reference fields from 3 km climate model simulations from HARMONIE to improve precipitation estimates in data-sparse regions (Lussana and others, 2019; Lussana, 2020). Several versions of seNorge (e.g. seNorge1.1: Mohr, 2008 and seNorge:_2018 Lussana and others, 2019) have previously been applied in SMB and runoff modelling of glacierized areas in Norway (e.g. Engelhardt and others, 2013, 2014; Li and others, 2015; Sjursen and others, 2023), and to correct downscaled climate model projections in assessments of climate-change impacts (Wong and others, 2016; Hanssen-Bauer and others, 2017). Overall, seNorge_2018 is considered to improve precipitation estimates compared to its predecessors. Nevertheless, the probability of large errors is considered greatest for precipitation in remote, mountainous regions with low station density (Lussana and others, 2019), where glaciers commonly reside.

3. Methods

3.1 SMB model

The SMB of a glacier over a given period (e.g. year, season) is the sum of accumulation and ablation on its surface (Cogley and others, 2011). Accumulation at Jostedalsbreen is mainly in the form of snowfall, while surface ablation is mainly melt of snow, firn and ice. We calculate the SMB of Jostedalsbreen on the 1 km resolution DEM of the seNorge dataset. We use the temperature-index model (see e.g. Hock, 2005) employed in Sjursen and others (2023), where melt of snow or ice in a grid cell i at the daily time step t, m _snow/ice,i,t (mm w.e. $^\circ$C⁻¹ d⁻¹), is computed using melt factors for snow and ice, MF _snow/ice, when the mean daily temperature in a gridcell T _i,t is above a melt threshold temperature ($T_{\rm m} = 0^\circ$C):

(1)$$m_{\rm{snow/ice}, i, t} = \left\{\matrix{ MF_{\rm{snow/ice}} ( T_{i, t} - T_{\rm m}) \hfill & \rm{if } \quad T_{i, t} > T_{\rm m},\; \hfill \cr 0 \hfill & \rm{if } \quad T_{i, t} \leq T_{\rm m}. \hfill }\right.$$

Firn melt is estimated as the average of daily melt of snow and ice, since the albedo of firn is typically between that of snow and ice (Cuffey and Paterson, 2010). To account for differences in albedo between snow and ice, we set MF _ice = MF _snow/0.7. Thereby we assume that the melt rate for clean snow is 70% of that of clean ice (e.g. Singh and others, 2000, in line with calibrated values for all glaciers in Norway from Engelhardt and others, 2013). Daily accumulation in a grid cell is computed as the fraction of the daily total precipitation in the cell falling as snow, assuming a linear decrease from entirely solid to liquid in a $\pm 1^\circ$C interval around 1$^\circ$C (Jennings and others, 2018). We evaluate mass changes of the ice cap over a hydrological year (1 October–30 September, with 30 April as end of accumulation season). In model calibration and validation we assess mass changes based on dates of maximum and minimum mass for a more accurate comparison to available observations of individual glaciers (e.g. glaciological SMB is measured for end of accumulation/melt seasons).

Due to the uncertainty in the meteorological forcing data, we add a temperature correction T _corr ($^\circ$C) to the daily mean temperature and multiplying the daily total precipitation by a precipitation correction factor P _corr (dimensionless). The unknown model parameters are thus MF _snow, P _corr and T _corr, whose values are constrained using the Bayesian framework described in Section 3.2.

3.2 Bayesian parameter estimation

We employ a Bayesian framework (see e.g. Gelman and others, 2014) to estimate probability distributions of the SMB model parameters and to quantify uncertainty in modelled SMB. Our procedure consists of two steps that leverage two different observational datasets and aims to estimate model parameters that constrain accumulation and ablation on Jostedalsbreen (step 1), while also providing accurate estimates of SMB for each individual glacier (step 2). These steps are first summarized below, before we describe the details of each step in Sections 3.2.1 and 3.2.2.

In step 1 we estimate a global parameter set θ = {P _corr,glob,  T _corr,glob,  MF _snow,glob} that most accurately represents the SMB of the entire ice cap, constrained by seasonal glaciological SMB observations from five glaciers (Table 1). We also estimate the distribution of a model error that is not accounted for by the model parameters and represents the structural model uncertainty (e.g. due to missing or simplified process representation) and uncertainty in the input data. This allows us to properly quantify the predictive uncertainty of the model and provide robust SMB estimates.

Table 1. Overview of glaciological SMB observations for glaciers of Jostedalsbreen used in this study

ID refers to the Norwegian Water Resources and Energy Directorate (NVE) glacier ID (Andreassen and others, 2022). Location refers to glacier location in the South (S), Central (C) or North (N) part of the ice cap. Area and elevation range (min–max) refer to the 2019 inventory (Andreassen and others, 2022). n _ann and n _seas are the respective number of mass-balance years with annual and seasonal glaciological SMB observations for each glacier over the study period 1960–2020 (e.g. Kjøllmoen, 2017; Andreassen and others, 2020; Kjøllmoen and others, 2022). Period refers to the time period covered by SMB observations (mass-balance years). Glaciological SMB measurements for Nigardsbreen includes ID2297, 2299 and 2311.

In step 2 we spatially adjust the precipitation and temperature correction parameters recovered in step 1, by estimating a set of glacier-specific precipitation and temperature correction parameters ϕ _j = {P _corr,j,  T _corr,j} for each glacier j. To this end, we employ two decadal geodetic mass-balance observations for each glacier (2000–2009 and 2010–2019) from Hugonnet and others (2021) since this dataset covers the entire ice cap. To estimate ϕ _j, we use the posterior distributions of P _corr,glob and T _corr,glob obtained from step 1 as the prior distributions in step 2. The geodetic observations have low temporal resolution and high uncertainty, and posterior estimates can therefore be expected to be strongly influenced by the choice of prior distribution (Sjursen and others, 2023). We mitigate this by using posterior distributions from step 1 as priors in step 2, which represent reliable estimates of accumulation and ablation.

In step 2 we fix MF _snow to the median of the posterior distribution of MF _snow,glob and thus choose to spatially adjust two parameters that each mainly controls either accumulation (precipitation correction) or ablation (temperature correction) (Réveillet and others, 2017). This is because the geodetic observations provide limited information to constrain strongly correlated parameters (e.g. temperature correction and melt factor; Rounce and others, 2020b; Sjursen and others, 2023). Furthermore, we expect that spatial patterns of temperature and precipitation may not be accurately represented in seNorge_2018 due to the complex topography of the region and significant local effects on weather patterns. In step 2 we thus address possible spatial biases in the meteorological forcing data over the ice cap by adjusting the well-constrained parameter values in step 1.

3.2.1 Step 1: estimation of global parameter set and model error

We formulate a deterministic model (Eqn (2)) that is similar to those of Rounce and others (2020b) and Sjursen and others (2023), but that also takes into account that the SMB model is an imperfect representation of an observed system. That is, instead of assuming that the model describes the observed system up to an observation error $\epsilon _{n}$, we include an additional unknown model error η _n:

(2)$$B_{{\rm obs}, n} = B_{{\rm mod}, n}( X_n,\; \theta) + \epsilon_{n} + \eta_{n} ,\; $$

where B _obs,n and B _mod,n(X _n,  θ) are observed and modelled SMB over n periods of mass change, respectively, and X _n is the set of model input data. Here, η _n is meant to represent any predictive uncertainty that is not accounted for by parameter uncertainty. This includes uncertainty in the model structure, for example, from missing or crudely parameterized physical processes, but also other sources of uncertainty that are not accounted for otherwise. We consider $\epsilon _{n}$ and η _n to be statistically independent since there is no physical relation between these errors. Further, we assume $\epsilon _{n}$ and η _n to be normally distributed (${\cal N}$) with means of zero and constant variances. The variance of the distribution of $\epsilon _{n}$ is given by the uncertainty in the SMB observation $\sigma _{B_{\rm obs}}^2$, while η _n has unknown variance $\sigma _{\eta }^{2}$:

(3)$$\epsilon_{n} \sim {\cal N}( 0,\; \, \sigma_{B_{\rm obs}}^{2}) ,\; $$

(4)$$\eta_{n} \sim {\cal N}( 0,\; \, \sigma_{\eta}^{2}) .$$

We employ Markov chain Monte Carlo (MCMC) simulations that make use of the following proportionality in Bayes’ theorem to estimate the joint posterior distribution of θ and σ _η given a set of mass-balance observations B _{obs,1 : N} and input data X _{1 : N}:

(5)$$p( \theta,\; \sigma_{\eta}\vert B_{{\rm obs}, 1\, \colon \, N},\; X_{1\, \colon \, N}) \, \propto \, p( \theta) \, p( \sigma_{\eta}) \, L_{B_{\rm obs}} ,\; $$

where p denotes probability and $L_{B_{\rm obs}} = p( B_{{\rm obs}, 1\, \colon \, N}\vert \theta ,\; \, \sigma _{\eta },\; \, X_{1\, \colon \, N})$ is termed the likelihood: the probability of observing the data B _{obs,1 : N} given our deterministic model (Eqn (2)). In MCMC simulations we employ the logarithm of the likelihood function $l_{B_{\rm obs}} = \ln ( L_{B_{\rm obs}})$ to ensure stability and efficient computation. Under the assumption of independent and normally distributed errors with constant variances, we formulate the log-likelihood as follows:

(6)$$\eqalign{l_{B_{\rm obs}} & = -{N\over 2} \ln ( 2 \pi) -{N\over 2} \ln( \sigma_{B_{\rm obs}}^2 + \sigma_{\eta}^2) \cr & \quad- {1\over 2 ( \sigma_{B_{\rm obs}}^2 + \sigma_{\eta}^2) }\sum\limits_{n = 1}^{N} ( B_{{\rm obs}, n} - B_{{\rm mod}, n}( X_n,\; \theta) ) ^2 .}$$

We employ seasonal SMB observations such that the modelled SMB B _mod,n over the period n is the modelled summer or winter SMB, and B _obs,n the SMB observation for the same period, with associated uncertainty $\sigma _{B_{\rm obs}}$. We assume that seasonal SMB observations are conditionally independent given our model, such that we can express the full log-likelihood function $l_{B_{\rm seas}}$ as the sum of the log-likelihood functions for each of winter and summer SMB ($l_{B_{\rm w}}$ and $l_{B_{\rm s}}$, respectively):

(7)$$l_{B_{\rm seas}} = l_{B_{\rm w}} + l_{B_{\rm s}} ,\; $$

where $l_{B_{\rm w}}$ and $l_{B_{\rm s}}$ are given by Equation (6).

For estimation of the global model parameter set θ = {P _corr,glob,  T _corr,glob,  MF _snow,glob} and the standard deviation (SD) in model error σ _η, we employ seasonal glacier-wide glaciological SMB observations (Table 1) for every other mass-balance year of the period 1962–2020 (even years, starting with mass-balance year 1962 and ending with 2020), totalling 56 mass-balance years of winter and summer SMB for the five glaciers. We use annual and seasonal SMB observations for the remaining 53 mass-balance years for validation of posterior predictive SMB. When comparing modelled SMB to observations from Nigardsbreen and Austdalsbreen, we employ the same time series of glacier outlines as used in glaciological records (for shorter records the 1966 outline is considered representative).

Following Sjursen and others (2023), we determine the uncertainty in seasonal glacier-wide glaciological SMB based on estimates from the reanalysis of the long-term glaciological SMB records in Norway (Andreassen and others, 2016), assuming that observations of summer and winter SMB are independent (Dyurgerov and Meier, 1999) and that the uncertainty in summer SMB accounts for two-thirds of the uncertainty of glacier-wide annual SMB (e.g. Kjøllmoen and others, 2022). Uncertainty in glacier-wide annual SMB is estimated at ±0.34 and ±0.30 m w.e. a⁻¹ for Nigardsbreen (1964–2013) and Austdalsbreen (1988–2009), respectively (Andreassen and others, 2016). Individual error estimates are lacking for the short-term glaciological SMB records on Jostedalsbreen, but are considered to be of similar magnitude as the long-term series (Kjøllmoen, 2017). For simplicity, we assume that the estimated uncertainty in glacier-wide annual SMB measurements for Nigardsbreen is representative for all glaciers; $\sigma _{B_{\rm w}} = 0.19$ m w.e. a⁻¹ and $\sigma _{B_{\rm s}} = 0.28$ m w.e. a⁻¹ for winter and summer SMB, respectively (Sjursen and others, 2023).

3.2.2 Step 2: estimation of glacier-specific precipitation and temperature correction

Our deterministic model for step 2 is similar to step 1 (Eqn (2)), but applied to each glacier individually with $\epsilon _{n, j} \sim {\cal N}( 0,\; \, \sigma _{B_{{\rm obs}, n, j}}^{2})$. For the model error we assign a decadal model uncertainty σ _η,10 yr from the posterior mean of σ _η, given that σ _η is the model uncertainty associated with a seasonal SMB prediction. The log-likelihood function for each glacier j in step 2 is thus:

(8)$$l_{B_{{\rm obs}, 10\, {\rm yr}, j}} = l_{B_{{\rm obs}, 00-09}, j} + l_{B_{{\rm obs}, 10-19}, j} ,\; $$

where $l_{B_{{\rm obs}, 10\, {\rm yr}, j}}$ is the combined log-likelihood for the geodetic mass-balance observations of glacier j over each of the periods 2000–09 (N = 00–09) and 2010–19 (N = 10–19), and is given by

(9)$$\eqalign{l_{B_{{\rm obs}, N, j}} & = -{1\over 2} \ln ( 2 \pi) -{1\over 2} \ln( \sigma_{B_{{\rm obs}, N, j}}^2 + \sigma_{\eta, 10\, {\rm yr}}^2) \cr & \quad- {1\over 2 ( \sigma_{B_{{\rm obs}, N, j}}^2 + \sigma_{\eta, 10\, {\rm yr}}^2) }( B_{{\rm obs}, N, j} - B_{{\rm mod}, N, j}( X_j,\; \phi_{\,j}) ) ^2 .}$$

For each decadal geodetic mass-balance observation we assign the uncertainty $\sigma _{B_{{\rm obs}}, \, N, \, j}$ reported by Hugonnet and others (2021) for a given period N and glacier j. Similar to step 1, we employ MCMC simulations to estimate the posterior of ϕ _j = {P _corr,j,  T _corr,j} for each glacier (see Appendix B for details).

3.2.3 Posterior predictive SMB simulations

We perform posterior predictive SMB simulations using posterior distributions of step 1, but with posterior means of global parameters P _corr,glob and T _corr,glob corrected to posterior means of P _corr,j and T _corr,j estimated in step 2. This allows us to run posterior predictive simulations for the deterministic model described in Equation (2) with posterior estimates of MF _snow,glob and η _n, but with spatial adjustment of posteriors of P _corr,glob and T _corr,glob. More specifically, for each glacier we adjust the posterior of the global precipitation correction P _corr,glob by $c_j = \mu _{P_{\rm corr}, \, j}-\mu _{P_{\rm corr}, \, {\rm glob}}$, the difference between the means of P _corr,j and P _corr,glob, respectively:

(10)$$P_{\rm corr, glob} + c_j \sim {\cal N}( \mu_{P_{\rm corr, glob}} + c_j,\; \, \sigma_{P_{\rm corr, glob}}^{2}) .$$

We perform a corresponding operation for temperature correction using T _corr,glob and T _corr,j. The underlying assumption in Equation (10) is that posterior distributions of P _corr,glob and T _corr,glob are approximately normal, which is demonstrated in Section 4.1.

In posterior predictive simulations we employ the set of outlines from 1966, 2006 and 2019, and follow the principle applied in homogenization of glaciological SMB records in Norway of using each outline for half of the period before and after its date (Andreassen and others, 2016).

3.2.4 Prior distributions for global parameter set

As prior distribution for P _corr,glob we choose a normal distribution truncated at zero (to ensure positivity), with a mean of 1.0 and SD of 0.25 as the prior distribution. Although previous evaluation reveals that precipitation estimates over outlet glaciers of Jostedalsbreen (Nigardsbreen and Austdalsbreen) may be underestimated in seNorge_2018 (Sjursen and others, 2023), we do not know if estimates of precipitation correction factors (based on observations from 1990 to 2009 for Austdalsbreen and Nigardsbreen) are representative for the whole ice cap over the period 1960–2020. Our choice of prior for P _corr,glob gives 95% confidence interval limits at ~0.5 and 1.5, meaning that we are confident that the under- or overestimation of precipitation sums over Jostedalsbreen do not exceed 50%. Since there are no indications of bias or large errors in daily mean temperature in seNorge_2018 (Lussana and others, 2019), we choose a normal distribution with mean ± SD of $0\pm 0.5^\circ$C (95% confidence interval limits at ~$\pm 1^\circ$C) for the prior distribution of T _corr,glob.

Similarly to Rounce and others (2020b) and Sjursen and others (2023) our choice of prior for MF _snow,glob is based on Braithwaite (2008), who found a value of 4.1 ± 1.5 mm w.e. $^\circ$C⁻¹ d⁻¹ for the melt factor for snow at the equilibrium line altitude (ELA) of 66 glaciers. However, in light of previous parameter estimates for outlet glaciers of Jostedalsbreen (Sjursen and others, 2023), we believe that values are more likely closer to the mean and therefore adopt a zero-truncated normal distribution with mean ± SD of 4.1 ± 1.0 mm w.e. $^\circ$C⁻¹ d⁻¹ (95% confidence interval at 2.1 and 6.1 mm w.e. $^\circ$C⁻¹ d⁻¹) as prior distribution for MF _snow,glob.

As the prior for the SD of the model error σ _η we choose a half-normal distribution since the SD is a positive number and because we believe that there is a high probability of small errors and a low probability of very large errors. Further, we choose a scale parameter of 0.67 for the half-normal distribution such that the model error is likely (95% confidence interval) within 1.5 m w.e., which reflects the reported error distribution in studies with similar SMB model set-up (Huss and Hock, 2015).

4. Results

In this section we show the posterior parameter distributions resulting from our two-step parameter estimation procedure (Section 4.1) and present the simulated SMB of Jostedalsbreen from 1960 to 2020, highlighting both the spatial and temporal variability (Section 4.2). Unless specified otherwise, reported SMB is based on simulations with 1000 posterior predictive samples of the posterior distribution shown in Section 4.1.

4.1 Posterior parameter estimates

We find that posterior distributions of θ = {P _corr,glob,  T _corr,glob,  MF _snow,glob} and σ _η in step 1 of parameter estimation are well constrained (Fig. 2). The posterior of the precipitation bias correction (P _corr,glob) has a mean/median ± SD of 1.25/1.25 ± 0.04, and naturally shows the lowest spread since it has limited correlation to other parameters and is informed by winter SMB observations. The corresponding statistics for the posterior of the melt factor for snow (MF _snow,glob) is 3.58/3.56 ± 0.25 mm w.e. $^\circ$C⁻¹ d⁻¹ and $-0.14/-0.14\pm 0.34^\circ$C for the temperature bias correction (T _corr,glob), both shifted towards slightly lower-ablation values compared to prior distributions.

Figure 2. Marginal prior (grey dashed lines) and posterior (blue solid lines) probability distributions of global parameter set: (a) precipitation correction factor P _corr,glob, (b) melt factor for snow MF _snow,glob, and (c) temperature bias correction T _corr,glob and (d) SD in model error σ _η.

The mean/median ± SD of the posterior distribution of model uncertainty (i.e. SD in model error, σ _η) is 0.32/0.32 ± 0.04 m w.e. The model uncertainty reflects the error in modelled glacier-wide seasonal SMB (Eqn (2)), and is slightly higher than the uncertainty in observed glacier-wide winter and summer SMB from glaciological records. The error in modelled glacier-wide annual SMB can be estimated as the sum of normally distributed errors (following our assumption in Eqn (4)), such that $\eta _{B_{{\rm mod}, a}} \sim {\cal N}( 0,\; \, 0.45^{2})$ m w.e.

Overall, posterior distributions of spatially corrected P _corr,j from step 2 display lower values (lower precipitation sums) compared to P _corr,glob, while posteriors of T _corr,j are mostly shifted towards higher values (higher temperature) compared to T _corr,glob (Fig. 3). The minimum/maximum values of the medians of the posteriors of P _corr,j and T _corr,j are 1.20/1.27 and $-1.06/0.98^\circ$C, respectively. Posteriors show spatial patterns across Jostedalsbreen, with higher values of P _corr,j (higher precipitation sums) on the southeastern side and in northern parts of the ice cap (Fig. 3a) and higher values of T _corr,j on the northwestern side and in the south (Fig. 3b). However, there are some local variations to these patterns, for example, the smaller glaciers in the central northwestern region that shows high values of P _corr,j and low values of T _corr,j.

Figure 3. Median values of marginal posterior probability distributions of (a) P _corr,j and (b) T _corr,j for each glacier j of Jostedalsbreen.

4.2 Mass balance of Jostedalsbreen 1960–2020

Overall, the modelled SMB of Jostedalsbreen is slightly negative over the period 1960–2020 (Fig. 4). We find a median cumulative SMB of −4.05 m w.e. (95% credible interval (CI): −12.52, 5.12 m w.e.), equivalent to an annual SMB rate of −0.07 m w.e. a⁻¹ (95% CI: −0.21, 0.08 m w.e. a⁻¹). The median summer and winter SMB rates over the model period are −2.10 m w.e. a⁻¹ (95% CI: −2.19, −2.00 m w.e. a⁻¹) and 2.02 m w.e. a⁻¹ (95% CI: 1.92, 2.14 m w.e. a⁻¹), respectively. Considering individual glaciers over the period 1960–2020, annual SMB rates are generally slightly positive for glaciers in the southwestern part of Jostedalsbreen, close to zero for glaciers in the central part and overall negative for glaciers in the northeast (Fig. 5a). Some smaller glaciers at the margins of the central and southwestern parts of the ice cap also display negative SMB. The largest outlet glaciers Tunsbergdalsbreen and Nigardsbreen show negative annual SMB rates of −0.43 m w.e. a⁻¹ (95% CI: −0.62, −0.23 m w.e. a⁻¹) and −0.09 m w.e. a⁻¹ (95% CI: −0.26, −0.12 m w.e. a⁻¹), respectively.

Figure 4. (a) Median glacier-wide annual (grey, whiskers represent 95% CI), winter (blue) and summer (red) SMB (m w.e.) of Jostedalsbreen over the period 1960–2020, based on 1000 posterior predictive samples. (b) Cumulative SMB for the ice cap from 1960 to 2020, based on median of 1000 posterior predictive samples (shaded area represents 95% CI).

Figure 5. (a) Glacier-wide annual average SMBs (median SMB in m w.e. a⁻¹) using 1000 posterior predictive samples and gridded (b) annual, (c) winter and (d) summer SMB rates over the period 1960–2020 based on median parameter values. Glaciers with glaciological SMB records are highlighted (Sup: Supphellebreen, Tun: Tunsbergdalsbreen, Nig: Nigardsbreen, Aus: Austdalsbreen, Ves: Vesledalsbreen).

Our model results reveal that Jostedalsbreen has experienced both periods of mass loss and gain over the past 60 years, with large temporal variability (Fig. 4 and Table 2). In the 1960s, the ice cap experienced significant mass loss, followed by a relatively stable period from the 1970s until the mid-1980s. From the end of the 1980s and throughout the 1990s the ice cap gained mass, followed by a mass deficit of similar magnitude from 2000 until the mid-2010s. From the mid-2010s, the mass of Jostedalsbreen has again been relatively stable according to our model results. Considering magnitudes of decadal variations in SMB of the ice cap (Table 2), the 1980s and 1990s are the only positive decades, with the 1990s showing the largest mass gain. The positive SMB over these two decades is driven by relatively low magnitude of summer SMB (85 and 92% of the average for 1960–2020, respectively) and higher-than-average magnitude of winter SMB in the 1990s (130%). The most negative decade in terms of annual SMB rate over the ice cap is the 1960s, followed by the most recent decades 2000–2009 and 2010–2019. The 1960s display average summer SMB, but very low winter SMB over the ice cap (73%). In contrast, the 2000s and 2010s show average winter SMB rates, with overall negative annual rates dominated by high magnitudes of summer SMB (118 and 109%, respectively). However, within the past decade there are relatively large interannual variations in SMB, with relatively high magnitudes both for positive and negative years (Fig. 4a).

Table 2. Overview of modelled annual and seasonal SMB rates for Jostedalsbreen for different decades

Values in m w.e. a⁻¹ given as: rate (95% CI limits). Percentages are relative to median rate 1960–2020.

Within the overall temporal trends there is significant variability in SMB between regions of the ice cap (Figs 5 and 6). We evaluate these trends on a decadal basis and for the regions North, Central and South (Fig. 1). The 1960s display negative annual SMB over most of the ice cap, with particularly negative rates in the North (Fig. 6b). While the 1970s and 1980s indicate near-balanced or positive rates for glaciers in the South and Central parts, SMB rates in the North remain negative. In the 1990s, SMB rates are overall positive for all three parts of the ice cap, with highest magnitude in the South. All regions display negative SMB rates in the 2000s and 2010s. The South and Central parts show similar annual SMB rates to the 1960s, but SMB rates for glaciers in the North parts are strikingly more negative in the 1960s and 1970s compared to the 2000s.

Figure 6. Distribution of (a) winter, (b) annual and (c) summer SMB (m w.e.) for different decades and regions in order: North (N), Central (C) and South (S). Horizontal lines in boxplots indicate SMB rate (mean) and boxes and whiskers extend to the IQR and minimum and maximum SMB, respectively.

Glaciers in the North and South of Jostedalsbreen generally display more negative summer SMB than glaciers in the Central part (Fig. 6c). The largest winter SMB rates are generally found for glaciers in the South of the ice cap (Figs 5a, 6a). Magnitudes of summer SMB show considerable temporal variability which is relatively uniform across regions (Fig. 6c). Differences in winter SMB, however, show both strong temporal and spatial variability between decades, with particularly large variability in the North and South (Fig. 6a).

5. Discussion

5.1 Meteorological drivers of temporal variability in SMB

Our results indicate that temporal trends in SMB on Jostedalsbreen were largely driven by winter accumulation variability between 1960 and 2000, while increasingly negative summer ablation dominates annual SMB after 2000 (Table 2, Figs 4 and 6). The mass gain during the late 1980s through the 1990s is in line with glaciological SMB records for Nigardsbreen and Austdalsbreen and subsequent advances of several outlet glaciers (e.g. Bergsetbreen, Bødalsbreen, Brenndalsbreen, Kjenndalsbreen, Nigardsbreen; NVE, 2022). This period of mass gain is found for glaciers in western Norway in general (e.g. Andreassen and others, 2005) and has been attributed to increased snow accumulation (Andreassen and others, 2005; Winkler and others, 2009) associated with transient changes in large scale weather patterns. In particular, this period has been shown to coincide with a period of strongly positive North Atlantic Oscillation (NAO) index (e.g. Nesje and others, 2000; Marzeion and Nesje, 2012; Trachsel and Nesje, 2015; Mutz and others, 2016), which is characterized by strong westerly winds and anomalously high winter precipitation over western Scandinavia. The magnitude of our modelled 1990s winter SMB anomaly (Table 2) is also in line with weather station records which show that winter precipitation (December–February) in western Norway was particularly high in the 1990s (Hanssen-Bauer, 2005; Konstali and Sorteberg, 2022), up to 30% higher than the 1900–2019 mean in the mid-1990s (Konstali and Sorteberg, 2022). While regional trends in temperature do not indicate particularly low summer temperatures in the 1980s and 1990s (Hanssen-Bauer, 2005), our results indicate that lower magnitudes of ablation have contributed to the overall mass gain of Jostedalsbreen around the 1990s. In contrast to the 1990s, our findings suggest that mass losses of the ice cap in the 1960s were primarily driven by lower than average winter SMB. This anomaly is also supported by weather station records from western Norway that show a significant negative winter (December–January) precipitation anomaly in the 1960s (~25% lower than the 1900–2019 mean in the mid-1960s; Konstali and Sorteberg, 2022).

Our results indicate that increasingly negative summer balances are driving increased mass loss of Jostedalsbreen since the early 2000s (Table 2, Figs 4 and 6), in line with previous findings for glaciers in Norway (e.g. Mutz and others, 2016). Overall, annual air temperature in western Norway has increased by ~0.7$^\circ$C from the start of the 20th century (1900–2014; Hanssen-Bauer and others, 2017) with the largest increases found in spring (March–May; 0.9$^\circ$C) and autumn (September–November; 0.8$^\circ$C). The modelled mass loss of Jostedalsbreen from around 2000 is in line with the overall negative trends in glaciological SMB records in Norway, with the 2000s as the most negative decade (Andreassen and others, 2020; Kjøllmoen and others, 2022). Our results show that strongly negative summer SMB rates from the early 2000s are to some degree counteracted by relatively high winter SMB rates (Table 2 and Figs 6a, c). Overall, precipitation in western Norway has increased by 18% between 1960 and 2019 (Konstali and Sorteberg, 2022), such that increased ablation due to higher temperatures may be compensated by increased winter accumulation. Andreassen and others (2020) found that the NAO index was positive for several years in the 2010s, and attributed part of the overall less negative SMB of glaciers in Norway in this decade to high winter precipitation rates. We find that summer SMB at Jostedalsbreen was less negative in the 2010s compared to the 2000s (Table 2), but with relatively large interannual variations in both summer and winter SMB (Fig. 4a). This increases the sensitivity of our decadal analysis to partitioning of years. Considering the decade 2011–2020 instead of 2010–2019 gives a positive SMB rate of 0.14 m w.e. a⁻¹ (compared to −0.29 m w.e. a⁻¹ for 2010–2019), with winter/summer SMB magnitudes 110/104% of the average of 1960–2020. Thus, controls on annual SMB in the 21st century seemingly vary between years. Large parts of Jostedalsbreen are located at high elevations; 72 (48)% of the area of the ice cap is located above 1500 (1600) m a.s.l. (2019 outline and DEM from 2020; the Norwegian Mapping Authority). This means that the ice cap currently has a relatively large accumulation area distributed over a small elevation range. However, with ongoing and future expected increases in temperatures and associated rising ELAs, as well as the feedback of surface lowering on SMB, it is unclear to what extent increased winter precipitation will continue to compensate for a shrinking accumulation area and stronger ablation.

Several studies (e.g. Nesje and others, 2000; Andreassen and others, 2005, 2020; Mernild and others, 2014; Trachsel and Nesje, 2015) have investigated the influence of winter and summer SMB on annual SMB for glaciers in different climatic settings in Norway and found the same overall relationship: the annual SMB of maritime glaciers and continental glaciers is mainly controlled by winter precipitation and summer temperatures, respectively. Following the approach of Andreassen and others (2005, 2020), we compared ratios of the SDs in winter and summer SMB to SDs in annual SMB (sBw/sBa and sBs/sBa, respectively) for each glacier of Jostedalsbreen. When computing sBw/sBa and sBs/sBa over the time series as a whole, ratios are relatively equal (e.g. for Nigardsbreen sBw/sBa and sBs/sBa is 0.69 and 0.52, respectively) and in line with Andreassen and others (2005). However, when evaluating sBw/sBa and sBs/sBa over 20 year rolling windows, ratios are not stationary (e.g. between 0.57 and 0.86 for sBw/sBa and 0.37 and 0.60 for sBs/sBa for Nigardsbreen; Fig. 11). This analysis indicates that the relative contribution of winter SMB to annual SMB was particularly high towards the end of the 20th century but is decreasing towards the present, along with a simultaneous increase in the relative importance of summer SMB (Fig. 11). These findings are in line with Trachsel and Nesje (2015) who found that for Scandinavian (including Nigardsbreen), variations in winter precipitation was more important than variations in summer temperature for annual SMB in the second half of the 20th century, but that the relative influence of summer temperature has increased in more recent years.

5.2 Variability in SMB between glaciers and regions

Modelled glacier-wide SMB rates on Jostedalsbreen show spatially varying signals with some distinct regional patterns (Fig. 5): overall slightly positive SMB in the south, near balance in the central part (but with overall negative SMB for large outlet glaciers), and relatively large negative SMB in the north. Following Andreassen and others (2023), we investigate topographic controls (statistics from 2019-inventory) on glacier-wide SMB rates over the period 1960–2020 and find the strongest correlation (−0.42, p ≤ 0.001) with hypsometric index (HI, calculated according to Jiskoot and others, 2009) and median elevation (0.39, p ≤ 0.001). The HI can be used to classify glaciers as very top heavy (HI < −1.5), top heavy ($-1.5\ < \ {\rm HI}\ < \ -1.2$), equidimensional (−1.2 < HI < 1.2), bottom heavy ($1.2\ < \ {\rm HI}\ < \ 1.5$) or very bottom heavy (HI > 1.5). Our results indicate that glaciers with higher HI (more bottom heavy) or lower median elevation generally have more negative SMB rates. This is not unexpected since the hypsometry influences the relative size of the accumulation and ablation areas, and therefore glacier sensitivity to winter versus summer SMB. Andreassen and others (2023) found that median elevation showed the strongest correlation with geodetic mass balance for the smaller sample of 49 glaciers on central and northern Jostedalsbreen over the shorter period 1966–2020, but lower correlation for HI. It should be noted that we omitted the detached tongue of Brenndalsbreen (ID2301, categorized as very bottom heavy) from the correlation analysis since it should be considered as an outlier following the assumption of normality underlying the Pearson correlation coefficient. We did not find strong correlations between annual SMB rates and other geometric variables (minimum elevation, maximum elevation, slope, aspect, length, area; the strongest of these is 0.23 for aspect, p ≤ 0.05).

To investigate potential topographical controls on regional patterns of SMB (Figs 5 and 6), we consider the HI of glaciers in different regions of the ice cap (North, Central, and South; Fig. 1). Most bottom-heavy glaciers (high HI) are located in the North region, consistent with the negative SMB rates found in this region. Meanwhile, nine of the 13 glaciers that can be characterized as equidimensional, bottom heavy or very bottom heavy, are located in North (e.g. Austdalsbreen and Vesledalsbreen), two are located in Central and two are located in South. However, glaciers in South generally have lower median elevations than glaciers in Central and North (98/69%, 85/50% and 63/11% of glaciers in Central, North and South, respectively, have median elevation >1500/1600 m a.s.l.). Lower median elevations in South and North compared to Central are in line with more negative summer SMB rates in these regions (Fig. 6c). However, topographic controls do not translate directly to regional patterns in annual SMB showing mostly balanced and positive SMB in South and negative SMB in North (Figs 5 and 6b).

In addition to topographical controls, the regional differences in SMB on Jostedalsbreen can be explained by spatial variability in winter precipitation on the ice cap. For example, the South receives more winter precipitation than the rest of the ice cap, which drives high winter SMB in this region and compensates for relatively large negative summer SMB (Figs 5c and 6). In addition, the North shows large temporal variability in winter SMB, with positive and negative anomalies of greater magnitude than the rest of the ice cap (55/152% of the 1960–2020 average in the 1960s/1990s). Jostedalsbreen is influenced by both frontal and orographic precipitation and precipitation amounts can show substantial local differences (e.g. Laute and Beylich, 2018). In this context it is interesting to note the magnitudes and spatial patterns of glacier-specific precipitation corrections (Fig. 3). Distributions of P _corr,j indicate that seNorge_2018 underestimates magnitudes of winter precipitation on Jostedalsbreen at different degrees, but particularly on the southeast facing and northern part of the ice cap. Due to the complex terrain around Jostedalsbreen, the ice cap's location in central western Norway and its large extent with a main ice divide stretching ~60 km from southwest to northeast, it is likely that precipitation amounts are influenced by variations in weather patterns, as well as local topographical effects. In addition to the orographic effect on precipitation, redistribution of snow by wind may play a role. However, we expect the latter to be mainly relevant on a sub-grid scale, that is, for the snow distribution across individual glacier units, but less important on a larger scale, when comparing individual glacier units or regions of the ice cap. Whether these combined effects are accurately captured in the meteorological dataset seNorge_2018 is an open question that should be subject to further investigation.

5.3 Model performance

We evaluate model performance using data that are not employed in calibration of the model: glacier-wide SMB for odd years of glaciological SMB records, point SMB from stake measurements for all available years, and geodetic mass balance for parts of the ice cap (Andreassen and others, 2023). Details of the model performance evaluation can be found in Appendix C. Modelled SMB is generally in good agreement with glacier-wide and point SMB from glaciological records (Figs 7 and 8). However, the comparison indicates that the magnitude of modelled ablation on the tongue of Nigardsbreen may be underestimated, but compensated by lower ablation at higher elevations such that modelled glacier-wide summer SMB agree well with glaciological records. This bias is supported by a slightly low value of the melt factor for ice compared to estimates from sonic ranger measurements on the tongue of Nigardsbreen in 2021 and 2022 (Appendix C). However, this bias should not be overemphasized since the estimated melt factors only reflect conditions over a narrow time frame and geographical area, while model parameter values inherently reflect average conditions. Still, underestimation of ablation on glacier tongues could result in positive biases that are exacerbated with increasing temperatures.

Modelled SMB is also in agreement with geodetic mass balance of 49 glaciers (73% of the ice-cap area) from 1966 to 2020 (Andreassen and others, 2023; Fig. 9). Considering individual glaciers, the geodetic mass-balance rate is within 1.5 times the interquartile range (IQR) of the modelled SMB rate for 34 of the 49 glaciers. Spatial variability in modelled SMB is generally in agreement with geodetic mass balance, which shows most pronounced thinning in the northeast and on low-elevation glacier tongues (Andreassen and others, 2023). However, modelled glacier-wide SMB is more negative than geodetic mass balance for glaciers in the northern part of the ice cap (e.g. ID2471, ID2474 Vesledalsbreen, ID2478 Austdalsbreen, ID2481; Fig. 9). The glacier that shows the largest discrepancy is the detached tongue of Brenndalsbreen (ID2301), where the median modelled SMB rate is −3.66 m w.e. a⁻¹, significantly more negative than the geodetic mass-balance rate of −0.54 m w.e. a⁻¹. Brenndalsbreen is fed by ice falls and avalanches from above (Engen and others, 2024) processes which are not accounted for in the SMB model. Glaciers with large positive discrepancies between modelled surface and geodetic mass balance is the upper part of Brenndalsbreen (ID2305) and Briksdalsbreen (ID2316), both located in the central western part, and Bergsetbreen (ID2318) in the central-east. Other glaciers with large positive or negative discrepancies are smaller glaciers on the margins of the ice cap (e.g. IDs 2285: west, 2258 and 2489: north, 2328 and 2333: east). It should, however, be noted that the comparison (Fig. 9) does not account for the difference in area used for calculating glacier-wide values (geodetic mass balance uses the average of the 1966 and 2019 areas). In addition, geodetic mass-balance estimates are converted from elevation to mass changes assuming a constant density, which may not reflect the spatial variability in snow, firn and ice densities across the ice cap.

It is important to mention that our model employs the seNorge_2018 DEM for the entire period 1960–2020, such that surface elevation changes are not accounted for. However, we consider the effect of surface lowering on mass balance to be negligible since surface elevation changes over Jostedalsbreen are limited over this period (Andreassen and others, 2023). The overall change in ice cap area over the modelling period is also relatively small, with a reduction of 5.2% (26.0 km²) from 1966 to 2006 and 3.4% (15.9 km²) from 2006 to 2019 (Andreassen and others, 2023). However, area changes vary between glaciers and periods, which means that for some glaciers modelled glacier-wide SMB estimates may be influenced by area changes not being properly accounted for.

5.4 Spatiotemporal variations in model parameters

Our obtained precipitation correction factors suggest that precipitation sums in seNorge_2018 are underestimated. This is in line with previous results for glaciers along the maritime-continental climate gradient in southern Norway (Sjursen and others, 2023), and corroborated by comparison of modelled accumulation to distributed snow water equivalent derived using snow depth from ground-penetrating radar measurements on the ice cap (Fig. 10, see Appendix C for details). It is not uncommon that reanalysis datasets show variable performance in capturing precipitation in mountainous regions with complex terrain (e.g. Zandler and others, 2019; Guidicelli and others, 2023), and different versions of seNorge also show discrepancies in precipitation amounts (Lussana and others, 2019). However, since we simultaneously estimate P _corr,j and T _corr,j, there are likely compensating effects of modelled ablation and accumulation on decadal SMB (Sjursen and others, 2023), such that care should be taken when interpreting the magnitude of biases. Nevertheless, the posterior of P _corr,j is unlikely to deviate strongly from its well-constrained prior (P _corr,glob). The advantage of this is that P _corr,j is informed by measurements of winter accumulation (through P _corr,glob). The disadvantage is that the spatial variability of P _corr,j will be somewhat limited. It should also be noted that parameter values are dependent on the values of fixed parameters, for example, threshold temperature for snow likely affects magnitudes of P _corr,j, as well as other data used. Observations used for parameter estimation could be afflicted with biases, for example, comparison of elevation differences from Hugonnet and others (2021) with repeat LiDAR surveys for Nigardsbreen and Austdalsbreen indicates significantly more negative geodetic mass balance using repeat LiDAR (Andreassen and others, 2023). However, there are large differences in uncertainty between the two estimates, which presents an additional argument for accounting for uncertainty in observations used to constrain models as done in this study.

Since we employ constant melt factor distributions over the ice cap in step 2 of parameter estimation, spatial patterns in T _corr,j could also reflect spatial variations in melt factors, for example, differences in solar radiation forcing between glaciers. If this was the case, we might expect more pronounced differences between predominantly north- or south-facing glaciers. Instead, variations in T _corr imply overall higher melt factors on the northwestern compared to the southeastern side of the ice cap. These differences could be explained by limitations in seNorge_2018, unresolved processes in the model and/or compensating effects of ablation and accumulation on decadal SMB. Nevertheless, constant melt factor distributions are a limitation of our model set-up, as melt factors have been shown to be transient (e.g. Gabbi and others, 2014; Ismail and others, 2023). However, we expect this temporal variability to at least partly be reflected in posterior distributions.

In relation to spatiotemporal variations in melt factors, it is interesting to compare the posterior estimate of MF _snow to values from Sjursen and others (2023), where posterior distributions of MF _snow were estimated for Nigardsbreen and Austdalsbreen individually, using seasonal glaciological SMB over the period 2000–19. Compared to this study their estimates are slightly lower for Austdalsbreen (3.53/3.51 ± 0.28 mm w.e. $^\circ$C⁻¹ d⁻¹) and somewhat higher for Nigardsbreen (4.21/4.21 ± 0.42 mm w.e. $^\circ$C⁻¹ d⁻¹), and with somewhat higher uncertainty (likely due to the combination of a larger set of observations and compensation by the model error estimated in this study). Melt factors are expected to decrease with an earlier onset of melt (Ismail and others, 2023), which may be occurring at Jostedalsbreen with the increase in spring temperatures over the past century (Hanssen-Bauer and others, 2017). However, comparison to melt factors for the recent period 2000–2019 do not indicate this. The difference between the estimates likely reflects both spatial and temporal variability, since melt factors in this study reflect variability over a longer time period and a spatial compromise between five glaciers. In addition, since posterior distributions here are jointly estimated (e.g. with T _corr), comparison of parameter values for individual glaciers should not be overemphasized.

Encompassing possible spatiotemporal variations in parameter uncertainty and model error highlights our argument for performing rigorous uncertainty estimation in SMB modelling, particularly in temperature-index models where melt processes are parameterized. This is particularly important since it is unclear if temperature-index models with constant parameter values are suitable for modelling SMB under changing climatic conditions, with studies showing contradicting evidence (e.g. Gabbi and others, 2014; Réveillet and others, 2018; Ismail and others, 2023). Energy-balance approaches have the advantage of constraining and explaining underlying physical processes. Increased availability of high-resolution climate products will alleviate their reliance on in situ meteorological data and/or downscaling of relatively coarse-resolution climate model input to the scale of the glacier. However, energy-balance models will still rely on site-specific assessment of a parameter space with significant model sensitivity (e.g. Zolles and others, 2019) that is currently more difficult to explore due to computational demands. Although the method demonstrated here may be more readily applied with temperature-index approaches due to their lower computational cost, novel methodological developments, for example, approximate Bayesian inference by using emulators to explore the relationship between parameters and observations (Cleary and others, 2021), could provide similar opportunities with more computationally expensive models.

Our parameter estimation set-up is similar to the empirical Bayesian approach of Rounce and others (2020a, 2020b, 2023), where regional prior distributions are first estimated empirically by aggregating optimized parameter values for each glacier in a region, followed by estimation of a posterior parameter distribution for each individual glacier using satellite-derived geodetic mass balances in a Bayesian model. An important novelty in this study is that we also employ a Bayesian approach to estimate the prior distribution for the glacier-specific parameter estimation using seasonal glaciological SMB (step 1), such that the prior distribution in step 2 is well-constrained and represents plausible local magnitudes of accumulation and ablation. We recognize the possibility of adopting a full Bayesian hierarchical approach (see e.g. Gelman and others, 2014) where global and glacier-specific parameters could be estimated simultaneously by assuming that glacier-specific parameters are drawn from a common population. However, this would incur significant additional computational cost and it is unclear if it would provide any additional benefits in terms of constraining model parameters and modelled SMB. Therefore, we believe that our two-step approach is sufficient in this respect and provides additional flexibility in terms of interpreting both global and glacier-specific parameter values.

With the increasing availability of satellite-borne datasets to inform SMB, we believe that SMB-modelling efforts should be adapted to take advantage of this new wealth of information to improve SMB estimates. We demonstrate one such method to leverage several observational datasets with complementary characteristics to provide robust spatiotemporal estimates of SMB over the relatively large and diverse region of Jostedalsbreen. Although seasonal glaciological SMB measurements to constrain accumulation and ablation, as used in this study, are not available in many regions of the world, it is likely that other datasets can be used to the same end, for example snowlines (e.g. Barandun and others, 2021; Geck and others, 2021) or higher resolution remote-sensing-based glacier mass balance (e.g. Belart and others, 2017; Pelto and others, 2019; Falaschi and others, 2023).

6. Conclusion

We modelled the SMB of Jostedalsbreen ice cap in western Norway over the period 1960–2020 using a temperature-index model with assimilation of both seasonal glaciological SMB observations (available for five glaciers of the ice cap) and satellite remote-sensing based decadal geodetic mass balance for the entire Jostedalsbreen. This procedure allows us to constrain winter accumulation and summer ablation, while accounting for local differences between glaciers. Overall, we found that Jostedalsbreen was nearly in balance over the past 60 years, with a small annual average mass loss of −0.07 m w.e. a⁻¹ (95% CI: −0.21, 0.08 m w.e. a⁻¹). In addition to large interannual variability in seasonal and annual SMBs, the model reveals decadal trends in SMB that can be attributed to anomalies in winter accumulation and/or summer ablation. The 1960s were characterized by mass loss, mainly attributed to low winter accumulation. In contrast, the 1990s show significant mass gains driven by high winter accumulation. Finally, substantial mass loss occurred in the 2000s, dominated by increased summer ablation due to warming air temperatures.

Our results thus suggest that SMB trends on Jostedalsbreen in the second half of the 20th century have largely been driven by variations in winter SMB due to positive and negative winter precipitation anomalies. From the start of the 21st century SMB is dominated by increased ablation due to higher temperatures, but with interannual variability influenced by variations in winter precipitation, which partly offset the effects of warming in several recent years. The SMB evolution of Jostedalsbreen stands in contrast to overall global trends that show persistently negative SMB for most glaciers.

We find that SMB varies spatially between glaciers and regions. The northern part of the ice cap and low-lying glacier tongues display the most negative rates, while the southern part shows overall positive rates. Our model reveals that spatiotemporal variations in winter accumulation and summer ablation drive SMB patterns across Jostedalsbreen. These are linked to climate variability and ongoing climate change, on the one hand, and local topographic controls, on the other hand. We expect such spatiotemporal differences in SMB-controls to have a significant influence on the future evolution of the ice cap.

Our Bayesian approach demonstrates a framework for leveraging the advantages of different information sources: the constraints on parameter values offered by in situ glaciological measurements and the unprecedented spatial coverage of satellite-derived geodetic observations that facilitate spatial adjustment of parameters to local conditions. The method allows for additional insights, such as revealing possible spatial biases in meteorological forcing data. Overall, parameter estimates indicate that winter precipitation in the seNorge_2018 meteorological dataset is underestimated over Jostedalsbreen, although possibly at different degrees both spatially and temporally.

We highlight the need for accurate mass-balance observations with sufficient temporal resolution and spatial coverage in order to constrain mass-balance models. Seasonal observations (such as provided by glaciological SMB measurements) allows the model to reproduce magnitudes of accumulation and ablation, while the spatial coverage offered by geodetic methods inform spatial variability. We therefore advocate employing complementary datasets that provide information about the spatiotemporal variability in glacier mass balance. The framework presented here illustrates an approach to utilize such datasets while simultaneously addressing the inherent uncertainties in the observations to generate robust estimates of SMB.