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, Reference Zemp, Hoelzle and Haeberli2009). 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, Reference Hock2005). However, the processes governing mass accumulation remain poorly constrained (Hock and others, Reference Hock, Hutchings and Lehning2017). 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, Reference Kuhn1995; Hoffman and others, Reference Hoffman, Fountain and Achuff2007; Florentine and others, Reference Florentine, Harper, Fagre, Moore and Peitzsch2018). 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, Reference Huss and Hock2018; Parkes and Marzeion, Reference Parkes and Marzeion2018).
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, Reference Jaedicke and Gauer2005; McGrath and others, Reference McGrath2018; Pramanik and others, Reference Pramanik, Kohler, Schuler, Van Pelt and Cohen2019). This variability is partially driven by orographically enhanced precipitation (e.g. Mott and others, Reference Mott2014; Vionnet and others, Reference Vionnet2017) and preferential deposition (Lehning and others, Reference Lehning, Löwe, Ryser and Raderschall2008). In turn, post-depositional wind and gravitationally driven snow transport have a considerable impact on snow-mass distribution (Kuhn, Reference Kuhn1995; Hoffman and others, Reference Hoffman, Fountain and Achuff2007; De Beer and Sharp, Reference De Beer and Sharp2009; Huss and Fischer, Reference Huss and Fischer2016; Mott and others, Reference Mott2019). However, quantifying respective contributions of these mechanisms to mass balance is hindered by difficulties to gather spatially detailed snow accumulation data (Machguth and others, Reference Machguth, Eisen, Paul and Hoelzle2006; McGrath and others, Reference McGrath2018).
Reliable modelling strategies towards simulating snow transport in mountainous terrain (e.g. Christen and others, Reference Christen, Kowalski and Bartelt2010; Dadic and others, Reference Dadic, Mott, Lehning and Burlando2010; Vionnet and others, Reference Vionnet2021) 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, Reference Reijmer and Hock2008; van Pelt and others, Reference van Pelt2012). Snow drift simulations have been coupled to mass-balance models to assess wind-driven snow erosion on ice sheets (Lenaerts and others, Reference Lenaerts2012), and there is recent work towards the development of a model simulating seasonal contributions from avalanching on mountain glaciers (Turchaninova and others, Reference Turchaninova2020). 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, Reference van Pelt2012, Reference van Pelt, Schuler, Pohjola and Pettersson2021) with a routine describing wind-driven snow (re)distribution (Winstral and others, Reference Winstral, Elder and Davis2002) and a routine describing gravitational snow transport (Gruber, Reference Gruber2007). 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, Reference Jansson and Pettersson2007).
Study area
Storglaciären is a small Scandinavian-type polythermal glacier (Pettersson and others, Reference Pettersson, Jansson and Holmlund2003) with a surface area of 3.1 km2 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, Reference Jonsell, Hock and Duguay2013). 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, Reference Jansson, Richardson and Jonsson1999). 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, Reference Jansson and Pettersson2007, 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, Reference Schytt1965); and Jansson and Pettersson (Reference Jansson and Pettersson2007) 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, Reference Mott, Vionnet and Grünewald2018).
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, Reference Jansson and Pettersson2003). (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, Reference Holmlund and Schytt1987). Within this grid, we aligned and concatenated a more recent 10 m × 10 m DEM based on GNSS surveying (Mercer, Reference Mercer2016), 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 (Reference Reijmer and Hock2008), a constant lapse rate of −0.007 K m−1 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, Reference Holmlund and Schneider1997). 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, Reference Evans, Essery and Lucas2008). 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, Reference Jansson, Linderholm, Pettersson, Karlin and Mörth2007). 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:
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, Reference Eriksson2014). Lewis and others (Reference Lewis, Mobbs and Lehning2008) 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, Reference Hersbach2018) 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 (Reference Eriksson2014), 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:
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, Reference Holmlund and Jansson1999). 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, Reference Mercer2016). The accuracy of measured values is evaluated critically in Jansson (Reference Jansson1999) 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, Reference Ligtenberg, Helsen and van den Broeke2011) and for seasonal snow (van Kampenhout and others, Reference van Kampenhout2017), 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, Reference van Pelt2012, Reference van Pelt2019). 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, Reference Cogley2011) is calculated for a unit of time and space as:
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, Reference Oerlemans and Grisogono2002). R denotes runoff that leaves the firn layer, meaning it is not further available as stored water or for refreezing (van Pelt and others, Reference van Pelt2012).
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, Reference Dadic, Mott, Lehning and Burlando2010). 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, Reference Liston and Sturm1998; Liston and others, Reference Liston2007). 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, Reference Mott and Lehning2010; Vionnet and others, Reference Vionnet2021). 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, Reference Marsh, Pomeroy, Spiteri and Wheater2020). A more parameterised approach, described in Winstral and others (Reference Winstral, Elder and Davis2002), 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 (Reference Winstral, Elder and Davis2002) 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:
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, Reference Lehning, Löwe, Ryser and Raderschall2008), 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 (Reference Winstral, Elder and Davis2002) approach within glaciology have been noted in Schirmer and others (Reference Schirmer, Wirz, Clifton and Lehning2011), 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, Reference van Pelt2014; McGrath and others, Reference McGrath2018).
Gravitational snow transport model
Avalanches are regularly cited as a non-negligible component of accumulation (e.g. Benn and Lehmkuhl, Reference Benn and Lehmkuhl2000), and the mass gain resulting from gravitational transport is known to be substantial notably for high mountain glaciers (e.g. Laha and others, Reference Laha2017) and on very small cirque glaciers (e.g. Mott and others, Reference Mott2019). Recent work quantified these contributions for a set of glaciers in the Caucasus mountains (Lazarev and others, Reference Lazarev, Turchaninova, Seliverstov, Komarov and Sokratov2018; Turchaninova and others, Reference Turchaninova2019). 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, Reference Christen, Kowalski and Bartelt2010) 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 (Reference Gruber2007), 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, Reference Bernhardt and Schulz2010).
The reader is referred to Gruber (Reference Gruber2007) 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, Reference van Pelt2012), 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 (Reference Winstral, Elder and Davis2002) and Gruber (Reference Gruber2007) 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 (SDmax), and the rain to snow transition temperature (T S/R). SDmax 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 (SDmax), 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 SDmax 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 SDmax 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, Reference McClung and Schaerer2006). 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, Reference Reijmer and Hock2008; van Pelt and others, Reference van Pelt2012). γ 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) γ P, which clearly has a dominant impact on direct accumulation but also on wind and gravitationally transported mass;
(2) SDmax, which determines the spatial and volumetric extent of wind deposition, and affects the gravitational transport component;
(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 2 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 2 with different values of the precipitation gradient γ P. The two other parameters (SDmax 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−1, and shows marginal decrease only above 200% 100 m−1. 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−1 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 2 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−1), SDmax and D lim both at 0, (b) plotted over values of SDmax (m) with γ P fixed at 40% 100 m−1 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−1 and SDmax fixed at 750 m. Optimum values used as calibrated model parameters indicated in red.
The second calibration parameter is the maximum sheltering distance SDmax, regulating the amount of wind-driven snow transport. The RMSE and R 2 of model runs over various values of SDmax 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 2 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 SDmax and D lim parameters, i.e. accounting for snow transport in the model, increases R 2 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−1, SDmax = 0 m, D lim = 0 m w.e. (b) $\gamma _{\rm P} = 40\percnt$
100 m−1, SDmax = 750 m, D lim = 0 m w.e. (c) $\gamma _{\rm P} = 40\percnt$
100 m−1, SDmax = 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, Reference Jonsell, Hock and Duguay2013). 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 2 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 2 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 × 105 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 SDmax. 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, Reference Jansson and Pettersson2007).
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, Reference Mercer2018), while ST-EBFM simulates climatic mass balance (Eqn (3); van Pelt and others, Reference van Pelt2012). Previous observations and modelling have reported substantial internal accumulation on Storglaciären (Schneider and Jansson, Reference Schneider and Jansson2004; Reijmer and Hock, Reference Reijmer and Hock2008) 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, Reference Rasmussen2012). In the absence of detailed data, an assumption was made on wind directions (Fig. 1d) based on summertime observations (Eriksson, Reference Eriksson2014) 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, Reference Jansson1999). 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, Reference Jansson1999). 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, Reference Hancock, Eckerstorfer, Prokop and Hendrikx2020) 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, Reference Machguth, Eisen, Paul and Hoelzle2006). Additionally, longer term geodetic mass-balance measurements can be used to calibrate systemic biases in point mass-balance data (Zemp and others, Reference Zemp2013; O'Neel and others, Reference O'Neel2019). While long-term observations of mass balance on Storglaciären agree well with geodetic estimates of mass change (Zemp and others, Reference Zemp2010), 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 SDmax (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 SDmax yields optimal performance at SDmax = 625 m, with R 2 = 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−1 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, Reference Jansson and Pettersson2007), 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 (Reference Winstral, Elder and Davis2002) 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, Reference Lewis, Mobbs and Lehning2008; Duine and others, Reference Duine, Hedde, Roubin and Durand2016). 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, Reference Mott, Vionnet and Grünewald2018), 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, Reference Mott, Schirmer, Bavay, Grünewald and Lehning2010), 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, Reference Jansson and Pettersson2007). While simulating near surface wind fields with a non-hydrostatic model would likely improve the spatial accuracy of deposition patterns (Xue and others, Reference Xue, Droegemeier and Wong2000; Dadic and others, Reference Dadic, Mott, Lehning and Burlando2010; Vionnet and others, Reference Vionnet2021), 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, Reference Schirmer, Wirz, Clifton and Lehning2011; Winstral and Marks, Reference Winstral and Marks2014). 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, Reference Louchet2021). 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, Reference Roth2018). 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, Reference Mercer2016; Morin and others, Reference Morin2016). 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, Reference McGrath2015; O'Neel and others, Reference O'Neel2019). Even when more extensive mass-balance measurements are conducted following the glaciological method, values are primarily collected along glacier centre lines (Kaser and others, Reference Kaser, Fountain, Jansson, Heucke and Knaus2003), 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, Reference van Pelt2019), 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, Reference Vore, Bartholomaus, Winberry, Walter and Amundson2019). 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, Reference Kuhn1995; Hoffman and others, Reference Hoffman, Fountain and Achuff2007; De Beer and Sharp, Reference De Beer and Sharp2009; Huss and Fischer, Reference Huss and Fischer2016; Mott and others, Reference Mott2019). 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, Reference van Pelt2014). 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.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.96
Acknowledgements
We thank Veijo Pohjola for helpful discussions and feedback on an earlier version of this work. Additionally, we are grateful to the Tarfala mass-balance programme and its volunteers for the wealth of available observations of winter balance.
Open access
