Spatial variability in winter mass balance on Storglaciären modelled with a terrain-based approach

Abstract

Although most processes governing the surface mass balance on mountain glaciers are well understood, the causes and extent of spatial variability in accumulation remain poorly constrained. In the present study, we couple an energy balance–snow and firn mass-balance model to terrain-based modelling routines estimating mass redistribution by snowdrift, preferential deposition and avalanching. We find this newly coupled model improves the spatial accuracy of winter balance simulations on Storglaciären, Sweden, while retaining versatility and a low computational cost. Accumulation on Storglaciären is primarily driven by direct precipitation, which is locally increased due to small-scale orographic effects. Wind-driven snow transport leads to substantial deposition in the accumulation zone and slight erosion in the ablation zone. Avalanching is the smallest contributor to winter balance, but cannot be neglected. The role of mass transporting processes in maintaining the current mass equilibrium on Storglaciären highlights the necessity to understand the links between climatic predictors and accumulation in order to accurately assess climate sensitivity.

Introduction

Surface mass balance is widely acknowledged as the essential metric of mountain glacier monitoring, providing a direct measure of an ice masses’ response to local climate (Zemp and others, 2009). The physical processes governing surface ablation are relatively well understood, with their accurate description largely driven by the advent of physically based and spatially distributed models (Hock, 2005). However, the processes governing mass accumulation remain poorly constrained (Hock and others, 2017). This gap in understanding is especially troublesome when considering small- and mid-sized mountain glaciers, as certain cirque glaciers are sustained entirely by locally enhanced accumulation rates, in areas below regional equilibrium lines (Kuhn, 1995; Hoffman and others, 2007; Florentine and others, 2018). Despite their size, the health of these small ice masses is relevant to regional hydrology and to short-term sea level rise (e.g. Huss and Hock, 2018; Parkes and Marzeion, 2018).

Spatial variability in mass accumulation on mountain glaciers has been ascribed to the variability in snow deposition in complex terrain (e.g. Jaedicke and Gauer, 2005; McGrath and others, 2018; Pramanik and others, 2019). This variability is partially driven by orographically enhanced precipitation (e.g. Mott and others, 2014; Vionnet and others, 2017) and preferential deposition (Lehning and others, 2008). In turn, post-depositional wind and gravitationally driven snow transport have a considerable impact on snow-mass distribution (Kuhn, 1995; Hoffman and others, 2007; De Beer and Sharp, 2009; Huss and Fischer, 2016; Mott and others, 2019). However, quantifying respective contributions of these mechanisms to mass balance is hindered by difficulties to gather spatially detailed snow accumulation data (Machguth and others, 2006; McGrath and others, 2018).

Reliable modelling strategies towards simulating snow transport in mountainous terrain (e.g. Christen and others, 2010; Dadic and others, 2010; Vionnet and others, 2021) have been in existence for several years and are increasingly sophisticated, but remain challenging to distribute at a spatio-temporal extent and resolution needed for their integration with distributed glacier mass balance models, which tend to have grid resolutions below 30 m and time steps below 6 hwhen simulating mountain glaciers (e.g. Reijmer and Hock, 2008; van Pelt and others, 2012). Snow drift simulations have been coupled to mass-balance models to assess wind-driven snow erosion on ice sheets (Lenaerts and others, 2012), and there is recent work towards the development of a model simulating seasonal contributions from avalanching on mountain glaciers (Turchaninova and others, 2020). To our knowledge, no effort has been made towards simulating both snow drift and avalanching on mountain glaciers at a spatial and temporal scale that allows dynamic coupling to current state of the art mass-balance modelling.

In this study we complement a distributed climatic mass-balance model (energy balance firn model – EBFM; van Pelt and others, 2012, 2021) with a routine describing wind-driven snow (re)distribution (Winstral and others, 2002) and a routine describing gravitational snow transport (Gruber, 2007). We apply this coupled snow transport-energy balance firn model (ST-EBFM) to simulate the climatic mass balance of Storglaciären, Sweden, where high-resolution winter balance measurements offer an excellent opportunity for model calibration and validation. Encouraging model accuracy shows ST-EBFM to be a numerically efficient method to account for complex and often ignored accumulation processes within glacier climatic mass-balance models. Here, our results help to assess the relative contributions of wind drift, avalanching and snowfall to the complex accumulation patterns on Storglaciären (Jansson and Pettersson, 2007).

Study area

Storglaciären is a small Scandinavian-type polythermal glacier (Pettersson and others, 2003) with a surface area of 3.1 km² and an altitudinal range between 1140 and 1700 m a.s.l. The glacier is located in the Tarfala valley in northern Sweden (67.90$^\circ$ N, 18.60$^\circ$ E), under a sub-arctic climate with temperatures consistently below freezing during the winter months (Fig. 1d; Jonsell and others, 2013). Its geometry is characterised by an upper cirque underlying the steep east face of Kebnekaise, Sweden's highest mountain, and a 2 km long valley component, also flanked by steep topography on both sides (Fig. 1b). The glacier has been subject to extensive mass-balance surveys, yielding spatially complete measurements of winter balance (Jansson and others, 1999). The spatial pattern in annual accumulation on Storglaciären is complex, with a low correlation between the yearly glacier wide winter balance (B _w) and specific winter balance (b _w) values (Jansson and Pettersson, 2007, their Fig. 2; Supplementary material S.1 in this study). Mass accumulation is especially high near the surrounding steep topography. In addition, several localised areas of enhanced and decreased accumulation occur in various places on the glacier surface. Such a pattern of winter accumulation hints at strong influences by wind and gravitationally driven processes of snow redistribution on the winter balance of Storglaciären. Possible mass contributions by avalanching have been noted in early work on the glacier (Schytt, 1965); and Jansson and Pettersson (2007) attribute the localised deviations between b _w and B _w to micro-scale wind-driven snow deposition and erosion. Both these processes are the result of an interaction between topography, wind and precipitation (Mott and others, 2018).

Fig. 1. (a) Location of study area within Scandinavia. (b) True colour ESA Sentinel-2 imagery acquired on 4th July 2021, overlain by contours spaced 50 m (Jansson and Pettersson, 2003). (c) Wind roses, obtained from ERA-5 re-analysis and from Tarfala AWS, and corrected reanalysis wind direction data used as model input. All wind rose angles are in degrees, and bins are based on frequency over the 1998–2003 period. (d) 1997–2003 average monthly temperature recorded at the Tarfala AWS and 1997–2003 average monthly cumulative precipitation estimates for Tarfala. The estimated precipitation values are obtained by taking 133% of the hourly cumulative values recorded at the Nikkaluokta AWS (Supplementary material S.2).

Data

Digital elevation model

We use a DEM of the entire Tarfala valley and Kebnekaise massif, digitised at a 15 m × 15 m grid resolution from a topographic map (Holmlund and Schytt, 1987). Within this grid, we aligned and concatenated a more recent 10 m × 10 m DEM based on GNSS surveying (Mercer, 2016), which covers the more limited spatial extent of the surface of Storglaciären itself. The resulting grid allows for detailed modelling and has sufficient spatial extent to consider reasonable ranges of snow transport.

Meteorological data

Air temperature forcing for the model consists of hourly resolution time series obtained at the Tarfala automatic weather station (AWS), located 1 km from the glacier terminus (SMHI, 2020c) at 1143 m a.s.l. (Table 1). Following Reijmer and Hock (2008), a constant lapse rate of −0.007 K m⁻¹ is applied to the measured temperature to correct for elevation.

Table 1. Overview of data used in this study

b _w, winter balance; b _s, summer balance; AWS, Automatic Weather Station; TRS, Tarfala Research Station; SMHI, Swedish Meteorological & Hydrological Institute.

Respective sources are referred to in text and included in the reference list.

The Tarfala region is subject to substantial spatial variation in precipitation, with cumulative rain and snowfall generally decreasing in the eastern direction (e.g. Holmlund and Schneider, 1997). Earlier work investigating snowfall on Storglaciären opted to use data acquired at Kråkmo, near the Norwegian coast and over 100 km away from Tarfala (Evans and others, 2008). The basis for this choice is that the snowfall close to the summit of Kebnekaise might have more in common with the western flank of the massif than with the rain shadow affected eastern side. However, individual precipitation events can be of higher intensity on alternatively the eastern and the western slope of the Scandes range, depending on the concomitant circulation pattern (Jansson and others, 2007). A comparison between temporally limited measurements made at the Tarfala AWS (SMHI, 2020c) and weather stations in Kråkmo and in Nikkaluokta (SMHI, 2020b), ~20 km eastwards of the glacier, generally yields higher correlations between the precipitation records at Tarfala and at Nikkaluokta (S.2). While precipitation at Nikkaluokta cannot also be considered fully representative of conditions on Storglaciären, the Nikkaluokta record is temporally detailed (1 h resolution) and complete. As such, it lends itself well to our purposes in this study, and we utilise the comparison with the sporadic observations at Tarfala to determine a linear correction factor of 133% between Nikkaluokta and Tarfala. As an example, Figure 1d shows a comparison of our estimated cumulative precipitation and observed values at the Tarfala AWS over the 1999 summer and fall. The Tarfala AWS hourly precipitation estimate P _tarfala is in turn distributed over the DEM while adjusting for orographic amplification and micro-scale variability with an elevation dependent linear precipitation gradient γ _P:

(1)$$P_{x, y} = P_{\rm tarfala} [ 1 + \gamma_{\rm P} ( Z_{x, y}-Z_{\rm tarfala}) ] $$

where P _x,y is the precipitation at cell x,  y on the DEM, Z _x,y is that cell's elevationand Z _tarfala is the elevation of the Tarfala weather station (1143 m a.s.l.).

The closest continuous observations of near surface wind speed and direction are acquired at the Tarfala AWS (SMHI, 2020b). Although wind speeds increase in the uphill direction due to orographically forced convergence, measurements at the AWS are thought to represent wind speeds on the glacier relatively well. Meanwhile, the prevalent wind direction from the northwest is clearly a down-valley flow, and differs from wind fields reported for the westeast-oriented Storglaciären valley (Eriksson, 2014). Lewis and others (2008) offer helpful insights towards the potential behaviour of near surface wind fields around the Kebnekaise massif: flow largely occurs longitudinally to the contours along valley bottoms, while it crosses ridges at perpendicular angles. This is in line with the observed wind directions at the AWS, but poses a problem as there are no other continuous wind direction observations in the area. An alternative data source is the ERA-5 re-analysis dataset (Hersbach and others, 2018) at the pressure level of 800 hPa, approximately equivalent to the elevation of Kebnekaise (Fig. 1b). The re-analysis data represent larger scale air flows, and is thus not influenced by small-scale topography. This means that some assumptions need to be made before it is considered representative for the surface of Storglaciären. Here we assume that wind directions within the Storglaciären valley primarily follow the westeast axis of the valley. Air flow is further likely to be directed along a similar westeast axis across the summit ridge of Kebnekaise, and along the westeast direction of the adjacent Rabots Glaciär valley. This pattern has been identified for near-surface winds during the summer season in Eriksson (2014), and we expect winter conditions to be similar. In order to simulate this, ERA-5 hourly winds are redirected to $95^\circ$ clockwise from north if the original direction was given between $10^\circ$ and $270^\circ$, and redirected to $265^\circ$ from north if the original direction was between $190^\circ$ and $330^\circ$ (Fig. 1c).

Relative humidity is measured as a percentage at the Tarfala AWS, and available at an hourly resolution (SMHI, 2020c). Cloud cover is measured at the Nikkaluokta weather station since 2004 and also available as hourly values (SMHI, 2020b). Prior to 2004, we use values measured at the SMHI weather station in Katterjåkk (SMHI, 2020a). The latter lies 60 km to the northeast of Storglaciären, and constitutes the geographically nearest available measurements. We expect the introduction of error resulting from discrepancies in cloud cover between the respective weather stations and our study site to be limited, as the radiative energy balance during the wintertime is low, and has little direct effect on mass accumulation. Finally, we correct air pressure measurements conducted at Nikkaluokta (SMHI, 2020b) for grid elevation as follows:

(2)$$P = P_{\rm SL} {\rm e}^{Z -1.2440\times10^{-4}}$$

where P _SL is sea level pressure and Z is the grid elevation in meters above sea level.

Mass-balance measurements

We calibrate and validate our model with surface mass-balance point measurements from 1995 to 2010 acquired annually in late April to early May through the Tarfala mass-balance programme (Table 1). The employed glaciological method measurement strategy for measuring winter balance consists of an extensive network of snow depth probings in a 100 m × 100 m resolution grid (Fig. 1b), combined with several snow pits and snow cores to estimate densities for the snow-depth to snow water equivalent (SWE) conversion (Holmlund and Jansson, 1999). Coverage of the probing network is near glacier wide, but certain areas were occasionally left un-sampled due to overhead hazards or crevasses. All snow-mass present on the glacier is included, regardless of suspected origin (Mercer, 2016). The accuracy of measured values is evaluated critically in Jansson (1999) at a precision of 0.1 m w.e.

Model set-up

Energy-balance firn model

The current study presents additions to a coupled energy balance–snow and firn model (EBFM), which solves the surface energy balance to calculate surface temperature and melt. EBFM includes a multi-layer subsurface component that computes snow and firn pack densities, temperature and water content. The multi-layer approach accounts for vertical liquid water motion through the snow and firn pack, as well as capillary retention and refreezing. The thermodynamic effects of this water transport are accounted for, as well as its effect upon subsurface densities. These densities are described separately for firn (Ligtenberg and others, 2011) and for seasonal snow (van Kampenhout and others, 2017), incorporating the effects of metamorphism, overburden packingand wind compaction. EBFM is commonly used to simulate glacier climatic mass balance, and extensively described in previous studies (van Pelt and others, 2012, 2019). The reader is referred to these works for details on model physics and numerical approaches used in the energy balance and melt calculations. The present model differs from EBFM in that it adds a snow transport (ST) component to the consideration of climatic mass balance, and is thus further referred to as ST-EBFM. In ST-EBFM, climatic mass balance (CMB, referring to mass change at and near to the glacier surface, Cogley and others, 2011) is calculated for a unit of time and space as:

(3)$$CMB = \int \left[P + {\rm SD} + {\rm AD} + {Q_{\rm LH}\over L_{\rm s, v}} - R\right]$$

where P is the precipitation, SD is themass from wind-driven snow transport and AD is deposits from gravitationally driven snow transport. SD can be both positive or negative while AD is always positive. The two terms distinguish Eqn (3) from earlier CMB computations. Q _LH/L _s, v is themass exchange with the atmosphere through sublimation and riming (L _s), evaporation and condensation (L _v). The formulations of these terms describe the sensible and latent heat fluxes as a function of turbulence driven by katabatic glacier wind (estimated without relying on local wind information) and near surface temperature and humidity gradients (Oerlemans and Grisogono, 2002). R denotes runoff that leaves the firn layer, meaning it is not further available as stored water or for refreezing (van Pelt and others, 2012).

Wind-driven snow transport model

Several studies have presented modelling approaches aiming at describing the influence of wind on the spatial distribution of snow-mass: notably, horizontal and vertical wind fields have been successfully modelled over glacierised terrain using a non-hydrostatic regional model (Dadic and others, 2010). Although the simulation results are compared with snow depth data, no quantification of transported volumes are made in the paper. A largely physically based model of snowdrift couples snow cover properties to wind speed and a terrain parameterisation based on curvature to simulate snow deposition over space (Liston and Sturm, 1998; Liston and others, 2007). More recent and increasingly sophisticated approaches model snowdrift coupled to snowpack models at high spatio-temporal scales by utilising downscaling techniques from localised weather models (Mott and Lehning, 2010; Vionnet and others, 2021). For the purpose of integrating wind-driven snow redistribution within a glacier mass-balance model however, the approaches mentioned above carry the significant downside of being computationally expensive and heavily reliant on high-quality input data (e.g. Marsh and others, 2020). A more parameterised approach, described in Winstral and others (2002), is preferred in this study.

Reasoning that snow transport by wind follows a continuity principle and arguing that terrain is the main driver of snow redistribution patterns, Winstral and others (2002) propose a modelling approach that solely relies on topographic analysis. The model strategy consists of quantifying the degree of likeliness that a wind speed deceleration occurs in a certain area, and that this deceleration leads to eddy formation and snow deposition. This quantification occurs through a sheltering index that centres around the maximum elevation difference to upwind topography:

(4)$$S_{x}( x_{i},\; y_{i}) = {\rm max}\left[\tan\left({Z( x_{\rm v},\; y_{\rm v}) - Z( x_{i},\; y_{i}) \over \sqrt{( x_{\rm v}-x_{i}) ^{2} + ( y_{\rm v}-y{i}) ^{2}}} \right)\right]$$

where S _x(x _i,  y _i) is a measure of suitability for snow deposition on any given cell on an xy grid, Z(x _i,  y _i) is that cell's elevation and Z(x _v,  y _v) is a vector containing the elevations of cells along the upwind direction. This sheltering index is converted to values of snow deposition through a distributed accumulation factor A _f(x _i,  y _i = f(S _x(x _i,  y _i),  P _s(x _i,  y _i)), where P _s accounts for the original elevation driven precipitation gradient. We include an overview of the model as well as its implementation within ST-EBFM in Supplementary material S.4.

Although Winstral and others's model predates the concept of preferential deposition (Lehning and others, 2008), we argue that the probabilistic nature of the terrain-based approach allows to account for wind-driven influence on both initial deposition (preferential deposition) and post-depositional redistribution (snowdrift). We note here that we include both of these processes within the category wind-driven snow redistribution in this text. Potential applications of the Winstral and others (2002) approach within glaciology have been noted in Schirmer and others (2011), and previous studies have made use of the model as a tool for discussion of the influence of topography on accumulation variability (van Pelt and others, 2014; McGrath and others, 2018).

Gravitational snow transport model

Avalanches are regularly cited as a non-negligible component of accumulation (e.g. Benn and Lehmkuhl, 2000), and the mass gain resulting from gravitational transport is known to be substantial notably for high mountain glaciers (e.g. Laha and others, 2017) and on very small cirque glaciers (e.g. Mott and others, 2019). Recent work quantified these contributions for a set of glaciers in the Caucasus mountains (Lazarev and others, 2018; Turchaninova and others, 2019). The approach considers terrain parameters to estimate potential return periods, and utilises the well-established RAMMS modellisation of Voellmy-SAM avalanche flow (Christen and others, 2010) to simulate run-out areas and seasonal deposit volumes. However, applying the strategy at higher temporal resolutions, a necessity when the aim is to couple gravitational transport to a mass-balance model, would be numerically intensive. Therefore, we use a parameterisation method proposed in Gruber (2007), where the complexity of modelling the dynamics of specific avalanches is circumvented by rather considering continuous snow redistribution. While the strategy fails to capture specific avalanching events and the flow dynamics and kinetic energy that go with it, it redistributes new snow according to maximum volume thresholds that depend on topography. A very similar method was coupled to a snow cover model in non-glacierised terrain, where the incorporation of gravitational redistribution strongly improves the model's ability to reproduce spatial snow distribution patterns (Bernhardt and Schulz, 2010).

The reader is referred to Gruber (2007) for a detailed description of the modelling approach. In addition, Supplementary material S.5 reiterates the principal components and provides a detailed description of additions made and of the model's implementation within ST-EBFM.

Numerical approach

An initialisation process computes the different sheltering indices for the set of wind directions, outlined in Supplementary material S.4, as well as the maximum deposition depth (D _max), outlined in S.5. At each time step, meteorological data inputs are extracted from their respective records. When precipitation is deemed to fall as snow because the local air temperature is higher than a set rain snow to snow transition temperature (T _S/R, van Pelt and others, 2012), wind and gravitationally driven snow transport are computed for both Storglaciären and for nearby terrain, adjusting the snow deposition following the terrain-derived parameters. All processes of snow redistribution and deposition are thus modelled simultaneously at each time step, and the spatial distribution of new snow is not subject to further transport in subsequent time steps. The surface energy balance is assessed for the cells within the glacier outline, with an iterative approach yielding surface temperature and, when at the melting point, energy available for melt. The newly computed mass input and energy balance are used in the snow model, with newly deposited snow, surface melt or erosion respectively entering or departing the uppermost layers. Finally, the net sum of added and subtracted mass yields the mass balance of the time step (Eqn (3)).

Figure 2 shows a broad outline of the model components and computation sequence. While we allow the model to iterate over every day of the considered study period, only the mass change occurring between 15th September of any given year to 15th May of the next year counts towards the winter balance to ensure simulations are representative of observed winter balance values. Every gridcell is considered at each time step. Wind and gravitationally driven mass redistribution are computed for the Tarfala wide DEM, while the energy balance, snow model and mass-balance calculations are made only for the area within the glacier outline.

Fig. 2. Flowchart of the principal model components in ST-EBFM and their sequence of operations. Purple components iterate over time steps; green components iterate over gridcells. Black components iterate over both. Components newly presented in this study are highlighted in grey.

We initialise our model domain with a spin-up run utilising climatic forcing data for the entirety of our 15 year study period, in order to achieve realistic surface and sub-surface conditions. Subsequent model runs utilise these realistic values as initial conditions. A more in depth description of the numerical routines are given in Supplementary materials S.4 and S.5.

It is important to note here a difference in approach between the terrain-based modelling of snow redistribution and the process-based modelling of the surface energy balance and subsurface components in the EBFM component of the model. The Winstral and others (2002) and Gruber (2007) terrain-based approaches to estimating snow redistribution are probabilistic: new snow-mass, incoming as precipitation, is reassigned to the location in the spatial domain where it is most likely to end up, based on parameters dictated by topography. This redistribution step occurs within the model iteration in which new snowfall enters the model. Because snow-mass is simply reassigned to its ultimate destination regardless of the pathway it takes there, the approach has the advantage of circumventing the modelling of distinct events, such as a windstorm that re-mobilises previously deposited snow or a large magnitude avalanche that fractures on a deep weak layer, that would require high detail information of snow and weather conditions, and that would be computationally expensive. We do note that we redistribute snow based on wind transport first at each model iteration, making the wind adjusted mass deposition available for gravitational transport. Avoiding the modelling of distinct events is simplistic, and models mass delivery to the glacier by wind and gravitational transport earlier and more gradually than it likely occurs in reality. The model error this could drive when dynamically coupled to the process-based surface energy balance and subsurface components of the model is considered further in the discussion.

Results

Parameter sensitivity

The parameters best suited for model calibration are identified through a brief consideration of the model response amplitude to various parameters. Table 2 column ΔSnowdrift shows the response of wind transported mass to two individual perturbations in the precipitation gradient (γ _P), the maximum sheltering distance (SD_max), and the rain to snow transition temperature (T _S/R). SD_max determines the area of the glacier surface that receives wind transported mass, and thus dominates the response. Increased precipitation allows for increased wind transport, and the resulting mass deposition increases substantially with it. Meanwhile, the response to variations in T _S/R is negligible.

Table 2. Overview of parameter sensitivity experiment results

Perturbations of the elevation-dependent precipitation gradient (γP), the rain/snow transition temperature (T _S/R), the maximum sheltering distance (SD_max), the maximum slope angle holding snow (β _lim), the mass deposition limit (D _lim) and the minimum runout angle (α _min). The units of the perturbations are identical to the parameter units. Each line represent an individually tested parameter perturbation value. Response values of wind transported snow mass (Δsnowdrift), gravitationally transported mass (ΔGrav. dep.) and total winter balance (Δb _w) shown in the three rightmost columns. All response values are glacier wide and temporal averages of 1998–2003 accumulation seasons. Notations n.t. indicate that the perturbation–response combination was not tested.

Table 2 column ΔGrav. dep. shows that the responses of gravitationally transported mass to T _S/R is negligible and to SD_max are considerably smaller than the responses to other parameters. In turn, γ _P has a strong impact on gravitationally transported mass: as with wind-driven transport, increases in precipitation are concentrated through gravitational mass transport and the mass response to γ _P can be enhanced locally. A decreased minimum runout angle (α _min), lengthening the maximum reach of deposits, generally increases mass deposition. A relatively similar response is driven by the maximum slope angle retaining any mass (β _lim). In turn, a lower β _lim leads to more gravitational transport along the slopes surrounding the glacier and more mass deposited on the flat glacier. Finally, lowering D _lim leads to shallower and more widespread deposits, yielding lower mass deposition values, while increasing the maximum deposit thickness drives a considerable increase in average deposited mass.

The previous paragraphs identify SD_max to have a determining role in the volume of wind transported mass. Besides the precipitation gradient, the D _lim parameter strongly affects the amount of gravitational transport. We note that the high relative range of perturbation applied to D _lim must partly explain its high impact on transported mass. Nevertheless, this wide range of variation in D _lim is realistic as the thickness of avalanche deposits strongly varies with terrain shape and snow conditions (e.g. McClung and Schaerer, 2006). Table 2 column Δb _w shows the response of the overall winter balance to these two parameters, along with the parameters traditionally used for accumulation calibration in studies that do not incorporate post-depositional snow transport: γ _P and T _S/R (e.g. Reijmer and Hock, 2008; van Pelt and others, 2012). γ _P clearly retains a critical impact on winter balance. The impact of T _S/R remains low, confirming that winter temperatures on Storglaciären generally fall well below the rain to snow threshold.

With the negligible response of winter balance to T _S/R, the three obvious remaining calibration parameters are:

1. (1) γ _P, which clearly has a dominant impact on direct accumulation but also on wind and gravitationally transported mass;

2. (2) SD_max, which determines the spatial and volumetric extent of wind deposition, and affects the gravitational transport component;

3. (3) D _lim, which has very localised impacts and can prove useful in obtaining spatially accurate results in the distributed model.

The above list includes a parameter for each of the studied contributors to winter balance (snowfall, snowdrift and avalanching), and its order follows the scale of the respective parameter perturbations’ impact on winter balance: from glacier wide for γ _P to very localised for D _lim. The order also reflects the sequence of snow through the model: the most simple consideration is direct deposition of precipitation on the glacier, while the longest considered itinerary of any given modelled snow particle would be precipitation, wind transport and finally gravitational transport. Following this sequence to tune the respective parameters largely avoids the need of iterative re-calibration or the consideration of the entire 3-D parameter space.

Model calibration

We evaluate model performance by comparing winter balance values measured at specific probing locations with model output values at these same locations and for 15th May of each considered year. Our calibration period considers the 1998–99 winter, and the 2000–03 winters (the 1999–2000 winter is excluded as reliable observations were not readily available). On average there are 212 probed measurements available for each year, yielding a total of 1059 compared points for the calibration period. We compute the RMSE and R ² as performance metrics. The first step in the calibration process is to tune the precipitation gradient, which affects all further components of accumulation. Figure 3a shows RMSE and R ² with different values of the precipitation gradient γ _P. The two other parameters (SD_max and D _lim) are set to zero to ensure the precipitation gradient is calibrated independently from redistribution processes. Model performance is unchanged with increasing γ _P above 35% 100 m⁻¹, and shows marginal decrease only above 200% 100 m⁻¹. The uneven response is driven by the large spatial spread of the measured winter balance: above the critical value, increases in accuracy continue to be made in areas of high winter balance, however, these begin to be cancelled out by overestimation of winter balance within areas of lower accumulation. The spatial gradient in winter balance is thus largely decoupled from elevation, meaning the γ _P parameter cannot entirely account for accumulation variability. We set γ _P at +40% 100 m⁻¹ in further calibration as this achieves the highest model performance while minimising overestimations of winter balance. The choice has implications for the processes involved in snow redistribution and how they are represented in a calibrated model, which we will examine further in our discussion of calibration uncertainty.

Fig. 3. RMSE and R ² between measured and modelled b _w on 15th May during the calibration period (winters 1998–2003 excl. 1999–2000 – a total of 1059 included comparison points). Markers indicate tested model parameters (a) plotted over values of γ _P (% 100 m⁻¹), SD_max and D _lim both at 0, (b) plotted over values of SD_max (m) with γ _P fixed at 40% 100 m⁻¹ and D _lim at 0 m w.e., (c) plotted over values of D _lim (m w.e.), with γ _P fixed at 40% 100 m⁻¹ and SD_max fixed at 750 m. Optimum values used as calibrated model parameters indicated in red.

The second calibration parameter is the maximum sheltering distance SD_max, regulating the amount of wind-driven snow transport. The RMSE and R ² of model runs over various values of SD_max are shown in Fig. 3b. The maximum considered sheltering distance yielding the best model performance is 750 m.

The final tuning parameter is the maximum deposit depth of individual avalanche events. Figure 3c shows model performance over various values of D _lim. The irregular variation in RMSE and R ² reflects the complex nature of the impact of avalanches on mass balance. Since the gravitational transport component depends so closely on topography, any change in D _lim will lead to improvements in some locations on the glacier surface but deterioration in others. Model response to variation in D _lim is much less uniform than for the previous two parameters, resulting in the non-monotonic nature of the aggregate effect shown in Fig. 3c. Nevertheless, an optimal value, allowing the least aggregate error, can be identified at D _lim = 0.05 m w.e.

The model performance metrics obtained with the optimal set of tuning parameters at the various stages of calibration are summarised in Table 3. Figure 4 shows the corresponding modelled vs measured values of winter balance. It becomes clear that the addition of wind and gravitational mass transport improves the model's ability to capture the spatial variability of winter balance: the inclusion of the SD_max and D _lim parameters, i.e. accounting for snow transport in the model, increases R ² and reduces the RMSE to values that are not achieved with just the γ _P tuning parameter (Table 3). The improvement in model performance confirms the added value of capturing mechanisms of snow redistribution in winter balance modelling. Nevertheless, the inclusion of these mechanisms also increases the spread of modelled values (Fig. 4). In many cases, the location of mass contributions does not align perfectly with the location of observed high accumulation: although systematic error is reduced, random error is frequently increased.

Fig. 4. Scatter plots of measured winter balance and modelled mass balance on 15th May with optimised tuning parameters: (a) γ _P = 40% 100 m⁻¹, SD_max = 0 m, D _lim = 0 m w.e. (b) $\gamma _{\rm P} = 40\percnt$ 100 m⁻¹, SD_max = 750 m, D _lim = 0 m w.e. (c) $\gamma _{\rm P} = 40\percnt$ 100 m⁻¹, SD_max = 750 m, D _lim = 0.05 m w.e. b _w indicates the spatially averaged modelled winter balance for each parameter combination. The measured mean b _w over 1059 probing points during the calibration period is 1.28.

Table 3. Overview of tuned model parameters and resulting model performance

ΔB _w is the modelled spatially averaged winter balance minus the observed spatially averaged winter balance.

Model validation

The model, calibrated over the 1998–2003 winters, is tested by comparing independent simulation results to observed winter balance over five accumulation seasons during the 2004–10 period, excluding the 2008–09 winter for data availability (Fig. 5a). Climatic winter conditions used for model input are relatively uniform between the two periods (S.3; Jonsell and others, 2013). The simulation yields satisfactory model performance that is quite similar to those obtained during the calibration process, with an RMSE of 0.58 m w.e. and an R ² of 0.65. Although the simulation of winter balance is not entirely accurate, the validation results hint that the achieved degree of accuracy is consistent between the two 5 year periods. Additionally, the performance improves when compared to the calibrated non-modified EBFM, which yields an RMSE of 0.72 m w.e. and an R ² of 0.54 over the same period (Fig. 5b).

Fig. 5. Scatter plot of measured winter balance and modelled CMB on 15th May for the 2005–10 winters (excluding 2009). (a) ST-EBFM. The average modelled B _w is 1.24 m w.e., while the average observed B _w is 1.43 m w.e. (b) Original EBFM model, for comparative purposes. The average modelled B _w is 0.81.

Figure 6 shows a comparison of observed and modelled winter balance, temporally averaged over the validation period. Supplementary material S.6 shows a year by year comparison over the same period. ST-EBFM reproduces the spatial variation in b _w reasonably well (Fig. 6c): the error on mass balance is within ±1 m w.e. over most of the glacier surface, and the sharp increase in b _w towards the Kebnekaise headwall is well accounted for. Nevertheless several areas of underestimation exist, notably across the glacier at 3.99 × 10⁵ UTM east, where the model misses the locally increased b _w evident in Fig. 6a. Additionally, a sudden increase in underestimation of b _w by the model is apparent in Fig. 6c, coinciding with the abrupt limit put upon wind-driven deposition by the sheltering distance SD_max. In turn, modelled contributions by gravitational transport seem too spatially concentrated, leading to localised model overestimations surrounded by underestimations.

Fig. 6. Comparison of observed and modelled winter balance, temporally averaged over the validation period (2005–10 winters, excluding 2009). (a) Observed winter balance, interpolated from probe network measurements. (b) Winter balance modelled with ST-EBFM. (c) Error on ST-EBFM b _w. The error is negative when the model underestimates b _w and positive when the model overestimates b _w.

Processes contributing to accumulation

Figure 7 shows the spatial distribution of the net contribution made by direct snowfall, wind transport and gravitational transport to the winter balance, temporally averaged over the validation period. In turn, Table 4 presents the net glacier wide contributions to the winter balance made by each component of accumulation over the validation period, and further shows these contributions averaged over the entire study period, covering the winters from 1997 to 2010, with the exclusion of 1999–2000 and 2008–09. Direct snowfall is the main mass contributor to winter balance (72.4%), but wind-driven snow transport contributes a substantial share (18.7%). The smallest contribution is made by gravitational transport, but the share is non-negligible (8.9%).

Fig. 7. ST-EBFM modelled mass contributions to specific winter balance temporally averaged over the calibration period (2005–10). (a) modelled snowfall. (b) modelled deposition of wind transported mass. (c) modelled deposition of gravitationally transported mass. Empty sections are a result of the consideration of only terrain exposed to avalanching (S.5). (d) Modelled total specific winter balance.

Table 4. Modelled components of winter accumulation

Direct snowfall includes riming, which represents ~0.1% of the total mass accumulation. Upper lines show mass contributions for the 2005–2010 validation period, and are also plotted in Fig. 7. Lower lines show the averages of mass contributions to winter balance values over the combined calibration and validation period.

Snowfall

Direct precipitation, adjusted only for elevation, provides close to three quarters of the modelled winter balance (Table 4). Figure 7a shows the spatial distribution of the modelled mass gain from snowfall. The entire area of the glacier gains mass from snowfall, with snowfall contribution values increasing linearly between a minimum of +0.4 m w.e. at the glacier toe and a maximum of +1.8 m w.e. near the headwall. This consistency yields high mass gain totals, underlining the role of snowfall as the principal component of mass gain. It also indicates that other components of the model must provide the observed spatial complexity (S.1; Jansson and Pettersson, 2007).

Wind-driven snow transport

On average, 18.7% of modelled mass gain consists of snow-mass affected by wind transport, including both post-depositional redistribution and preferential deposition. Figure 7b shows that this contribution is highly variable in space. Wind deposition of up to +2.5 m w.e. occurs on the westernmost edge of the glacier, which is sheltered from westerly winds. These values decrease rapidly with distance from the main ridge of the Kebnekaise massif, and most of the glacier surface falls outside of the large-scale sheltering index under the predominant westerly wind direction, and thus undergoes little increased deposition. Here, we note some spatial variability, with light erosion occurring in convex areas of the glacier and deposition occurring in concave areas. The area most affected by erosion (−1.04 m w.e.) occurs in the southeast corner, near the ridge between Storglaciären and Björlings Glaciär, where the sheltering index is low under both westerly and easterly wind directions (S.4).

Gravitational transport

With an average of just under 9% of the total added mass, gravitational transport is the smallest of the investigated contributors to winter balance. Figure 7c shows how this is due to the spatially restricted reach of avalanche deposits. Simultaneously, gravitational transport is the highest single gridcell mass contributor, reaching +3.1 m w.e. in proximity to the east face of Kebnekaise. In such locations, avalanching seems to be extremely influential on the specific mass balance.

Discussion

Uncertainty in terrain-based accumulation modelling

Figure 5 indicates satisfactory model performance, but also shows remaining error in the estimation of b _w. Here we address potential sources of error in our model, its implementation and its validation.

Data uncertainties

Firstly, the glaciological method employed in the Tarfala mass-balance programme measures surface mass balance (e.g. Mercer, 2018), while ST-EBFM simulates climatic mass balance (Eqn (3); van Pelt and others, 2012). Previous observations and modelling have reported substantial internal accumulation on Storglaciären (Schneider and Jansson, 2004; Reijmer and Hock, 2008) and the inclusion of internal accumulation in our model could affect our ability to compare simulations with the observed surface mass balance. However, we monitor mass contribution by internal accumulation throughout the model domain, and find values consistently below $< \!2\percnt$ of B _w, indicating that this error source is minor.

A second source of error lies in the data used for model input. Snowfall is notoriously difficult to measure at automated weather stations and thus is an unresolved source of systematic error (Rasmussen and others, 2012). In the absence of detailed data, an assumption was made on wind directions (Fig. 1d) based on summertime observations (Eriksson, 2014) and a conceptual view of wind fields around the Kebnekaise massif, but this remains a source of uncertainty. A final possible error source lies in the comparison with the surface mass-balance measurements used as control data. Some uncertainty exists even in the high-quality observations of winter balance conducted on Storglaciären (c. 0.1 m w.e.; Jansson, 1999). Furthermore, the probing network shown in Fig. 1b avoids the glacier margins that are most prone to enhanced accumulation through wind transport and especially avalanching. This could lead to underestimation between 5% and 10% of the glacier wide winter balance, and much larger error locally (Jansson, 1999). The uncertainty on the observed winter balance could at least partly explain some of the substantial model overestimations shown in Fig. 6c, as we simply interpolate the observations onto the model grid in order to facilitate visualisation. More broadly, the availability of observational data for calibration and validation of modelled winter balance could be a recurrent problem in future applications of ST-EBFM. Very few glaciers are subject to as dense a probing network as Storglaciären, and avalanche prone areas will remain under-sampled due to safety concerns. Here, additional observation types could be valuable additions to control data. Direct measurements of snow-mass displacement through avalanching (Hancock and others, 2020) could help constrain the gravitational transport component individually, and snow and ice penetrating radar surveys could test whether modelled deposition locations and rates are accurate (Machguth and others, 2006). Additionally, longer term geodetic mass-balance measurements can be used to calibrate systemic biases in point mass-balance data (Zemp and others, 2013; O'Neel and others, 2019). While long-term observations of mass balance on Storglaciären agree well with geodetic estimates of mass change (Zemp and others, 2010), the point measurements used for model testing in this work are un-calibrated.

Calibration uncertainties

During calibration of the post-depositional transport processes, we find a clear optimum value for the wind drift parameter SD_max (Fig. 3b), while model performance responds erratically to variation in the avalanche deposit thickness and reach parameter D _lim (Fig. 3c). This reflects the complexity of avalanching over the entire Storglaciären catchment, as changing conditions produce spatially varied behaviour that cannot be fully captured with a single parameter.

Meanwhile, the precipitation gradient parameter γ _P drives the largest response in winter balance (Table 2), and yet model performance remains unchanged over a wide range of values of $\gamma _{\rm P}( [ 40 - 200 \% \ 100\, {\rm m}^{-1}]$; Fig. 3a). At low values of γ _P, accumulation is estimated correctly in the ablation zone but underestimated at higher elevations. At high values of γ _P, accumulation is overestimated at the toe and correct at the upper edges of the glacier. As γ _P is the first parameter to be calibrated, the lower boundary is selected and the underestimation of b _w is largely compensated for by the addition of mass contribution from wind transport and avalanching. This approach invariably yields better model performance than when continuing the calibration process with the higher boundary of γ _P: with $\gamma _{\rm P} = 200\% \, {\rm m}^{-1}$, calibration of the wind sheltering parameter SD_max yields optimal performance at SD_max = 625 m, with R ² = 0.62, RMSE = 1.85 m w.e. and $\Delta \overline {b_{\rm w}} = + 1.6$ m w.e. This is an overestimation of the winter balance, meaning model performance with $\gamma _{\rm P} = 200\% \, {\rm m}^{-1}$ could only be worsened when including gravitational transport, as mass contributions from avalanching are always positive. As such, we suggest that with the lower γ _P of 40% 100 m⁻¹ we capture the impact of the actual precipitation gradient, without over-correcting for other accumulation processes.

We note here that our calibration process only explores a rather limited number of the possible combinations of values for the three parameters. As a result, the sequence of parameter tuning may be of influence on the optimum parameter values, and consequently could also influence our results. it would be interesting for further work employing terrain-based modelling to explore the parameter space more extensively in pursuit of better model performance.

Model uncertainties

The ablation zone on Storglaciären undergoes considerable wind erosion (Jansson and Pettersson, 2007), which might be underestimated by the current model. While our simulations do produce slight erosion in large parts of the ablation zone, the Winstral and others (2002) wind transport routine considers the effects of topography on wind speed in only one dimension, that of the considered wind direction at any given time step. This means that topographical features normal to the considered wind direction are omitted when estimating erosion and deposition. Meanwhile, in instances of down-glacier air flow, tunnelling effects in narrow valleys tend to increase near surface wind velocities (e.g. Lewis and others, 2008; Duine and others, 2016). Such an effect could increase wind velocities over the lower ablation area of Storglaciären, perhaps driving increased snow erosion and explaining the low winter balance values in the area. Additionally, while we do account for sublimation of non-drifting snow on the glacier surface, we do not account for blowing snow sublimation, which could be a significant contributor to mass loss through wind erosion (Mott and others, 2018), especially in the lower ablation area where tunnelling effects could lead to higher wind speeds and an increased volume of snow in suspension. Increased snow removal in the ablation zone would allow a calibrated model to retain a higher γ _P value, enhancing accumulation zone performance without overestimating ablation area winter balance. This suggests that ST-EBFM performance could be further improved by including a parameter regulating the extent of snow erosion, but such a parameter would have little physical basis in the absence of spatially and temporally detailed information on near surface wind fields.

Similarly, the consideration of a single wind direction at any given time step means the implemented wind-driven snow transport routine does not capture localised reversals of the main flow direction in rotors and eddies. As lee side turbulence generally results in lower wind velocities and snow deposition within the affected area (e.g. Mott and others, 2010), the overall estimation of wind-driven mass contribution is likely robust to the omission. However, it could give rise to local errors, and account for the failure to simulate eddy-like depositional patterns (Jansson and Pettersson, 2007). While simulating near surface wind fields with a non-hydrostatic model would likely improve the spatial accuracy of deposition patterns (Xue and others, 2000; Dadic and others, 2010; Vionnet and others, 2021), such a model applied at the temporal resolution of ST-EBFM would be computationally intensive.

Thirdly, fresh snow is redistributed following the wind sheltering index and through the gravitational transport model at each time step: new snow-mass is instantaneously re-assigned to its most likely final deposition location, and the simulated transport flux occurs incrementally with each new snowfall rather than through more realistic abrupt events. We suggest this only marginally affects our estimation of seasonally cumulative wind-driven deposition, as the aggregate pattern of wind transport is repeatedly found to be very consistent from year to year (S.1; Schirmer and others, 2011; Winstral and Marks, 2014). Conversely, gravitational transport, and specifically slab avalanching and cornice falls, are much more intermittent as they are the result of sequences of very specific conditions (e.g. Louchet, 2021). Return periods of large avalanches are often longer than a single year, as reflected by the high inter-annual variability found near the western headwall in Supplementary material S.1. This variability is missed by the continuous gravitational transport routine.

Future applications

The approach of terrain-based accumulation modelling also has significant advantages. Because of its computationally inexpensive nature, ST-EBFM can be implemented at very high spatial resolutions or over extensive spatial domains.

Furthermore, model complexity is largely irrelevant if meteorological data limitations are not addressed first. Due to the inherent difficulties of obtaining meteorological data during wintertime in alpine environments, the availability of distributed near surface wind and snowfall data is likely to remain a constraint for most studies investigating accumulation, and physically based models of wind transport and avalanching would suffer from the lack of accurate forcing data (e.g. Roth and others, 2018). Meanwhile, topographic data are obtained through one-time acquisition, and is often readily available from previous work or through remote-sensing products (e.g. Mercer, 2016; Morin and others, 2016). With terrain as the primary forcing, the snow transport model is less reliant on high-resolution input data, and sparse weather station data or downscaling routines can be used more confidently.

Contrary to the intensively studied Storglaciären, the mass balance of many of the world's less accessible glaciers is determined from rather sparse ground truthing, often relying on extrapolation methods that generate substantial uncertainty especially of winter balance (McGrath and others, 2015; O'Neel and others, 2019). Even when more extensive mass-balance measurements are conducted following the glaciological method, values are primarily collected along glacier centre lines (Kaser and others, 2003), and are thus likely to miss high accumulation values resulting from avalanching and much of the wind-driven deposition closer to valley walls. Here, the easily implementable ST-EBFM could be a valuable tool to complement direct observations and provide estimations of winter and net balance that account for spatial variability in accumulation. The value of snow surface energy-balance modelling in generating snow-mass datasets has been demonstrated before at regional scales (van Pelt and others, 2019), and similar applications to individual glaciers might benefit from accounting for snow transport to ensure spatially accurate winter balance estimations. Finally, studies investigating the impact of meltwater pulses on glacial hydrology and flow dynamics frequently rely on spatially distributed modelling of surface melt (e.g. Vore and others, 2019). As progress is made on observing the exact location and nature of subglacial drainage, accurately simulating the location of melt and snow-mass available for melt seems increasingly important.

Several studies have noted the persistence of small glaciers below regional equilibrium line altitudes, suggesting post-depositional mass redistribution as a mechanism through which locally enhanced accumulation could sustain glacier ice (Kuhn, 1995; Hoffman and others, 2007; De Beer and Sharp, 2009; Huss and Fischer, 2016; Mott and others, 2019). Our model results are in line with this suggestion: in order to simulate the spatial variability in winter balance, we model 27.4% of annual winter balance to originate from wind and gravitationally driven snow redistribution, implying that Storglaciären's net mass balance would be substantially different in the absence of these processes. ST-EBFM could be used as a tool to quantify mass contributions from topographically driven snow redistribution for small glaciers in un-supportive climates, and to gain further insight in their climate sensitivity.

Finally, accounting for post-depositional snow transport could improve long-term simulations of mass balance. Table 4 indicates that a considerable share of mass added during winter could be missed in traditional mass-balance models. While this deficit is likely largely accounted for by increasing snowfall during calibration, the failure to simulate the spatial distribution of added mass might lead to underestimations of summer melt rates: irregular snow cover leads to exposure of patches of lower albedo firn and bare ice earlier in the melt season, enhancing melt rates (van Pelt and others, 2014). While a more in depth evaluation of the improvements to net balance simulations would be needed, it seems likely that accounting for snow transport would benefit both reconstructions of past mass balance of small glaciers and predictions of their future mass-balance evolution under specific climate scenarios.

Conclusion

This study presented a newly coupled ST-EBFM approach and its implementation towards simulating winter balance on Storglaciären. The terrain-based modelling of wind and gravitationally driven snow redistribution in the landscape improves our model's ability to reproduce observed spatial patterns in accumulation. Our simulations suggest that over 25% of Storglaciären's average winter mass gain could stem from wind and gravitational transport, underlining these processes’ role in maintaining the current mass equilibrium of the glacier. Remaining error can be attributed to simplifications around snow transport processes, uncertainties in the calibration process and uncertainties in meteorological forcing data. Nevertheless, the encouraging model performance results suggests ST-EBFM is a viable and versatile tool to assess specific winter balance, and could have applications beyond the present study.