# Retrieving snow depth distribution by downscaling ERA5 Reanalysis with ICESat-2 laser altimetry

Zhihao Liu, Simon Filhol, Désirée Treichler

► To cite this version:

Zhihao Liu, Simon Filhol, Désirée Treichler. Retrieving snow depth distribution by downscaling ERA5 Reanalysis with ICESat-2 laser altimetry. Cold Regions Science and Technology, 2025, 239, pp.104580. 10.1016/j.coldregions.2025.104580. hal-05443274

HAL Id: hal-05443274

https://hal.science/hal-05443274v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Retrieving snow depth distribution by downscaling ERA5 Reanalysis with ICESat-2 laser altimetry

Zhihao Liu<sup>a, ID, *</sup>, Simon Filhol<sup>a, b</sup>, Désirée Treichler<sup>a, ID</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

Dataset link: https://doi.org/10.5281/zenodo.10048875, https://github.com/GlacioHack/xdem, https://github.com/liuh886/subgrid_snow, https://doi.org/10.5281/zenodo.2572731

**Keywords:**

- Snow depth
- Subgrid variability
- Laser altimetry
- Statistical downscaling
- ICESat-2
- XGBoost

## ABSTRACT

Estimating the variability of snow depth in remote areas poses significant challenges due to limited spatial and temporal data availability. This study uses snow depth measurements from the ICESat-2 satellite laser altimeter, which are sparse in both space and time, and incorporates them with climate reanalysis data into a downscaling-calibration scheme to produce monthly gridded snow depth maps at microscale (10 m). Snow surface elevation measurements from ICESat-2 along profiles are compared to a digital elevation model to determine snow depth at each point. To efficiently turn sparse measurements into snow depth maps, a regression model is fitted to establish a relationship between the retrieved snow depth and the corresponding ERA5 Land snow depth. This relationship, referred to as subgrid variability, is then applied to downscale the monthly ERA5 Land snow depth data. The method can provide timeseries of monthly snow depth maps for the entire ERA5 time range (since 1950). We observe that the generic output should be calibrated by a small number of localized control points from a one-time field survey to reproduce the full snow depth patterns. Results show that snow depth prediction achieved a R2 model fit value of 0.81 (post-calibration) at an intermediate scale (100 m x 500 m) using datasets from airborne laser scanning (ALS) in the Hardangervidda region of southern Norway, with still good results at microscale (R2 0.34, RMSE 1.28 m, post-calibration). Bias is greatest for extremes, with very high/low snow depths being under- and overestimated, respectively. Modeled snow depth time series at the site level have a slightly smaller RMSE than ERA5 Land data, but are still consistently biased compared to measurements from meteorological stations. Despite such localized bias and a tendency towards average snow depths the model reproduces the relative snow distribution pattern very accurately, both for peak snow (Spearman's ρ 0.77) and patchy snow meltout in late spring (Matthews correlation coefficient 0.35). The method relies on globally available data and is applicable to other snow regions above the treeline. Though requiring area-specific calibration, our approach has the potential to provide snow depth maps in areas where no such data exist and can be used to extrapolate existing snow surveys in time and over larger areas. With this, it can offer valuable input data for hydrological, ecological or permafrost modeling tasks.

## 1. Introduction

In a warming world, understanding the spatio-temporal variations of seasonal snow is increasingly vital for climate impact assessments, meltwater supply (Immerzeel et al., 2020; Livneh and Badger, 2020), permafrost modeling (Gisnås et al., 2016) and ecological responses (Callaghan et al., 2011). Seasonal snow accumulates and melts away once a year in a number of climatic zones from, for instance, the forested regions of the taiga (largest terrestrial ecosystem), the open tundra of the Arctic and many high mountain ranges throughout the world (Sturm and Liston, 2021). Being white, snow is easily observable from space. However, its thin nature (typically less than 2 m) makes

it very challenging to measure snow thickness or mass from space at a large scale and with high repeatability. The complex processes driving snow metamorphism and precipitation estimation are hurdles too to rely on models for assessing seasonal snow dynamics. Observing and modeling seasonal snow at a large scale therefore remains a challenge (Tsang et al., 2022; Mudryk et al., 2020). This is exacerbated in remote and complex terrain with limited data availability (Bormann et al., 2018).

A major persistent gap is observing or estimating snow depth and water equivalent over mountain regions. In-situ point-based snow stations are generally located in gentle terrain at lower or mid-elevation,

<footnote><sup>a</sup> Department of Geosciences, University of Oslo, Oslo, Norway</footnote>

<footnote><sup>b</sup> Météo-France - CNRS, CNRM UMR 3589, Centre d'Etude de la Neige, Grenoble, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: zhihaol@geo.uio.com (Z. Liu), simon.filhol@geo.uio.no (S. Filhol), desiree.treichler@geo.uio.no (D. Treichler).</footnote>
often under-representing rugged and higher elevation (Fassnacht et al., 2018). Employing sensors on airplanes or UAVs (unmanned aircraft vehicles) offers very high-resolution data (Bühler et al., 2016), but are costly, of limited spatial extent and potentially with legal flight restrictions. However this can be an effective approach for watershed-scale snow depth mapping (Deems et al., 2013).

Otherwise, space-borne technologies can cover systematically large to global extents at regular intervals. Two main technologies exist to retrieve snow depth: radar and laser altimetry. Radar has several issues with footprint size, penetration depth, snow wetness, and modeling effects due to scattering (Dietz et al., 2012). For instance, passive microwave radar with a footprint size of up to 25 km is insufficient to capture the heterogeneous snow depth in mountainous terrain Tsang et al. (2022), Mudryk et al. (2020). Despite recent technology development in C-band radar enabling snow depth retrieval at resolutions ranging from 500 meters to 1 kilometer under dry snow conditions (Lievens et al., 2019, 2022), finer-scale snow depth remains poorly addressed (e.g. 100 m resolution, Grünewald et al. 2010, Mott et al. 2018). Another approach for high-resolution snow depth retrieval is to use laser scanning or photogrammetry. Both methods provide snow depth by collecting surface height data pre- and post-snowfall, enabling snow depth mapping by differencing elevations. Recently, there has been a significant effort to retrieve snow depth by combining snow-free digital elevation models (DEMs) with space-borne lidar measurements of the snow surface from ICESat Treichler and Kääb (2017) and its successor ICESat-2 (Neuenschwander and Pitts, 2019). This cutting-edge satellite offers high-resolution, accurate lidar elevation profile measurements of the Earth's surface, including snow-covered terrains.

Deschamps-Berger et al. (2023) derived snow depth from the ICESat-2 ATL06 products and reported an accuracy of 0.2 m (bias) and a precision (normalized median absolute deviation; NMAD) of 0.5 m for low slopes and 1.2 m for steeper areas over the upper Tuolumne basin, California, USA. Enderlin et al. (2022) found that snow depth estimates based on ICESat-2 data had a median absolute deviation (MAD) ranging from 0.2 m for slopes < 5° to over 1 m for slopes > 20°. Besso et al. (2024) questioned the varying accuracy of the ICESat-2 ATL08 product and developed a self-defined processed elevation product, which achieved a MAD of 0.14 m to 0.20 m and root mean square error (RMSE) of 0.18 m to 0.33 m for the Tuolumne Basin and Methow Valley, USA. These studies have uncovered ICESat-2 as an emerging and cost-efficient data source for snow depth and also brought attention to challenges associated with the data and elevation differencing workflow, primarily stemming from discrepancies and spatially/temporally varying inconsistencies between reference DEMs and ICESat-2. Therefore, implementing and improving this workflow requires careful co-registration and bias correction on DEMs & ICESat-2. Additionally, the sparse nature of ICESat-2 measurements presents another significant challenge to deriving comprehensive snow-depth maps: how can we extrapolate both spatially and temporally to areas outside of the ICESat-2 measurement profiles?

Another approach to understanding snow dynamics is through snow modeling. Researchers primarily use two modeling strategies to study these dynamics, namely process-based and statistical approaches. Process-based models (Lehning et al., 2006; Liston and Elder, 2006; Kim et al., 2021) incorporated physical processes, which are driven by meteorological forcing data and yield gridded snow depth products. Mazzolini et al. (2024) combined snow depth transects from the high-resolution ICESat-2 ATL03 product with snow modeling in a data assimilation framework. They spatially propagated sparse ICESat-2 snow profile information using an abstract distance measured in a feature space defined by topographical parameters and snow melt-out climatology. For the 1 km² Izas catchment in the Spanish Pyrenees, they show that adding snow depth information in addition to the traditionally used fractional snow-covered area observations improves the model skill score by 22%. However, these models are computationally

costly and thus often struggle to cover large areas or provide fine resolution. Models are hindered by complex near-surface atmospheric processes and limited data on precipitation and wind fields (Freudiger et al., 2017), leading to a new question: how to quantify the subgrid variability of snow depth (Clark et al., 2011)? The distribution patterns of snow exhibit a notable resemblance year after year due to their dependence on topography, vegetation, and consistent synoptic weather patterns (Sturm and Wagner, 2010; Parr et al., 2020). The consistent recurrence of this pattern supports the use of computationally efficient statistical approaches. Many studies seek to establish parameterizations for subgrid variability, such as snow depletion curves, snow depth elevation gradients, snow probability distribution (Mendoza et al., 2020b), subgrid snow depth coefficient of variation (Liston, 2004; He et al., 2019; Gisnås et al., 2016) or topographic correlations (Helbig and van Herwijnen, 2017; Mazzolini et al., 2024). Learning and reproducing a snow depth map at fine scales typically involves pattern recognition. Multiple-linear regression (Grünewald et al., 2013; Dvornikov et al., 2015), binary regression trees (Revuelto et al., 2014), random forests (Revuelto et al., 2020) or a convolutional neural network (Daudt et al., 2023) have been used to predict snow distribution patterns with varied performance (R² of 0.25–0.91). However, these statistical models typically require substantial training data from terrestrial or airborne sensors. Therefore, most models can hardly be transferred to other catchments or seasons (Grünewald et al., 2013; Revuelto et al., 2020). Another category of statistical models capable of generalizing subgrid variability is commonly known as downscaling models. These models are designed to refine data from coarse, broad-scale grids to localized subgrid levels (Maraun, 2019). When applying these models retrospectively or into the future, an important assumption is made: the statistical relationships remain constant over time, a condition known as stationarity. Currently, there is a limited number of statistical downscaling models applied to snow depth (Helbig et al., 2024; Tryhorn and DeGaetano, 2013). The primary obstacles involve obtaining sufficient snow depth measurements for training and testing, meteorological forcing data in high resolution and accurately recognizing variability through informative features.

In light of these advancements and challenges, we present a method using ICESat-2 data in conjunction with high-resolution DEMs and ERA5-land climate reanalysis data, to effectively generate comprehensive snow depth maps at the hillslope scale. The objectives and workflow of this study are as follows:

1. Retrieving snow depth from ICESat-2 laser altimetry data across Norway.

2. Using this data to train a machine learning-based downscaling model that accommodates spatial and temporal variations of snow depth in mountain environments.

3. Applying local scaling calibration to downscaled snow depth for a validation area, and validating output snow depth maps at different scales with in-situ observations, gridded snow model products, and meteorological stations.

4. Discussing the challenges encountered in snow depth retrieval and snow depth downscaling.

The study area of Norway/Hardangervidda mountain plateau is chosen based on the availability of validation data. Still, the workflow is designed to be globally applicable where an accurate DEM and a proper calibration dataset are available. To retrieve snow depth measurements from ICESat-2 data, we used nationally and globally available DEMs acquired during snow-free conditions. Subsequently, we statistically downscaled ERA5 Land from its ~9 km native resolution into 10 m using ICESat-2 snow depth measurements and a machine learning algorithm. To our knowledge, this study marks the first attempt to use ICESat-2 data to downscale ERA5 Land data, and the first attempt to propagate ICESat-2 snow depths in time and space using a statistical approach.
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Map of the study area, Hardangervidda mountain plateau. A photo of the typical summer landscape is displayed at the top (a, looking south towards the Finse area and Hardangervidda. Photographer: Simon Filhol). The main panel (b) shows ICESat-2 tracks, validation data from the 2008/2009 snow survey and available weather stations. ICESat-2 ATL08 snow-on data from October 2018 to October 2022 are displayed in red. The majority of tracks was surveyed only once and the few repeated tracks appear darker due to higher density. The blue ticks at the bottom highlight the ground tracks acquired during April (peak snow). The right panel shows the location of the validation area and the total number of available snow-on data segments for mainland Norway (c).</figcaption></figure>

## 2. Study area and data setting

Norway, located on the western side of the Scandinavian Peninsula in northwestern Europe, spans latitudes from 57° 58′N to 71° 11′ N and longitudes from 4° 40′ E to 30° 58′ E. The country features a diverse topography, ranging from lowland valleys to high mountains (highest peak at 2469 m a.s.l.). In winter, westerly winds bring moisture, resulting in substantial snowfall from the coast to inland areas. This snow acts as a crucial reservoir for hydropower, emphasizing the importance of estimating snow mass in mountain environments.

Hardangervidda, our validation area, is the largest mountain plateau in northern Europe, approximately 6500 km². Its comparably flat terrain is nevertheless covered by hills and troughs exposed to high winds and heavy snowfall (Fig. 1). The plateau predominantly lies above 1000 meters above sea level (m a.s.l.), featuring a low alpine ecosystem with grass heaths, dwarf shrubs, and higher areas with bare rock or lichen marsh tundra. The eastern region is characterized by its open terrain with numerous lakes and streams. The western and southern areas, reach up to 1700 m a.s.l. and act as significant orographic barriers to the prevailing westerly wind flow. As moist air masses encounter this mountain range, they are lifted and cooled, leading to increased

precipitation on the windward slopes and a subsequent decrease of precipitation on the leeward side. Snow accumulation typically begins in mid to late September at higher elevations, peaking around late April. Mean annual precipitation ranges from 750 mm to as much as 3000 mm over relatively short distances, with approximately 50%–60% of this annual precipitation falling as snow (Ketzler et al., 2021).

### 2.1. ICESat-2 ATL08 elevation data

Launched in September 2018, ICESat-2 is equipped with the Advanced Topographic Laser Altimeter System (ATLAS), which provides photon-counting lidar measurements at a global scale (Neuenschwander and Pitts, 2019). ATLAS comprises three parallel beam pairs separated by 3.3 km on the ground. Each beam pair includes a strong and weak beam separated by 90 m. The ATLAS system emits a pulse every 0.7 m along the track, covering a circular footprint with a diameter of ~15 m. At mid-latitudes, ICESat-2 ground tracks are generally not repeated but shifted for each overpass to maximize global spatial coverage. The ICESat-2 ATL08 product (level L3 A, version 5, Neuenschwander et al., 2021) offers elevation data in fixed segment
sizes of 100 meters along the ground track (Neuenschwander et al., 2022). For each segment, there are five geolocations (subsegments) in 20 m intervals. Instead of using the mean elevation of the 100 m segment (h_tte mean, e.g., Enderlin et al., 2022), we used the subsegment height referring to best-fit terrain elevation at the midpoint location of the segment (h_tte best fit 20m 2). The sub-segment heights are determined through polynomial fitting to terrain photons with slope correction and weighting (Neuenschwander et al., 2022). Norway has a total of 3968 ATL08 data granules available for analysis from 14 October 2018, to 12 October 2022. After removing invalid data, and excluding permanent ice and inland water, our dataset consists of 13,197,376 segments, including 4,778,904 snow-free segments on land and 8,418,472 segments with snow cover over land. The coverage of snow-on segments is displayed in Fig. 1 with blue ticks representing the ground tracks from April of one of the four years (as ground tracks are not repeated).

## 2.2. Snow-off elevation data

As reference ground for snow depth retrieval, we employ the Norwegian DTM1 elevation model (DTM1, Kartverket, 2022), a 1 m lidar-based product acquired by Kartverket between 2016 and 2022. As a sensitivity test of DEM resolution, we also utilize the 10-meter resolution variant from the same data provider, DTM10 (Kartverket, 2022). To demonstrate the workflow's applicability in areas without lidar-based elevation products, we incorporate global DEMs such as Copernicus GLO-30 (European Space Agency, 2021), hereafter referred to as COP30, and FABDEM (Hawker et al., 2022, Forest And Buildings Removed Copernicus DEM, hereafter referred to as FAB) as reference ground. COP30 is a 30-meter-resolution Digital Surface Model (DSM) acquired between December 2010 and January 2015 through synthetic aperture radar interferometry (TanDEM-X mission). FAB, a variant of COP30, eliminates buildings and trees using the random forest algorithm, enhancing accuracy. FAB serves as a reference for comparison with COP30.

## 2.3. Large-scale reanalysis data

ERA5 Land hourly data (version 5) (Muñoz Sabater, 2021a) is an ECMWF (European Centre for Medium-Range Weather Forecasts) reanalysis product covering the period from 1950 to the present. It describes water and energy cycles over global land areas with over 50 variables at a spatial resolution of approximately 9 km (Muñoz-Sabater et al., 2021). This reanalysis data supplies the necessary forcing data for the downscaling model to generate sub-grid products while also accounting for input errors in the model (Günther et al., 2019; Pflug et al., 2021). ERA5 Land's snow depth data (sde) represents the instantaneous snow thickness on the ground for the elevation of each grid cell, excluding snow on vegetation canopy (Muñoz Sabater, 2021a). Additionally, the ERA5 Land monthly (Muñoz Sabater, 2021b) dataset contributes instantaneous wind fields (u10, v10) at 10 m above the land surface.

## 2.4. Validation data

Our validation methods include ALS surveys, seNorge snow model data, Sentinel-2 satellite imagery, and meteorological station data. These diverse sources offer both spatial and temporal contexts for evaluating model performance:

* The ALS survey: The survey by Melvold and Skaugen (2013) provides 2 m gridded snow depths data for two winters over Hardangervidda. The survey encompassed six flight lines apart in 10 km intervals, each extending 80 km in an east-west direction with a crossline scanning width of 500 m (Fig. 1). The data were collected between 3–21 April 2008, 21–24 April 2009, and

21 September 2008 (snow-free reference). During the autumn collection period, the ground was in nearly bare condition except for few perennial snow patches (Melvold and Skaugen, 2013). The snow depth maps were regridded (averaged) to 10 m spatial resolution for this study. Figs. 5 and 6 show snow depth data for April 2008 for parts of flight line b (indicated in Fig. 1).

* The seNorge data (www.senorge.no) employs a snow model that predicts snow depth based on interpolated precipitation and temperature station observations (seNorge2018 v23.09) (Saloranta, 2012, 2016). It offers daily snow depth maps at a 1 km × 1 km grid resolution and is available via the public archive service Thredds<fnref n="1" /> of the Norwegian Meteorological Institute (MET Norway). We aggregated the daily snow depth from seNorge into monthly average values.

* Meteorological stations: this study compared the snow depth time series with three available meteorological stations in the region (Fig. 1). The weather station Sandhaug is located 50 m north of one of the ALS flight lines at an elevation of 1250 m above sea level (a.s.l.). The other station Mogen (954 m a.s.l.) is situated directly along one of the flight line. Additionally, Haukeliseter (990 m a.s.l.) is positioned between two flight lines. Monthly mean snow depth data were retrieved from MET Norway's Frost API<fnref n="2" />. Due to harsh observing conditions in our validation area, all station observations carry a median confidence level indicated by a quality flag of 2.

* Sentinel-2 satellite imagery (L2 A; Sentinel-2, 2022) was used for visual data quality checks and validation of the presence/absence of snow in the Lake Møsvatn area (Fig. 1).

## 2.5. Calibration data

For calibration purposes, we identified 7103 points where ICESat-2 tracks overlapped with ALS strips collected in April 2008 (Fig. 1). These points, distributed across the entire study area, were selected as control points to represent the study area while the rest of the area remains unseen to the downscaling model. The outputs from the downscaling model are then calibrated against these selected control points to ensure accuracy.

# 3. Methodology

## 3.1. ICESat-2 snow depth retrieval

Snow depth ($SD_{IS2}$) was derived from ICESat-2 high-resolution elevation measurements through an elevation differencing workflow. ICESat-2 ATL08 data were categorized into snow-on ($IS_{snow}$) and snow-free segments based on attributes and flags present in the ATL08 data: the snow mask (segment_snowcover) from the National Oceanic and Atmospheric Administration (NOAA) daily snow cover product (Neuenschwander et al., 2022), as well as the (brightness_flag). Snow-free segments were used for DEM co-registration and subsequent bias correction (Fig. 2, see 3.1.1, 3.1.2). The correction involved estimating the discrepancy ($\Delta_h$) between the reference DEM ($DEM_{ref}$) and ICESat-2 snow-free measurements:

$$SD_{IS2} = IS_{snow} - (DEM_{ref} + \Delta_h)$$

Moving terrain such as water surfaces and permanent ice were excluded from the analysis using the segment_landcover reference mask from the ATL08 product (based on the Copernicus Global Land Cover dataset at 100-meter spatial resolution, Neuenschwander et al., 2022).

<footnote n="1">¹ MET Norway's Thredds API: https://thredds.met.no. Last access: 11 Sep 2023.</footnote>

<footnote n="2">² Frost API, MET Norway's archive of historical weather and climate data: https://frost.met.no. Last access: 11 Sep 2023.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Flow chart of the snow depth (SD) retrieval and downscaling-calibration scheme. First, snow depth is extracted by elevation differencing. This step includes co-registration and bias correction of the DEM (blue box). Consequently, a tree-structure-based regressor for downscaling is trained and implemented to predict local variability of snow depth in any location and at any time (green box). The resulting snow depth map time series are validated over time and in space (orange box). For full definitions of acronyms, see Section 3. Satellite graphics source: NASA.gov.</figcaption></figure>

The snow depth retrieved at each point (Easting, Northing, time) provides spatially and temporally incomplete information on snow dynamics. To overcome this, we obtained daily maximum snow depth from ERA5 Land and interpolated linearly in space for each desired point (E, N, t). The ICESat-2 snow depth is then subtracted from the interpolated ERA5 Land snow depth to create a localized dependent variable, denoted as subgrid variability (Y) (Fig. 2). This variable represents the deviation from the aggregated mean snow depth at each point, used in the subsequent downscaling model (Section 3.2). The downscaling model predicts this subgrid variability, applying it to the ERA5 Land monthly snow depth data to generate snow depth time series. Due to the inherent nature of the downscaling model, raw outputs are biased towards average values while extreme values are underrepresented (i.e. conditional bias). Thus, a calibration step (Section 3.2.3) was added to better represent the full snow depth distribution. Data processing was done in Python with custom scripts that are available in a public GitHub repository (see Data availability). The processing relies on the libraries Xarray (Hoyer and Hamman, 2017) and Pandas (The pandas development team, 2024). The retrieved snow depth for Norway, spanning from October 2018 to October 2022 is available on Zenodo (DOI: https://doi.org/10.5281/zenodo.10048875).

### 3.1.1. Co-registration

Co-registration is a crucial step to align elevation datasets, with ICESat-2 snow-free data serving as a highly precise and spatially consistent reference. We used a computationally efficient gradient descent-based co-registration algorithm (for details see Supplementary Section A). The process, facilitated by the open-source xDEM tool (Xdem contributors, 2021), was applied to each DEM tile, ensuring accurate alignment across all datasets.

### 3.1.2. Bias correction

Bias correction estimates vertical discrepancies ($\Delta_h$) between DEMs and ICESat-2 data. Given that DEMs are often patched together from multiple datasets and various sensors, captured in different seasons and at different resolutions, they must be cautiously used as a reference ground surface (Hugonnet et al., 2022). Magruder et al. (2021) used ICESat-2 elevations to correct DEMs, taking canopy and slope into account. On the other hand, Tian and Shan (2021) and Enderlin et al. (2022) found that ICESat-2 ATL08 data underestimates terrain height when compared to the reference DEMs over steep terrain, and thus proposed a slope-dependent bias correction. Our study does not assert which dataset represents ground truth most accurately but focuses on quantifying the discrepancies so that we can exclude it from snow depth
elevation differencing. Essentially, estimating the discrepancies ($\hat{\Delta}_h$) between two datasets is a regression problem that can be described with physiographic parameters, vegetation conditions, and quality metrics from ICESat-2. We employ the XGBoost<fnref n="3" /> regression model (Chen and Guestrin, 2016), a gradient-boosted decision tree (GBDT) algorithm. The regression model trained based on ICESat-2 snow-free measurements, is later used to predict elevation discrepancies ($\hat{\Delta}_h$) for all other DEM grid cells where ICESat-2 snow-free data are not available (Supplementary Section B). This results in bias-corrected snow depth measurements for all ICESat-2 snow-on data points. The bias correction successfully removed a significant negative skewness in the elevation difference histograms for snow-free data points (Figure S.4). In addition, considerable dependence on first and second-order DEM derivatives (slope and curvature) was detected and removed. After bias correction, NMAD values between ICESat-2 snow-free data points and DEM datasets are reduced from 0.66 m to 0.48 m (DTM1), and from 1.87 m to less than 0.62 m (FAB) (Figure S.4).

## 3.2. Snow depth downscaling-calibration

A second XGBoost regression model is employed to downscale ERA5 Land snow depths. XGBoost has demonstrated its effectiveness in downscaling tasks, such as total water storage anomaly from satellite gravimeter (Ali et al., 2023), precipitation (Zhu et al., 2023) or wind speed (Hu et al., 2023). We train the XGBoost regression model on (bias-corrected) ICESat-2 snow depth measurements and use a comprehensive set of topo-climatic features, including snow depth from ERA5 Land (sde_era), east, north, elevation (h_te_best_fit), slope, aspect, topographic position index (TPI, Weiss 2001; see Section 3.2.1), curvature, planform curvature (planc), profile curvature (profc), cumulative wind-aspect factor ($W_{uf}$, see Section 3.2.2) and month of the year. These predictors offer valuable information into the physical driver of snowpack dynamics. For instance, slope and curvatures are basic metrics governing snow accumulation (Filhol and Sturm, 2019). We compute these terrain attributes using xDEM (Xdem contributors, 2021) based on DTM10 with Zevenbergen & Thorne algorithm (Zevenbergen and Thorne, 1987) at 10 m resolution, consistent with the model's output resolution. The time-varying wind fields are extracted from ERA5 Land monthly data (see Section 3.2.2). We assume ERA5 land products to be stationary, allowing us to apply the downscaling model trained on snow observations from 2018 to 2022 to other periods.

The XGBoost regression model uses decision trees in parallel structure to capture nonlinear relationships between snow depth subgrid variability and topo-climatic features. During the training, the structure and splits of the trees are guided by the goal of minimizing the prediction error on given loss functions. Two types of loss functions are employed in separate model versions with slightly different purposes: (1) Square error (reg: squareerror) as a loss function to estimate the conditional mean of the target variable, this was the main method used to generate the spatially distributed maps, and (2) Quantile regression (reg: quantileloss) to give probabilistic predictions, such as Q50 (median), Q25 and Q75, which are used for point-based downscaling at weather station locations, to gain insights into the uncertainties of downscaling over time by validating predictions against in-situ observation. The adoption of quantile regression is inspired by its successful application in probabilistic forecasting (Meinshausen, 2006; Zhang et al., 2018), for which no preknowledge of the target variable distribution(e.g. a normal distribution) is required while being robust to outliers. Quantile regression is therefore an ideal choice for our downscaling task as it captures the full range of possible snow depth values (e.g. extreme values) under varying climatic conditions.

### 3.2.1. Topographic position index (TPI)

The TPI is a metric used to access slope position and classify various landforms. It quantifies the difference between the elevation of a central pixel and the average elevation of its neighboring pixels (3 × 3 pixels). A TPI value of zero or near zero indicates a flat or nearly continuous slope. Positive TPI values suggest that the central pixel is significantly higher than the surrounding areas, forming a ridge or hill. Inversely, negative TPI values indicate that the central pixel is notably lower than its neighboring areas, signifying a valley. TPI has proven effective in predicting snow distribution in alpine environments (Revuelto et al., 2014; Cristea et al., 2017). To represent landforms at different scales, we used two additional indices: tpi_9 (calculated in 9 × 9-pixel windows, equivalent to 90 m × 90 m) and tpi_27 (270 m × 270 m).

### 3.2.2. Cumulative wind-aspect factor

The wind-aspect factor ($W_f$) (Bennett et al., 2022; Dvornikov et al., 2015) serves as a proxy for snow accumulation and erosion on topographic obstacles. It assigns positive values on the leeward side and negative values on the windward side of these features. We formulated the relationship between wind and aspect by a cosine function that ranges from -1 to 1 for any prevailing direction (see Fig. 3):

$$W_f = -\cos(\text{aspect} - \text{dir}_{\text{wind}})$$

where $dir_{wind}$ is the direction of the wind origin with northerly wind (blowing from north to south) referred to as 0°. This study further divided $W_f$ into leeward $W_{f_{positive}}$ and windward factors $W_{f_{negative}}$, multiplied by wind speed to the power of three (Fig. 3 b) to capture the cumulative effect of wind redistribution for each water year period.

$$W_{uf_{positive}} = \sum W_{f_{positive}} u_{wind}^3$$

$$W_{uf_{negative}} = \sum W_{f_{negative}} u_{wind}^3$$

where $u_{wind}$ is the monthly average wind speed from ERA5 Land, linearly interpolated to 10 meters resolution. The accumulation begins in September from zero until the next August when the value reaches its maximum. The value does not accumulate when the monthly average snow depth falls below 0.1 m during the annual cycle.

### 3.2.3. Calibration

In our study, the XGBoost downscaling model is tasked with predicting snow depth under a variety of conditions. We observed that XGBoost tends to produce a conservative estimate close to the mean, likely because the feature set does not consistently explain what critical conditions lead to extreme snow depths. This conservative tendency is a natural outcome of the model's objective to minimize overall prediction error, often resulting in a distribution that skews towards average conditions and under-represents the extremes (known as scaling bias, further discussed in Section 5.2). The calibration step scales the modeled snow depth values so that their distribution matches the locally observed snow depth distribution, using quantile mapping (Cannon et al., 2015; Li et al., 2010). Thereby, while preserving the relative snow depth predictions of the downscaling model, the calibration removes the scaling bias for each specified quantile by a scaling factor calculated from the control points, ensuring a good representation of local snow depth distribution:

$$x_{adj} = F_O^{-1}(F_M(x))$$

$$\Delta = \frac{x_{adj}}{x}$$

Here, $x_{adj}$ denotes the calibrated prediction. For a given percentage, as determined by the cumulative distribution function (CDF) of the XGBoost downscaling model output $F_M(x)$ and that of the observational control points $F_O(x_{adj})$, the discrepancy between two distribution's quantiles is encapsulated by the scaling factor ($\Delta$) for a given ratio.

<footnote n="3">³ The XGBoost (version 2.0.0) library can be accessed at https://xgboost.readthedocs.io/</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3. Quantifying the relationship between wind, aspect and snow redistribution. (a) The prevailing wind, e.g. from SE, results in a negative value on the windward side (snow erosion, $W_{u/aspect}$, brown) and positive values on the lee side (snow deposition, $W_{u/aspect}$, blue). Eight cardinal directions are plotted but the function works on any wind direction. (b) Research suggests that as wind speed increases, the transport rate increases exponentially with a power of 3 (adapted from Li et al. 2018).</figcaption></figure>

### 3.3. Model interpretation and assessment

To interpret the contribution of topo-climatic features in our tree-based models, we employed the SHAP values (Lundberg et al., 2020), a metric derived from game theory facilitating the understanding of the relative contribution of the model predictors (a.k.a. features). The SHAP value corresponds to the contributions of each feature to individual predictions. A feature can either play for or against the prediction (positive or negative SHAP value), and the magnitude of its SHAP value shows the significance of the feature role. The sum of all these contributions, plus the base prediction, provides the final prediction.

For fair comparisons, we aggregated all high-resolution validation datasets to the same resolutions. Model validation was conducted across two different scales. At the micro to site scale (respectively 100 to 10 m), we captured snow distribution over typical micro-terrain features such as hills and gullies using the original model output resolution of 10 m. At the mesoscale, we aggregated data into a 100 m × 500 m grid, reflecting the 500 m width of the ALS survey swath, with 100 m intervals in the transect direction to obtain a good characterization of snow depth variability across different aspects of hills and ridges. Successful prediction would result in a (near-) perfect data match in probability distribution and ranking correlation. To quantify the downscaling performance, we employ four key statistical metrics. RMSE and R2 scores evaluate the overall accuracy and fitness of the model. R2 score is computed using the standard implementation of the Python library Scikit-learn:

$$R2 = 1 - \frac{SS_{res}}{SS_{total}}$$

where $SS_{res}$ represents the variation in the data that the fitted model does not explain, expressed by the sum of squared residuals between the model output and measured data, and $SS_{total}$ is the total variation in the data, i.e., the sum of squared residuals with regard to the mean. The R2 value typically ranges from 1 (perfect fit) to 0, but can be negative if the model is evaluated on different data than used for training (as in our case) and for nonlinear models typically used with the machine learning approach (e.g. XGBoost, as used in this study). A negative R2 value means that the model performs worse than a constant function that always predicts the mean. As R2 is sensitive to the presence of bias, Spearman's rank correlation coefficient ($\rho$) is used as a fidelity metric, with a high $\rho$ indicating good similarity in spatial distribution. The Kolmogorov-Smirnov D statistic (KSD) quantifies the degree of probability distribution matching, with KSD = 0 indicating a perfect match. To compare the patchy snow distribution during melt-out season with binary snow cover data from satellite imagery, we use the area under curve (AUC) and Matthews correlation coefficient (MCC). For these metrics, a value of 1 indicates a perfect match, whereas values of 0.5 (AUC) and 0 (MCC) correspond to random guessing.

Additionally, we used variograms to quantify the model's ability to capture the spatial heterogeneity of snow depth. The semi-variance ($\gamma$)

is a measure of spatial variability, calculated for pairs of observations as half the average squared difference between values separated by a specific lag distance ($l$) (Oliver and Webster, 2014):

$$\gamma(l) = \frac{1}{2N(l)} \sum_{i=1}^{N(l)} (z(x_i) - z(x_i + l))^2$$

where $z(x_i)$ represents the snow depth at location $x_i$. The variogram indicates the rate at which correlation decreased with distance. By fitting variograms to the sum of the spherical model and Gaussian model for short and long ranges, respectively (following the method of Rolstad et al., 2009; Hugonnet et al., 2022), we identify spatial correlation of snow distribution at different scales. The variograms are computed using the xDEM tool (Xdem contributors, 2021).

## 4. Results

### 4.1. Mesoscale snow depth variability

Fig. 4 shows the snow depth maps for April 2008. Model input data, i.e., linearly interpolated snow depth from ERA5 Land (a), only represents large-scale variability. After downscaling, the fine-scale snow depth variability aligns with the topography both at the microscale (10 m, b) and aggregated to 1 km (d). In comparison, seNorge data (c) at the same spatial resolution appear smoother with less spatial variability and overestimate snow depth in the western mountains. The differences between the two data sets are smaller for April 2009 where there was generally less snow in Southern Norway (see Supplementary Section E). A comparison of April 2008 snow depths at mesoscale for flight line b is shown as a transect in Fig. 5. The ALS snow survey, downscaled output (after calibration), seNorge, ERA5 Land data and elevation are aggregated (averaged) to a resolution of 100 m × 500 m, owing the ALS survey's 500 m transect width. All datasets follow a decreasing trend from snow depths exceeding 6 m to the West (close to the coast) to around 2 m to the East (far from the coast).

The downscaled snow depths reveal impressive details in snow depth spatial variability at the mesoscale that corresponds well to the observed snow depth magnitude and spatial variability captured by the ALS snow survey. The goodness of fit is similar for the other flight strips (Figures S.7, S.8). Across the six flight lines, the calibrated downscaling models score an R2 of 0.81 and an RMSE of 0.53–0.57 m (Table 1). We achieve nearly as good performance at mesoscale when using the global DEMs COP30 and FAB for ICESat-2 snow depth sample retrieval rather than the Norwegian high-resolution datasets DTM1 and DTM10 (Figure S.7). There are only minimal differences between output snow depth maps and we did not detect any systematic biases, no matter which DEM was used for snow depth sample retrieval. The results are similarly good for the 2009 ALS snow survey (Figure S.8), with slightly lower RMSE values of 0.50 m to 0.55 m that can be attributed to generally lower snow depths in 2009.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Spatial distribution of snow depth in the Hardangervidda area, April 2008. (a) Downscaling model input, linearly interpolated snow depth from ERA5 Land. (b) Downscaled snow depth output at 10 m resolution. (c) seNorge snow depth at 1 km resolution. (d) Downscaled snow depth aggregated to 1 km resolution. Map coordinates are in meters UTM 33N.</figcaption></figure>

In comparison, R2 values of the (coarser-resolution) seNorge model data compared to the ALS survey are lower (0.49 and 0.65) and RMSE values higher (1.12 m and 0.60 m). The seNorge model shows little spatial variability, is less sensitive to topographic features and orographic barriers, and shows about 40% overestimation on the western mountain ridges but matches ALS data in the East. Melvold and Skaugen (2013) attributed the overestimation of the seNorge model to its strong reliance on weather stations located in low-elevation areas not representative for high mountains.

Upon closer inspection, ALS ground truth data exhibit many spikes on east-facing slopes (Fig. 5), likely resulting from wind redistribution and gravity processes. Our downscaling model acknowledges the importance of the aspect and wind-driven snow accumulation factors (see also Fig. 9 c,g) but does not always reproduce corresponding snow distribution in the correct places. ERA5 land snow depth input data exhibit a north-south systematic bias, shifting from overestimation to underestimation relative to the ALS data (Fig. 9 d). The model overcorrects snow depths in the Eastern part of flight line b, resulting in snow depth underestimation compared to the ALS survey.

## 4.2. Microscale snow depth variability

Fig. 6 compares the 2008 ALS snow survey with the downscaled snow depth at a 10-meter resolution. The area shown, located in the western part of flight line 2 (indicated in Fig. 1), exhibits distinct microscale landforms, with sheltered depressions hosting thick snow patches (>8 m) and wind-exposed hilltops featuring thin snow covers (close to 0 m). This area was also shown in Melvold and Skaugen (2013, Figure 3) to address the effect of spatial resolution on snow depth representation. The terrain features align with the spatial lag that exhibits significant autocorrelation of snow depth in the variogram in Fig. 10(e), corresponding to ca. 100 m and 1.3 km.

The transect line (marked in white) in the downscaled output across this varied terrain visually captures most of the observed pattern very

<caption>**Table 1**</caption>

Statistical comparison of snow depth estimates from different methods (downscaled output with snow retrievals based on different DEMs, seNorge, and ERA5 Land) against ALS snow survey data for all six flight lines (a–f) combined, for the years 2008 and 2009, respectively. R2, KSD and Spearman's ρ are unit-less, RMSE is in m.

<table><thead><tr><th rowspan="2">Dataset</th><th colspan="4">April 2008</th><th colspan="4">April 2009</th></tr><tr><th>R2</th><th>KSD</th><th>ρ</th><th>RMSE</th><th>R2</th><th>KSD</th><th>ρ</th><th>RMSE</th></tr></thead><tbody><tr><td>DTM1</td><td>0.81</td><td>0.09</td><td>0.88</td><td>0.53</td><td>0.78</td><td>0.07</td><td>0.87</td><td>0.50</td></tr><tr><td>DTM10</td><td>0.82</td><td>0.12</td><td>0.88</td><td>0.56</td><td>0.79</td><td>0.08</td><td>0.88</td><td>0.513</td></tr><tr><td>COP30</td><td>0.80</td><td>0.10</td><td>0.88</td><td>0.57</td><td>0.75</td><td>0.07</td><td>0.84</td><td>0.554</td></tr><tr><td>FAB</td><td>0.81</td><td>0.11</td><td>0.89</td><td>0.56</td><td>0.77</td><td>0.07</td><td>0.86</td><td>0.54</td></tr><tr><td>SeNorge</td><td>0.49</td><td>0.21</td><td>0.84</td><td>1.12</td><td>0.65</td><td>0.13</td><td>0.86</td><td>0.60</td></tr><tr><td>ERA5 Land</td><td>-0.79</td><td>0.27</td><td>0.71</td><td>0.82</td><td>-0.34</td><td>0.24</td><td>0.72</td><td>0.68</td></tr></tbody></table>

well. Extreme deep/shallow snow depths tend to be biased towards average values, but less so than if scaling calibration is not applied (Figure S.6). While Spearman's ρ of 0.77 indicates a strong rank correlation, the relatively low R2 value of 0.34 and high RMSE of 1.33 m (panel d) suggests a reduced statistical agreement compared to the mesoscale analyzes. The KSD of 0.12 for the transect line and a mean deviation of 0.28 m for the difference map indicate that the scaling calibration of the downscaled snow depths, a single function applied to the entire study area and not the shown sub-region specifically, does not fully reproduce extreme values and results in a slight overestimation of average snow depths in this sub-region. The residual differences appear to be correlated with terrain features and match areas with remaining snow patches in the DTM1, which we detect by comparing the DTM1 with the ALS snow-off data (Figure S.6 e).

## 4.3. Temporal variability of snow depth

The downscaling model is also able to propagate information in time rather than only space. For this purpose, we used quantile regression
<figure><img src="image_8.png" /><figcaption>Fig. 5. Snow depth profile of flight line b across Hardangervidda as shown in Fig. 1). Data is shown for April 2008 at mesoscale, i.e. averaged to 100 m × 500 m (flight strip width) cells along the flight strip length of 80 km. Model outputs of the different DEMs are nearly identical. The map strips illustrate snow depth from ERA5 Land, seNorge, ALS snow survey and downscaled model output using the DTM1 for the marked section of the datasets, before aggregation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Microscale (10 m) snow depth comparison in the west Hardangervidda, April 2008. The validated area is located along flight line b, as shown in the study area map in Fig. 1 and mesoscale snow depth profile in Fig. 5. (a) Transect line, marked white in panels b/c, (b) ALS snow survey validation strip with significant snow depth variations, (c) downscaled snow depth, (d) differences between the ALS data and model output (background: DTM1 DEM). Common horizontal scales of landforms (95 m and 1.3 km) are symbolized by black arrows in panel a.</figcaption></figure>

for the downscaling model. The uncalibrated model output is shown
in Fig. 7, to enable analysis of the local bias and its evolution over
time. Time series data are shown for the data cell (10 m pixel size)
corresponding to the locations of three stations within the Hardan-
gervidda area: Sandhaug (a), Mogen (b), and Haukeliseter (c). The
predicted interquartile range Q25–Q75 (IQR, shown in yellow) provides
an estimate of snow depth predictions uncertainty. Q50 (blue line)
represents the median prediction. Visually, Q75 shows the best match
with the weather station data whereas input ERA5 data over- and Q50
model data underestimate measured snow depths (Table 2). (R² of

0.71, 0.60, 0.76) The downscaled model (both Q50 and Q75) generally
performs similarly to the original ERA5 Land data across all stations.
R2 values range from -0.43 (Q50, Mogen) to 0.86 (Q50, Haukeliseter),
and do not describe performance well due to significant biases at the
individual site level. While the (uncalibrated) model exhibits bias in
predicting the absolute snow depths, the high Spearman's ρ (0.91 to
0.95, ERA5: 0.92 to 0.96) across all three stations suggests that the
downscaling successfully maintains the relative ordering of snow depth
time series. The better agreement of Q75 with observations (compared
to Q50) could reflect either (1) a systematic 25% underestimation bias
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Time-series validation of snow depth for three weather stations Sandhaug (a), Mogen (b), and Haukeliseter (c). Blue Dot: measured monthly mean and minimum/maximum values (light blue), Blue line: median downscaled snow depth (Q50) and interquartile range (IQR) for 25th–75th quantiles (yellow). Note: Haukeliseter and Sandhaug data availability starts from 2015 and 2014, respectively, and some years have missing data and incorrect measurements.</figcaption></figure>

<caption>**Table 2**</caption>

Statistical comparison of monthly snow depth estimates from different methods (Downscaled output and ERA5 Land) against weather station snow depth. RMSE is in meters.

<table><thead><tr><th rowspan="2">Station</th><th colspan="4">Model — Q50</th><th colspan="4">Model — Q75</th><th colspan="4">ERA5 Land</th></tr><tr><th>R2</th><th>KSD</th><th>ρ</th><th>RMSE</th><th>R2</th><th>KSD</th><th>ρ</th><th>RMSE</th><th>R2</th><th>KSD</th><th>ρ</th><th>RMSE</th></tr></thead><tbody><tr><td>Sandhaug</td><td>0.29</td><td>0.20</td><td>0.94</td><td>0.39</td><td>0.71</td><td>0.17</td><td>0.95</td><td>0.29</td><td>0.73</td><td>0.21</td><td>0.96</td><td>0.44</td></tr><tr><td>Mogen</td><td>-0.43</td><td>0.22</td><td>0.91</td><td>0.28</td><td>0.60</td><td>0.19</td><td>0.91</td><td>0.23</td><td>0.47</td><td>0.31</td><td>0.92</td><td>0.46</td></tr><tr><td>Haukeliseter</td><td>0.86</td><td>0.10</td><td>0.94</td><td>0.26</td><td>0.76</td><td>0.23</td><td>0.94</td><td>0.38</td><td>0.37</td><td>0.36</td><td>0.94</td><td>0.79</td></tr></tbody></table>

in the model, or (2) the local station measurements corresponding to the Q75 percentile of the model's predictions for similar conditions. Bias is primarily observed for high snow depth values during the peak snow season whereas the snow-free season is captured accurately. The peak snow bias is different for the three sites but has a consistent magnitude over time at each site. There are no indications of a better fit for the training period (late 2018 to late 2022) compared to earlier years. This suggests that the bias could be reduced or removed entirely by a scaling calibration such as applied to the maps resulting from the spatial propagation presented above.

## 4.4. Validating snow occurrence

Fig. 8 provides a visual comparison between the downscaling model output and Sentinel-2 imagery for June 2020 in the lake Møsvatn area (Fig. 1). The figure shows the model's accuracy in predicting the occurrence of snow during the rapid melt period. The downscaled snow map from our model visually shows a high level of agreement with the satellite snow extent on June 24, 2020, capturing the remaining snow patches aligning with topographical features (Fig. 8c, d). Panel (e) and (f) show that east-facing slopes retain more snow compared to west-facing slopes, a pattern that our model successfully captures. Closer inspection reveals that the model retains a thin snow layer in

most areas as well as minor discrepancies in snow distribution. Some of these differences can be explained by the binary nature of Sentinel-2 snow cover and a temporal mismatch, with the Sentinel-2 image taken on a specific day and the model output representing a monthly mean snow depth for June. In particular, the model tends to systematically overestimate snow coverage for specific terrain (indicated by the white circle), such as steep slopes, and underestimates snow coverage for certain terrain features (white circle). To express the match quantitatively, we converted the model output to snow occurrence/probability using a threshold of 0.3 m for full snow cover, and the Sentinel-2 data to a binary reference by computing the Normalized Snow Difference Index (using a threshold of 0.42, as recommended by the Sentinel-2 data provider Copernicus). For the area shown in Fig. 8 (extent of panels a, b), we receive a Matthews Correlation Coefficient of 0.35 and an Area Under the Curve (AUC) value of 0.76, both indicating acceptable model performance. For a dynamic view of snow occurrence variability, please refer to the video in the supplementary materials.

## 4.5. SHAP analysis and variogram assessment

The most informative relationships between input features and predicted subgrid variability are listed in descending order. Elevation
<figure><img src="image_11.png" /><figcaption>Fig. 8. Validation of snow occurrence using Sentinel-2 imagery for the Lake Møsvatn area (annotated in Fig. 1). (a) modeled snow cover distribution for June 2020, (b) Sentinel-2 image on June 24th, 2020. Panels (c) and (d) show local discrepancies in snow distribution, (e) and (f) show correctly modeled snow retention on east-facing versus west-facing slopes. The overestimation on steep slopes and underestimation in flatter terrain are annotated with white circles.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. SHAP interpretation for the downscaling model. Panel (a) ranks features by their impact on snow depth prediction. Subsequent panels (b–g) depict SHAP dependence plots, illustrating how each feature influences model predictions for the ALS survey areas in the Hardangervidda area. A high SHAP value positively influences predicted snow depths, and a low SHAP value causes lower snow depth predictions. Each dot represents an individual 10 m grid cell, and the light grey histograms on the x-axis show the relative distributions of feature values.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Spatial correlation from semivariograms based on six ALS validation strips (2008) and downscaling model outputs at 10 m resolution. Note the limited capture of variance in the transverse direction due to the ALS strips' constrained width (500 m). The vertical bar on crosses represents one standard deviation where data are sufficient to compute the statistic.</figcaption></figure>

$(h_{te\_best\_fit})$ emerges as the most significant factor, with higher elevations positively influencing subgrid snow depth (Fig. 9a). Snow depth from ERA5 Land also has a high impact, however, both very shallow snow depth and deep localized snow depth exist (Fig. 9a, b) where $(sde\_era)$ is high, owing to the low spatial resolution of the data not capturing the variable topography. Notably, relatively deeper snow (indicated by high SHAP values in Fig. 9b, c) is estimated by the downscaling model in concave terrain (negative TPI) with an east-facing slope. In contrast, thinner snow (lower SHAP values) is associated with convex (positive TPI) and west-facing slopes. The model estimates higher snow depths for the two northernmost ALS flight strips of the study area than the southern part (Figs. 9d, 5). The positive cumulative wind aspect factor (wuf\_positive) contributes to preferential snow accumulation and shows a stronger influence than the negative wind aspect factor (wuf\_negative, Fig. 9f, g). Most wind aspect factor data points have small values indicative of low wind speeds. These show much variability in their SHAP value, indicating unexplainable variance. The SHAP method cannot distinguish the contributions from correlated features such as the apparent importance and positive influence of slope, which might partly be caused by correlated features such as elevation and curvature. The phenomenon of having less snow on very steep slopes (e.g. >50°) is not observed (Fig. 9a), rather, the model associates high slopes with a positive SHAP value that indicates greater snow depths.

The empirical semivariogram in Fig. 10 shows that the downscaling models based on all DEM versions adequately capture the spatial variability and autocorrelation of snow depths in the Hardangervidda area, but suggests stronger spatial dependence at short ranges for ALS data than downscaled snow depths (ca. 80% and 50% of the total variance at 100 m lag). The distance of ca. 100 m corresponds to the wavelength of depression features commonly observed in the area (Fig. 6). Overall, the downscaling model exhibits ca. 20% less variance than the ALS data. At a distance of ca. 1300 m, the semi-variance for all datasets approaches the sill. This distance aligns closely with the typical size of hillslopes in this region, as depicted in Fig. 6. The interpolated SeNorge product, based on precipitation gradients, temperature–elevation lapse-rates and disregarding wind processes, has a resolution too coarse to resolve the hillslope-scale snow depth variability.

# 5. Discussion

## 5.1. Bias in elevation datasets

This study shows that ICESat2 ATL08 data is not just capable to retrieve snow depths, thus confirming previous findings by other

authors (Deschamps-Berger et al., 2023; Enderlin et al., 2022; Tian and Shan, 2021; Besso et al., 2024), while also improving snow depth retrieval through co-registration and rectification of DEMs. Previous studies have in common that retrieved snow depths are biased by remaining discrepancies and errors in the DEM/ICESat-2 data such as a slope-dependent bias. The bias correction introduced in our workflow can reduce errors that depend on first and second-order derivatives (slope/curvature). Several studies pointed to ICESat-2 as the primary source of this bias, as ICESat-2 ATL08 data tends to underestimate surface height under certain conditions. Moudrý et al. (2022) attributed the error to the presence of clouds and corresponding increasing atmospheric scattering effects leading to an increased photon travel time and, consequently, underestimating terrain height. Furthermore, the 100 m segment length of ATL08 (e.g. $h_{te\_mean}$ used by Enderlin et al. 2022) is considered insufficient for accurately mapping steep and rugged terrain Besso et al. (2024). Therefore, Enderlin et al. (2022) recommended using ATL08 only in areas with relatively low slopes and sparse vegetation cover. We also observed that bias is higher in the Hardangervidda area than in other, flatter parts of Norway, such as Finnmark (not shown in this paper). In our data sample, a large number of negative snow depth measurements were present before bias correction. Directly excluding these negative snow depths (without further bias correction) would result in unbalanced training samples. A correction solely based on the slope, as proposed by Enderlin et al. (2022) and Tian and Shan (2021), did not remove bias sufficiently for the ICESat-2 measurements to be used as model input for the downscaling. We find resolution-dependent biasing effects for the tested DEMs with regard to curvature that are essential to correct for accurate snow depth retrieval (Supplementary Section B). For Norway, we found that the spatial resolution of the ATL08 data is better than COP30 (30 m) when using $h_{te\_best\_fit}$ elevation values fitted to the middle 20 m of the 100-m-ATL08 segment (Supplementary Section B).

Also the snow-off DEMs hold potential limitations for ICESat-2-derived snow depths in certain terrain, especially high mountain or Arctic areas where the snow-free season is short and some snow patches remain during summer. Typically, some snow at the highest elevations is of little concern for general-purpose DEM data distributors who need to balance requirements of numerous DEM applications. The DTM1 was acquired during June/July according to metadata, under summer conditions. From comparison with ALS snow-off data (Figure S.6), we find that the Hardangervidda area was not entirely snow-free for any of the four tested DEMs nor the Arctic DEM (Porter et al., 2022). Remaining snow on the supposedly snow-free reference results in an underestimation of snow depth in parts of these areas, with typical topographic signatures (e.g., shaded depressions). Additionally, in areas with long snow cover duration there are fewer snow-free ICESat-2 segments available. Consequently, the bias correction introduced in this study may not be able to correctly capture and remove this bias, leading to underestimation of high snow depths. When comparing our downscaled snow depths with the ALS validation data, we find that on the microscale, the residuals are highly associated with a preferential snow deposition pattern. The systematic underestimation of deep snow depths is even more prominent in the downscaled output if the calibration step is not applied (Figures S.8, S.9). The bias correction and calibration steps only indirectly correct for the biasing effect of remaining snow patches on retrieved/downscaled snow depths. They also address other bias effects (see Section 5.2). Future work should thus include a targeted correction of the DEM data for remaining snow patches. This would improve our results and potentially reduce the need for the currently calibration step.

Enderlin et al. (2022) deemed COP30 data too imprecise for estimating snow depth. Encouragingly, our results demonstrate that the COP30 and FAB DEMs improved after bias correction, and performed comparably with DTM1 and DTM10 at mesoscale in treeless areas (Figure S.7). We attribute this to the bias correction step in our workflow which
successfully removes bias for all four tested DEMs (Supplementary material). We are thus confident that the regression-based bias correction is transferable to other regions for the same (global COP30 and FAB DEMs) or different DEM products. As the training area and sample were very large and diverse in this study (entire Norway), further research on the bias-correction method could examine the applicability to different landscapes/topography or the influence of the size of the area and the ICESat-2 snow-off sample. We note that the remaining elevation error of 0.48–0.62 m (NMAD after bias correction) for DTM1 and FAB, respectively, may be acceptable for deep snow or a snowpack highly affected by wind redistribution (e.g. this study), but can exceed the thickness of shallow snow packs. The influence of uncertain individual samples is mitigated by the approach of this study, where the ICESat-2 measurements are used as a training sample for a downscaling algorithm rather than on their own. The uncertainty may be reduced for regions with less snow depth variability, lowering de-facto the absolute magnitude of discrepancies between observations and model outputs. Further research is needed to better understand the effect of the snow depth sample uncertainty on the downscaled snow depth maps, in particular the poor representation of extremes. It should also be noted that in this study, terrain parameters used for downscaling are based on the DTM10. The spatial resolution of the COP30 and FAB DEMs is coarser (30 m) than our model output resolution (10 m) and they were used for snow depth retrieval only, not for downscaling. When applying our workflow to other regions the spatial resolution of the output snow depth maps will therefore be limited by the available DEM.

This study only validates the downscaled snow depths above the tree line and further work is needed to assess the performance of the proposed workflow in forested areas where ALS validation data is available. We expect that more, or different predictors may be needed to capture information on snow-canopy processes (Mazzotti et al., 2023) influencing the accumulation and ablation of snow. High-resolution lidar data allows for the removal of vegetation, but global-scale DEMs are usually based on photogrammetry or radar interferometry, which may not exclude vegetation cover. There are efforts to create global-scale vegetation-free DEMs like the FAB DEM included in this study, and further developments are expected from machine-learning algorithms. Currently, for the high elevation accuracy required for snow depth retrieval, we recommend to be cautious about using COP30 and FAB in forested areas, as vegetation is hardly (fully) removed.

## 5.2. Snow depth downscaling

Tree-structure models perform well in capturing nonlinear relationships. The model used in this study provides reasonable predictions for snow depth. However, unless calibration is applied, our approach leaves a significant amount of variability unaccounted for, particularly for periods/areas with deep snow depths as observed in the time series analysis (Fig. 7) or in the snow depth maps (Figures S.8, S.9). The underestimation of spatial variability arises from several factors: i) the regression model’s inherent limitations, ii) a likely imbalance of the training samples, and (iii) the limitation of the predictors to capture physical processes driving the redistribution of snow.

Inherently (i), deterministic regression models primarily yield conditional means and may not adequately represent extreme events, especially when key features do not fully explain variations. This often results in underfitting as the models strive to minimize overall prediction errors. Additionally (ii), the spatial distribution of snow depth is inherently scale-dependent (Melvold and Skaugen, 2013; Mott et al., 2018), meaning a sensor’s spatial resolution is critical in determining captured variability. Our training samples, derived from the ICESat-2 ATL08 elevation product, correspond to an area of ca. 20 m ×15 m, which is better than previous study cases but still averages out for finer-scale variability. ICESat-2’s sparse spatial sampling pattern with several kilometers gap between sampled snow depth profiles, combined with a 92-day revisit period of the satellite and possible cloud cover leads

to substantial gaps in seasonal data coverage. Data points with high snow depth measurements are thus likely to be relatively few. ICESat-2’s elevation data, derived from a ground-finding algorithm that uses a probability distribution function (PDF) of reflected photons (Neuen-schwander et al., 2022), can vary in quality depending on the terrain. In rugged terrain where snow depth is typically the deepest, noise or unbalanced sampling can lead to misrepresentation of extreme values.

The shortcomings of the model and data sample could be mitigated if model predictors were able to fully capture spatial variability of snow depth and its underlying physical processes (iii), which is not the case. ERA5 Land input snow depth values may provide incorrect prior information as the dataset seems to contain interpolation effects that are possibly a downscaling artefact from the production of ERA5 Land from the coarser ERA5 data (Hersbach et al., 2020). For example, snow depth values are much higher for ERA5 Land grid cells on or adjacent to large glaciers like the Jostedalsbreen ice cap in Southern Norway, and lower for grid cells on or adjacent to large water area. Consequently, although we have excluded all glaciers and lakes, neighboring ERA5 Land grid cells still contain these biasing effects. There, the predicted subgrid variability is not solely a reflection of inherent snow conditions but is significantly influenced by this spatially-dependent bias that the model needs to detect and correct for based on spatial information, i.e coordinate features (N, E). This highlights the importance of dense, localized sampling of ICESat-2 snow depths, which may be challenging for smaller study areas and for lower latitudes, given ICESat-2’s sparse sampling pattern. By expanding the training area, as done in this study, the model benefits from more samples but might also learn a more generalized, averaged representation of snow depth.

Among the predictor features, elevation emerges as the most significant, as expected from the temperature–elevation relationship, orographic and precipitation shadowing effects (Mott et al., 2018; Parr et al., 2020) where also aspect plays an important role. The TPI at a 90-meter scale is a strong predictor (Fig. 9) and corresponds to the size of terrain features in our study area. However, these predictors represent general patterns only. They are not sufficient to represent the locally different snow depth variability. Here, the cumulative effect of multiple snow transport events over an entire winter season plays a crucial role. By refining the wind-aspect factor (Bennett et al., 2022; Dvornikov et al., 2015) into new cumulative values to describe wind redistribution (i.e. accumulation or erosion), our model is able to explain a considerable amount of the observed preferential snow deposition (Fig. 9g). Erosion is not well captured (Fig. 9f) by the downscaling model. There is potential to improve wind-related predictors, as the spatial–temporal resolution of wind fields from ERA5 Land monthly data is likely insufficient for capturing complex wind-topography interactions. Downscaling of wind fields (Fiddes et al., 2022; Toumelin et al., 2023) could enhance the model’s ability to capture wind-driven variability. Additionally, the model does not account for the impact of wind on the energy balance at the snow-atmosphere interface, involving sensible and latent heat exchanges (Mott et al., 2017), a factor crucial to the dynamics of a snowpack. The choice of predictors in this study is based on our expertise and corresponds to commonly used predictors in statistical modeling of snow distribution. This choice inherently introduces bias as the predictors are not able to explain the entire variability. Future research could benefit from incorporating more features that are related to physical snow processes or the integration of snow melt-out information from fractional snow cover data from optical satellite imagery (Gascoin et al., 2019) commonly used in snow data assimilation schemes with process-based models (Margulis et al., 2016; Mazzolini et al., 2024). Used as an exploratory tool, the model might assist in learning what processes are at play and how to best represent them statistically, ideally leading to a simplification of the workflow.

The calibration step uses quantile mapping based on the observed cumulative snow depth distribution function during 2009 peak snow conditions in the Hardangervidda area. Our selected control points
cover the entire study area (~1.5 points per square kilometer), ensuring a comprehensive correction across all regions analyzed. While greatly improving our results, especially for the reproduction of extremes values, the need of a scaling calibration is a disadvantage for an extrapolation of our method to global applications. The slight underestimation of snow depths for the presented western sub-region in the validation area (Fig. 6) suggests that the single scaling function obtained from the Hardangervidda area does not fit the sub-region perfectly, and a locally adjusted scaling might result in a better fit.

Future research is needed to determine optimal training/model area sizes or transferable scaling calibration functions for areas with different topography and climate. However, extensive, other publicly available snow depth reference datasets like the Hardangervidda ALS dataset are currently not available for Norway and are generally very rare. In the western United States, the Airborne Snow Observatory data (ASO; Painter et al., 2016) includes multi-temporal snow depth maps of several large catchments that could be suitable to examine the local dependency of our downscaling approach and further develop the method. Similarly, the lidar dataset collected in Northern Alaska by Parr et al. (2020) could be used to train for an Arctic snowpack. As an alternative to using extensive validation data for model result calibration, future research could explore the use of small representative snow depth distributions such as strategically chosen, smaller snow depth maps from drone-based surveys, data from snow courses, or snow depth measurements from meteorological stations.

### 5.3. Validation

Snow depth distribution in alpine environment is influenced by distinct processes occurring at various spatial and temporal scales. Our model demonstrates high fidelity of snow distribution patterns even at the microscale (Fig. 6c), though with lower performance metrics compared to the mesoscale. Our findings complement the insights from Mott et al. (2018), which recommend a resolution finer than 50 meters to capture the spatial variability inherent to wind-driven processes. For example, Trujillo et al. (2007), Mott et al. (2011) and Mendoza et al. (2020a) report a distinct 'scale break' at around 100 m, with a stronger spatial autocorrelation of snow depths below the scale break than beyond. Such spatial correlation can be attributed to the wavelength of wind-driven snow accumulation features (Mott et al., 2018). The ALS data and our model output identified the presence of scale breaks and dependence on terrain features (95 m accounting for 79% variance for ALS data and 108 m for 60% variance for our model, Fig. 10), demonstrating the effectiveness of our scheme in capturing a high level of spatial details and the impact of landscape patterns. Notably, we also found a scale break at 1300 m by fitting the empirical variograms at two different ranges, which corresponds to the wavelength of ridges in this area (Fig. 10). These scale breaks underscore the critical resolution and minimal survey scope required to accurately capture spatial variability over the hillslope. Such insights are pivotal for designing snow surveys or evaluating snow models/products.

RMSE or R2 are standard metrics for evaluating regression model performance, focusing primarily on overall accuracy and fit. However, they may not adequately capture the variability and structural patterns in snow depth distribution, especially at the microscale. Figs. 6 and 7 show high signal fidelity (high Spearman rank correlation) however with poor R2. To enhance the evaluation of regression outcomes for snow depth subgrid variability, it is essential to consider the scales involved and incorporate additional metrics. As an alternative to the Spearman rank correlation used in this study, structural similarity index measure (SSIM, Parr et al. 2020) might be a suitable metric to capture relative similarities.

### 5.4. Application

Current limitations of ICESat-2 data, including the month-long data availability lag and significant spatial and temporal gaps, restrict its utility (Deschamps-Berger et al., 2023). Spatio-temporal propagation of ICESat-2 snow depths by use of modeling is thus the logical next step to produce timeseries of snow depth maps and make this data useful for wider applications. Hereby, standard statistical models that solely depend on topographic features for snow spatial distribution fall short on temporal variability. Our workflow resolves many of these constraints by incorporating temporal variability from ERA5 Land. Other emerging approaches to create snow depth maps from ICESat-2 data include combining ICESat-2 snow depths with snow cover information and a process-based model within a data assimilation framework (Mazzolini et al. 2024, using higher-resolution snow depths from ATL03 data) or a combination of deep learning and probabilistic data assimilation methods (Guidicelli et al. 2023, using snow depth tracks synthesized from snow depth maps to mimic ICESat-2 snow depths). Both approaches yield promising results and also use ERA5 Land data as forcing, thus in principle globally applicable, but focusing on smaller catchments and coarser spatial resolution than presented in this work. Future research could focus on comparing or merging these different approaches.

Notably, we find that the performance of aggregated downscaled snow depth is not sensitive to the absolute precision of the reference DEMs (Supplementary Section D), indicating that our approach is not limited to areas with high-resolution DEMs but can be applied using global DEMs in treeless conditions globally. However, careful consideration is required when applying this workflow to new study areas or applications:

1. Global DEMs like the Copernicus GLO-30 or FABDEM are a patchwork of data from satellite-derived DEMs and local/national data sources and may not have equally good quality in different regions,

2. Care has to be taken to ensure including enough ICESat-2 samples in/around the target area. Areas further south will have larger coverage gaps due to ICESat-2's sampling pattern that is densest at the poles. Upcoming ICESat-2 data expected for several more years will mitigate this

3. Applying calibration with the most representative field survey

4. Accounting for regional climatic variability that some area are heavily influenced by specific weather events or patterns.

The workflow described herein is applicable to a range of applications in need to assess snow distribution patterns where sparse snow observations exist. For example, our workflow can (1) mitigate data gaps in remote areas spatially and temporally, e.g. by interpolating and extrapolating weather station observations (Fig. 6, 7), (2) help to plan and reduce cost of new snow surveys (e.g. coordinating surveys with ICESat-2 coverage), (3) improve local estimates of snow water equivalent in large alpine watershed, and (4) be combined with a classical snow modeling approach. Suppose a snow model provides not just average value snow depth, but also a probability density distribution representing the subgrid variability (Gisnås et al., 2016), such a curve can then be translated into a snow depth map using our snow downscaling scheme.

In the current implementation, we observe some limitations in extrapolating the method in time, or simulating timeseries for a single point. First, our current implementation assumes stationarity of ERA5 Land data. For instance, Kouki et al. (2023) found discontinuities in the ERA5 Land snow products around 2004 due to the introduction of new satellite products in the assimilation scheme. Other stationarity issues may arise from the method used to generate ERA5 Land with successive spin-up periods (1949, 1981, 2001; Muñoz-Sabater et al., 2021) or due to changed climate conditions in the future. An alternative to expand the temporal extrapolation and improve stationarity issues
would be to employ an energy-mass balance snow model like Crocus or Snowpack (Vionnet et al., 2012; Bartelt and Lehning, 2002) using ERA5 as forcing. Such snow reanalysis approaches were found to be better at estimating snow water equivalent (SWE) globally than methods based on passive microwave observations (Mortimer et al., 2020). Our statistical approach performs well compared to the seNorge product, but there is currently no physically-based snow reanalysis available at a comparable spatial resolution in Norway. Any subsequent application of the proposed approach will require careful assessment of uncertainties and the quality of input data, including ERA5 Land stationarity. Second, simulating snow depth at the point scale as shown in Section 4.3 shows poor agreement with station measurements (Table 2), and the application of a simple correction factor to ERA5 Land timeseries would likely yield better results. This is partly due to various local processes that may affect snow distribution not captured by the model (e.g. localized wind field in respect to terrain, radiation, etc.) or locally disturbing the weather station, but also the model capabilities and design with mostly spatially constant input features. Further studies should consider adding predictor features that consider time explicitly or in a summarized way over the course of the snow season, such as a radiation or heat index (Cristea et al., 2017) or annual snow melt-out dates from fractional snow cover data.

Given careful calibration and information on snow density, the presented cost-efficient subgrid parameterization for snow depth could be used to estimate SWE or to correct precipitation/snowfall bias in snow models (Girotto et al., 2024), to eventually serving as input for discharge modeling (Helbig and van Herwijnen, 2017). Microscale snow depth maps could benefit studies of ecosystems in snow-covered regions, from habitat availability for wildlife (Liston et al., 2016) to plant phenology under the snowpack (Niittynen and Luoto, 2018). Snow depth parametrizations/distributions at the mesoscale for large areas or detailed snow depth maps at the microscale are crucial for local permafrost studies (Gisnås et al., 2016).

## 6. Conclusion

This study introduces a workflow for snow depth retrieval from ICESat-2 ATL08 and DEM data to downscale ERA5 Land snow depth data using XGBoost tree-structure machine learning models. The two datasets, ICESat-2 and ERA5 Land, have complementary resolutions in space and time that allow for the generation of accurate monthly snow depth maps at the hillslope scale. Hereby, ERA5 Land data primarily provides the temporal variability, and ICESat-2 the spatial variability in snow depths, propagated in space using terrain features and other relevant predictors to train the downscaling model. Advanced bias correction and calibration are part of the workflow to address inherent systematic errors present in the data and correct for residual bias.

(1) There are few snow depth observations available in remote areas, and to this day, no inexpensive ways to map small-scale variability exist. Here, ICESat-2 ATL08 data as presented in our workflow of co-registration and bias-correction stands out as a valuable data source.

(2) The downscaling-calibration scheme's performance to predict peak snow for 6 × 2 ALS flight strips in the Hardangervidda is very good at mesoscale (100 × 500 m, R2 values ranging from 0.74 to 0.88). At microscale (10 m), the spatial snow depth pattern is captured very well but absolute values are less so. The model is also able to represent the spatial pattern of snow melt-out during late spring as visible from snow cover satellite data.

(3) We introduce a new cumulative wind-aspect factor in the downscaling model that estimates snow wind re-distribution from ERA5 Land monthly wind fields in a cumulative way. This factor has a high predictive strength for the spatial distribution of snow depth at micro- and mesoscale in the downscaling model.

(4) The downscaling model is sensitive to systematic bias in the elevation data, like slope- and curvature-dependent bias, which is more critical in global DEMs. Our bias correction demonstrates significant improvements in such DEMs. Therefore, similar results can be obtained when using the global DEM Copernicus GLO-30 (30 m spatial resolution) compared to the Norwegian national DEM (DTM1) at 1 m spatial resolution.

(5) Post-calibration of downscaling model results is currently necessary to compensate for the model's under-representation of extreme values. This under-representation is likely caused by inherent model behavior (tending towards average values), remaining bias and the nature of the training data. Biasing factors include persistent snow patches on supposedly snow-free reference DEMs, simplified representation of wind processes, under-sampling and underestimation of high snow depths in ICESat-2 data, and temporal non-stationarity in ERA5 Land data. Future research and more adaptable calibration methods for varied scenarios may improve results and remove the need of post-calibration.

The result of this work is a scalable and explainable downscaling model. It provides a heuristic data-driven solution to model snow depth spatio-temporal variability, especially in mountainous regions. While the validation is specific to Hardangervidda in southern Norway the workflow could be applied to other non-vegetated, snowy regions of the world given an existing local calibration dataset.

## CRediT authorship contribution statement

**Zhihao Liu:** Writing – original draft, Visualization, Software, Methodology, Investigation, Data curation, Conceptualization. **Simon Filhol:** Writing – review & editing, Methodology. **Désirée Treichler:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Funding acquisition, Conceptualization.

## Funding

This work was supported by the SNOWDEPTH project, Research Council of Norway (Grant number 325519) and the research project ClimLAND, an EEA/EU collaboration grant between Norway and Romania 2014–2021, project code RO-NO-2019-0415,1290 contract no. 30/2020. It is a contribution to the strategic research initiative LATICE (University of Oslo, projectUiO/GEO103920) and also got some support from the French Ministry of Ecological Transition (convention 261168).

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

We express our gratitude to NASA for generously providing free access to the ICESat-2 data and to ECMWF for the ERA5 Reanalysis dataset. We also extend our appreciation to the Norwegian mapping agency for granting access to the national DEM, and to Kjetil Melvold/Thomas Skaugen from the Norwegian Water Resources and Energy Directorate (NVE) for the valuable ALS dataset.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.coldregions.2025.104580.
## Data availability

The training dataset, which includes retrieved snow depth data from ICESat-2 ATLO8 for Mainland Norway spanning the period from October 2018 to October 2022, is openly accessible through Zenodo with DOI: https://doi.org/10.5281/zenodo.10048875. The Gradient Descent Coregistration method is open source at https://github.com/GlacioHack/xdem. The data processing script can be accessed at https://github.com/liuh886/subgrid_snow. The validation ALS dataset is from the Norwegian Water Resources and Energy Directorate (NVE) at DOI: https://doi.org/10.5281/zenodo.2572731. The rest of the datasets are publicly accessible (data sources available in the text).

## References

Ali, S., Khorrami, B., Jehanzaib, M., Tariq, A., Ajmal, M., Arshad, A., Shafeeqe, M., Dilawar, A., Basit, I., Zhang, L., Sadri, S., Niaz, M.A., Jamil, A., Khan, S.N., 2023. Spatial downscaling of GRACE data based on XGBoost Model for improved understanding of hydrological droughts in the Indus Basin Irrigation System (IBIS). Remote. Sens. 15 (4), 873. http://dx.doi.org/10.3390/rs15040873.

Bartelt, P., Lehning, M., 2002. A physical SNOWPACK model for the swiss avalanche warning: Part I: numerical model. Cold Reg. Sci. & Technol. 35 (3), 123-145. http://dx.doi.org/10.1016/S0165-232X(02)00074-5. URL https://www.sciencedirect.com/science/article/pii/S0165232X02000745.

Bennett, K.E., Miller, G., Busey, R., Chen, M., Lathrop, E.R., Dann, J.B., Nutt, M., Crumley, R., Dillard, S.L., Dafflon, B., Kumar, J., Bolton, W.R., Wilson, C.J., Iversen, C.M., Wullschleger, S.D., 2022. Spatial patterns of snow distribution in the sub-arctic. Cryosphere 16 (8), 3269-3293. http://dx.doi.org/10.5194/tc-16-3269-2022.

Besso, H., Shean, D., Lundquist, J.D., 2024. Mountain snow depth retrievals from customized processing of ICESat-2 satellite laser altimetry. Remote Sens. Environ. 300, 113843. http://dx.doi.org/10.1016/j.rse.2023.113843.

Bormann, K.J., Brown, R.D., Derksen, C., Painter, T.H., 2018. Estimating snow-cover trends from space. Nat. Clim. Chang. 8 (11), 924-928. http://dx.doi.org/10.1038/s41558-018-0318-3.

Bühler, Y., Adams, M.S., Bösch, R., Stoffel, A., 2016. Mapping snow depth in alpine terrain with unmanned aerial systems (UASs): Potential and limitations. Cryosphere 10 (3), 1075-1088. http://dx.doi.org/10.5194/tc-10-1075-2016.

Callaghan, T.V., Johansson, M., Brown, R.D., Groisman, P.Y., Labba, N., Radionov, V., Bradley, R.S., Blangy, S., Bulygina, O.N., Christensen, T.R., Colman, J.E., Essery, R.L.H., Forbes, B.C., Forchhammer, M.C., Golubev, V.N., Honrath, R.E., Juday, G.P., Meshcherskaya, A.V., Phoenix, G.K., Pomeroy, J., Rautio, A., Robinson, D.A., Schmidt, N.M., Serreze, M.C., Shevchenko, V.P., Shiklomanov, A.I., Shmakin, A.B., Sköld, P., Sturm, M., Woo, M.-k., Wood, E.F., 2011. Multiple effects of changes in Arctic Snow Cover. Ambio 40 (S1), 32-45. http://dx.doi.org/10.1007/s13280-011-0213-x.

Cannon, A.J., Sobie, S.R., Murdock, T.Q., 2015. Bias correction of GCM precipitation by Quantile Mapping: how well do methods preserve changes in quantiles and extremes? J. Clim. 28 (17), 6938-6959. http://dx.doi.org/10.1175/JCLI-D-14-00754.1.

Chen, T., Guestrin, C., 2016. XGBoost. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, San Francisco California USA, pp. 785-794. http://dx.doi.org/10.1145/2939672.2939785.

Clark, M.P., Hendrikx, J., Slater, A.G., Kavetski, D., Anderson, B., Cullen, N.J., Kerr, T., Örn Hreinsson, E., Woods, R.A., 2011. Representing spatial variability of snow water equivalent in hydrologic and land-surface models: A review. Water Resour. Res. 47 (7), http://dx.doi.org/10.1029/2011WR010745.

Cristea, N.C., Breckheimer, I., Raleigh, M.S., HilleRisLambers, J., Lundquist, J.D., 2017. An evaluation of terrain-based downscaling of fractional Snow Water Covered area data sets based on LiDAR-derived snow depth and orthoimagery. Water Resour. Res. 53 (8), 6802-6820. http://dx.doi.org/10.1002/2017WR020799.

Daudt, R.C., Wulf, H., Hafner, E.D., Bühler, Y., Schindler, K., Wegner, J.D., 2023. Snow depth estimation at country-scale with high spatial and temporal resolution. ISPRS J. Photogramm. Remote Sens. 197, 105-121. http://dx.doi.org/10.1016/j.isprsjprs.2023.01.017.

Deems, J.S., Painter, T.H., Finnegan, D.C., 2013. Lidar measurement of snow depth: A review. J. Glaciol. 59 (215), 467-479. http://dx.doi.org/10.3189/2013JoG12J154.

Deschamps-Berger, C., Gascoin, S., Shean, D., Besso, H., Guiot, A., López-Moreno, J.I., 2023. Evaluation of snow depth retrievals from ICESat-2 using airborne laser-scanning data. Cryosphere 17 (7), 2779-2792. http://dx.doi.org/10.5194/tc-17-2779-2023.

Dietz, A.J., Kuenzer, C., Gessner, U., Dech, S., 2012. Remote sensing of snow – a review of available methods. Int. J. Remote Sens. 33 (13), 4094-4134. http://dx.doi.org/10.1080/01431161.2011.640964.

Dvornikov, Y., Khomutov, A., Mullanurov, D., Ermokhina, K., Gubarkov, A., Leibman, M., 2015. GIS and field data-based modelling of snow water equivalent in Shrub Tundra. Fennia 193, 53-65. http://dx.doi.org/10.11143/46363.

Enderlin, E.M., Elkin, C.M., Gendreau, M., Marshall, H., O'Neel, S., McNeil, C., Florentine, C., Sass, L., 2022. Uncertainty of ICESat-2 ATLO6- and ATLO8-derived snow depths for glaciated and Vegetated Mountain Regions. Remote Sens. Environ. 283, 113307. http://dx.doi.org/10.1016/j.rse.2022.113307.

European Space Agency, 2021. Copernicus Global Digital Elevation Model. Distributed by OpenTopography [Dataset]. Distributed by OpenTopography, http://dx.doi.org/10.5069/G9028PQB.

Fassnacht, S.R., Brown, K.S.J., Blumberg, E.J., López Moreno, J.I., Covino, T.P., Kappas, M., Huang, Y., Leone, V., Kashipazha, A.H., 2018. Distribution of snow depth variability. Front. Earth Sci. 12 (4), 683-692. http://dx.doi.org/10.1007/s11707-018-0714-z.

Fiddes, J., Aalstad, K., Lehning, M., 2022. TopoCLIM: Rapid topography-based downscaling of regional climate model output in complex terrain v1.1. Geosci. Model. Dev. 15 (4), 1753-1768. http://dx.doi.org/10.5194/gmd-15-1753-2022.

Filhol, S., Sturm, M., 2019. The smoothing of landscapes during snowfall with no wind. J. Glaciol. 65 (250), 173-187. http://dx.doi.org/10.1017/jog.2018.104.

Freudiger, D., Kohn, I., Seibert, J., Stahl, K., Weiler, M., 2017. Snow redistribution for the hydrological modeling of Alpine Catchments. WIREs Water 4 (5), e1232. http://dx.doi.org/10.1002/wat2.1232.

Gascoin, S., Grizonnet, M., Bouchet, M., Salgues, G., Hagolle, O., 2019. Theia snow collection: High-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data. Earth Syst. Sci. Data 11 (2), 493-514. http://dx.doi.org/10.5194/essd-11-493-2019.

Girootto, M., Formetta, G., Azimi, S., Bachand, C., Cowherd, M., De Lannoy, G., Lievens, H., Modanesi, S., Raleigh, M.S., Rigon, R., Massari, C., 2024. Identifying snowfall elevation patterns by assimilating satellite-based snow depth retrievals. Sci. Total Environ. 906, 167312. http://dx.doi.org/10.1016/j.scitotenv.2023.167312.

Gisnás, K., Westermann, S., Schuler, T.V., Melvold, K., Etzelmüller, B., 2016. Small-scale variation of snow in a regional permafrost model. Cryosphere 10 (3), 1201-1215. http://dx.doi.org/10.5194/tc-10-1201-2016.

Grünewald, T., Schirmer, M., Mott, R., Lehning, M., 2010. Spatial and temporal variability of snow depth and ablation rates in a small mountain catchment. Cryosphere 4 (2), 215-225. http://dx.doi.org/10.5194/tc-4-215-2010.

Grünewald, T., Stötter, J., Pomeroy, J.W., Dadic, R., Moreno Baños, I., Marturià, J., Spross, M., Hopkinson, C., Burlando, P., Lehning, M., 2013. Statistical modelling of the snow depth distribution in open alpine terrain. Hydrol. Earth Syst. Sci. 17 (8), 3005-3021. http://dx.doi.org/10.5194/hess-17-3005-2013.

Guidicelli, M., Aalstad, K., Treichler, D., Salzmann, N., 2023. A combined data assimilation and deep learning approach for continuous spatio-temporal reconstruction from Sparse Ground Tracks. http://dx.doi.org/10.2139/ssrn.4489553.

Günther, D., Marke, T., Essery, R., Strasser, U., 2019. Uncertainties in Snowpack Simulations—Assessing the impact of model structure, parameter choice, and forcing data error on point-scale energy balance snow model performance. Water Resour. Res. 55 (4), 2779-2800. http://dx.doi.org/10.1029/2018WR023403.

Hawker, L., Uhe, P., Paulo, L., Sosa, J., Savage, J., Sampson, C., Neal, J., 2022. A 30 m global map of elevation with forests and buildings removed. Env. Res. Lett. 17 (2), 024016. http://dx.doi.org/10.1088/1748-9326/ac4d4f.

He, S., Ohara, N., Miller, S.N., 2019. Understanding subgrid variability of snow depth at 1-km scale using Lidar measurements. Hydrol. Process. 33 (11), 1525-1537. http://dx.doi.org/10.1002/hyp.13415.

Helbig, N., Mott, R., Bühler, Y., Le Toumelin, L., Lehning, M., 2024. Snowfall deposition in mountainous terrain: A statistical downscaling scheme from high-resolution model data on simulated topographies. Front. Earth Sci. 11, http://dx.doi.org/10.3389/feart.2023.1308269.

Helbig, N., van Herwijnen, A., 2017. Subgrid parameterization for snow depth over mountainous terrain from flat field snow depth. Water Resour. Res. 53 (2), 1444-1456. http://dx.doi.org/10.1002/2016WR019872.

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, M., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., Thépaut, J.-N., 2020. The ERAS global reanalysis. Q. J. R. Meteorol. Soc. 146 (730), 1999-2049. http://dx.doi.org/10.1002/qj.3803.

Hoyer, S., Hamman, J., 2017. Xarray: N-D labeled arrays and datasets in Python. JORS 5 (1), 10. http://dx.doi.org/10.5334/jors.148.

Hu, W., Solz, Y., Yeligeti, M., von Bremen, L., Deng, Y., 2023. Downscaling ERA5 wind speed data: A machine learning approach considering topographic influences. Env. Res. Lett. 18 (9), 094007. http://dx.doi.org/10.1088/1748-9326/aceb0a.

Hugonnet, R., Brun, F., Berthier, E., Dehecq, A., Mannerfelt, E.S., Eckert, N., Farinotti, D., 2022. Uncertainty analysis of digital elevation models by spatial inference from stable terrains. IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens. 15, 6456-6472. http://dx.doi.org/10.1109/JSTARS.2022.3188922.

Immerzeel, W.W., Lutz, A.F., Andrade, M., Bahl, A., Biemans, H., Bolch, T., Hyde, S., Brumby, S., Davies, B.J., Elmore, A.C., Emmer, A., Feng, M., Fernández, A., Haritashya, U., Kargel, J.S., Koppes, M., Kraaijenbrink, P.D.A., Kulkarni, A.V., Mayewski, P.A., Nepal, S., Pacheco, P., Painter, T.H., Pellicciotti, F., Rajaram, H., Rupper, S., Sinisalo, A., Shrestha, A.B., Viviroli, D., Wada, Y., Xiao, C., Yao, T., Baillie, J.E.M., 2020. Importance and vulnerability of the world's water towers. Nature 577 (7790), 364-369. http://dx.doi.org/10.1038/s41586-019-1822-y.
Kartverket, 2022. Høgdedata og djupnedata [dataset].

Ketzler, G., Römer, W., Beylich, A.A., 2021. The climate of Norway. In: Beylich, A.A. (Ed.), Landscapes and Landforms of Norway. In: World Geomorphological Landscapes, Springer International Publishing, Cham, pp. 7–29.

Kim, R.S., Kumar, S., Vuyovich, C., Houser, P., Lundquist, J., Mudryk, L., Durand, M., Barros, A., Kim, E.J., Forman, B.A., Gutmann, E.D., Wrzesien, M.L., Garnaud, C., Sandells, M., Marshall, H.-P., Cristea, N., Pflug, J.M., Johnston, J., Cao, Y., Mocko, D., Wang, S., 2021. Snow ensemble uncertainty project (SEUP): Quantification of snow water equivalent uncertainty across North America via ensemble land surface modeling. Cryosphere 15 (2), 771–791. http://dx.doi.org/10.5194/tc-15-771-2021.

Kouki, K., Luojus, K., Riihelä, A., 2023. Evaluation of snow cover properties in ERA5 and ERA5-land with several satellite-based datasets in the northern hemisphere in spring 1982–2018. Cryosphere 17 (12), 5007–5026. http://dx.doi.org/10.5194/tc-17-5007-2023, URL https://tc.copernicus.org/articles/17/5007/2023/.

Lehning, M., Völksch, I., Gustafsson, D., Nguyen, T.A., Stähli, M., Zappa, M., 2006. ALPINE3D: A detailed model of mountain surface processes and its application to snow hydrology. Hydrol. Process. 20 (10), 2111–2128. http://dx.doi.org/10.1002/hyp.6204.

Li, H., Sheffield, J., Wood, E.F., 2010. Bias correction of monthly precipitation and temperature fields from intergovernmental panel on climate change AR4 models using equidistant quantile matching. J. Geophys. Res. 115 (D10), http://dx.doi.org/10.1029/2009JD012882.

Li, G., Wang, Z.S., Huang, N., 2018. A snow distribution model based on snowfall and snow drifting simulations in mountain area. JGR Atmospheres 123 (14), 7193–7203. http://dx.doi.org/10.1029/2018JD028434.

Lievens, H., Brangers, I., Marshall, H.-P., Jonas, T., Olefs, M., De Lannoy, G., 2022. Sentinel-1 snow depth retrieval at sub-kilometer resolution over the European Alps. Cryosphere 16 (1), 159–177. http://dx.doi.org/10.5194/tc-16-159-2022.

Lievens, H., Demuzere, M., Marshall, H.-P., Reichle, R.H., Brucker, L., Brangers, I., de Rosnay, P., Dumont, M., Girotto, M., Immerzeel, W.W., Jonas, T., Kim, E.J., Koch, I., Marty, C., Saloranta, T., Schöber, J., De Lannoy, G.J.M., 2019. Snow depth variability in the Northern Hemisphere mountains observed from space. Nat. Commun. 10 (1), 4629. http://dx.doi.org/10.1038/s41467-019-12566-y.

Liston, G.E., 2004. Representing subgrid snow cover heterogeneities in regional and global models. J. Clim. 17 (6), 1381–1397. http://dx.doi.org/10.1175/1520-0442(2004)017<1381:RSSCHI>2.0;CO;2.

Liston, G.E., Elder, K., 2006. A distributed snow-evolution modeling system (SnowModel). J. Hydrometeorol. 7 (6), 1259–1276. http://dx.doi.org/10.1175/JHM548.1.

Liston, G.E., Perham, C.J., Shideler, R.T., Cheuvront, A.N., 2016. Modeling snowdrift habitat for polar bear dens. Ecol. Model. 320, 114–134. http://dx.doi.org/10.1016/j.ecolmodel.2015.09.010, URL https://www.sciencedirect.com/science/article/pii/S0304380015004329.

Livneh, B., Badger, A.M., 2020. Drought less predictable under declining future snowpack. Nat. Clim. Chang. 10 (5), 452–458. http://dx.doi.org/10.1038/s41558-020-0754-8.

Lundberg, S.M., Erion, G., Chen, H., DeGrave, A., Prutkin, J.M., Nair, B., Katz, R., Himmelfarb, J., Bansal, N., Lee, S.-I., 2020. From local explanations to global understanding with explainable AI for trees. Nat Mach Intell 2 (1), 56–67. http://dx.doi.org/10.1038/s42256-019-0138-9.

Magruder, L., Neuenschwander, A., Klotz, B., 2021. Digital terrain model elevation corrections using space-based imagery and ICESat-2 laser altimetry. Remote Sens. Environ. 264, 112621. http://dx.doi.org/10.1016/j.rse.2021.112621.

Maraun, D., 2019. Statistical downscaling for climate science. Oxf. Res. Encycl. Clim. Sci. http://dx.doi.org/10.1093/acrefore/9780190228620.013.712.

Margulis, S.A., Cortés, G., Girotto, M., Durand, M., 2016. A Landsat-Era Sierra Nevada Snow Reanalysis (1985–2015). J. Hydrometeorol. 17 (4), 1203–1221. http://dx.doi.org/10.1175/JHM-D-15-0177.1.

Mazzolini, M., Aalstad, K., Alonso-González, E., Westermann, S., Treichler, D., 2024. Spatio-temporal snow data assimilation with the ICESat-2 laser altimeter. EGUSphere 1–29. http://dx.doi.org/10.5194/egusphere-2024-1404.

Mazzotti, G., Webster, C., Quéno, L., Cluzet, B., Jonas, T., 2023. Canopy structure, topography, and weather are equally important drivers of small-scale snow cover dynamics in sub-alpine forests. Hydrol. Earth Syst. Sci. 27 (11), 2099–2121. http://dx.doi.org/10.5194/hess-27-2099-2023.

Meinshausen, N., 2006. Quantile regression forests. J. Mach. Learn. Res. 7 (35), 983–999.

Melvold, K., Skaugen, T., 2013. Multiscale spatial variability of lidar-derived and modeled snow depth on Hardangervidda, Norway. Ann. Glaciol. 54 (62), 273–281. http://dx.doi.org/10.3189/2013AoG62A161.

Mendoza, P.A., Musselman, K.N., Revuelto, J., Deems, J.S., López-Moreno, J.J., McPhee, J., 2020a. Interannual and seasonal variability of snow depth scaling behavior in a Subalpine Catchment. Water Resour. Res. 56 (7), http://dx.doi.org/10.1029/2020WR027343.

Mendoza, P.A., Shaw, T.E., McPhee, J., Musselman, K.N., Revuelto, J., MacDonell, S., 2020b. Spatial distribution and scaling properties of Lidar-derived snow depth in the Extratropical Andes. Water Resour. Res. 56 (12), http://dx.doi.org/10.1029/2020WR028480.

Mortimer, C., Mudryk, L., Derksen, C., Luojus, K., Brown, R., Kelly, R., Tedesco, M., 2020. Evaluation of long-term northern hemisphere snow water equivalent products. Cryosphere 14 (5), 1579–1594. http://dx.doi.org/10.5194/tc-14-1579-2020, URL https://tc.copernicus.org/articles/14/1579/2020/.

Mott, R., Schirmer, M., Lehning, M., 2011. Scaling properties of wind and snow depth distribution in an Alpine catchment. J. Geophys. Res. 116 (D6), http://dx.doi.org/10.1029/2010JD014886.

Mott, R., Schlögl, S., Dirks, L., Lehning, M., 2017. Impact of extreme land surface heterogeneity on micrometeorology over spring snow cover. J. Hydrometeorol. 18 (10), 2705–2722. http://dx.doi.org/10.1175/JHM-D-17-0074.1.

Mott, R., Vionnet, V., Grünewald, T., 2018. The seasonal snow cover dynamics: review on wind-driven coupling processes. Front. Earth Sci. 6, http://dx.doi.org/10.3389/feart.2018.00197.

Moudrý, V., Gdulová, K., Gábor, L., Šárovcová, E., Barták, V., Leroy, F., Špatenková, O., Rocchini, D., Prošek, J., 2022. Effects of environmental conditions on ICESat-2 terrain and canopy heights retrievals in Central European mountains. Remote Sens. Environ. 279, 113112. http://dx.doi.org/10.1016/j.rse.2022.113112.

Muñoz-Sabater, J., Dutra, E., Agustí-Panareda, A., Albergel, C., Arduini, G., Balsamo, G., Boussetta, S., Choulga, M., Harrigan, S., Hersbach, H., Martens, B., Miralles, D.G., Piles, M., Rodríguez-Fernández, N.J., Zsoter, E., Buontempo, C., Thépaut, J.-N., 2021. ERA5-Land: A state-of-the-art global reanalysis dataset for land applications. Earth Syst. Sci. Data 13 (9), 4349–4383. http://dx.doi.org/10.5194/essd-13-4349-2021.

Mudryk, L., Santalaria-Otín, M., Krinner, G., Ménégoz, M., Derksen, C., Brutel-Vuilmet, C., Brady, M., Essery, R., 2020. Historical Northern Hemisphere snow cover trends and projected changes in the CMIP6 multi-model ensemble. Cryosphere 14 (7), 2495–2514. http://dx.doi.org/10.5194/tc-14-2495-2020.

Muñoz Sabater, J., 2021a. ERA5-Land hourly data from 1950 to present, copernicus climate change service (C3S) climate data store (CDS) [dataset]. http://dx.doi.org/10.24381/cds.e2161bac.

Muñoz Sabater, J., 2021b. ERA5-Land monthly data from 1950 to present, copernicus climate change service (C3S) climate data store (CDS) [dataset]. http://dx.doi.org/10.24381/cds.68d2bb30.

Muñoz-Sabater, J., Dutra, E., Agustí-Panareda, A., Albergel, C., Arduini, G., Balsamo, G., Boussetta, S., Choulga, M., Harrigan, S., Hersbach, H., Martens, B., Miralles, D.G., Piles, M., Rodríguez-Fernández, N.J., Zsoter, E., Buontempo, C., Thépaut, J.-N., 2021. ERA5-land: a state-of-the-art global reanalysis dataset for land applications. Earth Syst. Sci. Data 13 (9), 4349–4383. http://dx.doi.org/10.5194/essd-13-4349-2021, URL https://essd.copernicus.org/articles/13/4349/2021/.

Neuenschwander, A., Pitts, K., 2019. The ATL08 land and vegetation product for the ICESat-2 Mission. Remote Sens. Environ. 221, 247–259. http://dx.doi.org/10.1016/j.rse.2018.11.005.

Neuenschwander, A.L., Pitts, K., Jelley, B., Robbins, J., Klotz, B., Popescu, S.C., Nelson, R.F., Harding, D., Pederson, D., Sheridan, R., 2021. ATLAS/ICESat-2 L3A land and vegetation height, version 5.0, Boulder, Color. USA NASA Natl. Snow Ice Data Cent. Distrib. Act. Arch. Cent. http://dx.doi.org/10.5067/ATLAS/ATL08.005.

Neuenschwander, A., Pitts, K., Jelley, B., Robbins, J., Markel, J., Popescu, S., Nelson, R., Harding, D., Klotz, B., Sheridan, R., 2022. Ice, cloud, and Land Elevation Satellite 2 (ICESat-2) algorithm theoretical basis documents (ATBD) for land-vegetation along-track products (ATL08).

Niittynen, P., Luoto, M., 2018. The importance of snow in species distribution models of arctic vegetation. Ecography 41 (6), 1024–1037. http://dx.doi.org/10.1111/ecog.03348, URL arXiv:https://nsojournals.onlinelibrary.wiley.com/doi/pdf/10.1111/ecog.03348 https://nsojournals.onlinelibrary.wiley.com/doi/abs/10.1111/ecog.03348.

Oliver, M., Webster, R., 2014. A tutorial guide to geostatistics: Computing and modelling variograms and kriging. Catena 113, 56–69. http://dx.doi.org/10.1016/j.catena.2013.09.006.

Painter, T.H., Berisford, D.F., Boardman, J.W., Bormann, K.J., Deems, J.S., Gehrke, F., Hedrick, A., Joyce, M., Laidlaw, R., Marks, D., Mattmann, C., McGurk, B., Ramirez, P., Richardson, M., Skiles, S.M., Seidel, F.C., Winstral, A., 2016. The airborne snow observatory: Fusion of scanning lidar, imaging spectrometer, and physically-based modeling for mapping snow water equivalent and Snow Albedo. Remote Sens. Environ. 184, 139–152. http://dx.doi.org/10.1016/j.rse.2016.06.018.

Parr, C., Sturm, M., Larsen, C., 2020. Snowdrift landscape patterns: an arctic investigation. Water Resour. Res. 56 (12), http://dx.doi.org/10.1029/2020WR027823, e2020WR027823.

Pflug, J.M., Hughes, M., Lundquist, J.D., 2021. Downscaling snow deposition using historic snow depth patterns: diagnosing limitations from snowfall biases, winter snow losses, and interannual snow pattern repeatability. Water Resour. Res. 57 (8), http://dx.doi.org/10.1029/2021WR029999, e2021WR029999.

Porter, C., Howat, I., Noh, M.-J., Husby, E., Khuvis, S., Danish, E., Tomko, K., Gardiner, J., Negrete, A., Yadav, B., Klassen, J., Kelleher, C., Cloutier, M., Bakker, J., Enos, J., Arnold, G., Bauer, G., Morin, P., 2022. ArcticDEM - Strips, Version 4.1 [Dataset]. Harvard Dataverse, http://dx.doi.org/10.7910/dvn/c98dvs.

Revuelto, J., Billecocq, P., Tuzet, F., Cluzet, B., Lamare, M., Larue, F., Dumont, M., 2020. Random forests as a tool to understand the snow depth distribution and its evolution in mountain areas. Hydrol. Process. 34 (26), 5384–5401. http://dx.doi.org/10.1002/hyp.13951.
Revuelto, J., López-Moreno, J.I., Azorin-Molina, C., Vicente-Serrano, S.M., 2014. Topographic control of snowpack distribution in a small catchment in the central Spanish Pyrenees: Intra- and inter-annual persistence. Cryosphere 8 (5), 1989–2006. http://dx.doi.org/10.5194/tc-8-1989-2014.

Rolstad, C., Haug, T., Denby, B., 2009. Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: Application to the Western Svartisen ice cap, Norway. J. Glaciol. 55 (192), 666–680. http://dx.doi.org/10.3189/002214309789470950.

Saloranta, T.M., 2012. Simulating snow maps for Norway: Description and statistical evaluation of the seNorge snow model. Cryosphere 6 (6), 1323–1337. http://dx.doi.org/10.5194/tc-6-1323-2012.

Saloranta, T.M., 2016. Operational snow mapping with simplified data assimilation using the seNorge snow model. J. Hydrol. 538, 314–325. http://dx.doi.org/10.1016/j.jhydrol.2016.03.061.

Sentinel-2, C., 2022. Sentinel-2 MSI Level-2A BOA Reflectance [Dataset]. European Space Agency, http://dx.doi.org/10.5270/S2-znk9xsj.

Sturm, M., Liston, G.E., 2021. Revisiting the global seasonal snow classification: an updated dataset for earth system applications. J. Hydrometeorol. 22 (11), 2917–2938. http://dx.doi.org/10.1175/JHM-D-21-0070.1.

Sturm, M., Wagner, A.M., 2010. Using repeated patterns in snow distribution modeling: an arctic example. Water Resour. Res. 46 (12), http://dx.doi.org/10.1029/2010WR009434.

The pandas development team, 2024. Pandas-Dev/Pandas: Pandas v2.2.2 [Software]. Zenodo, http://dx.doi.org/10.5281/zenodo.10537285.

Tian, X., Shan, J., 2021. Comprehensive evaluation of the ICESat-2 ATL08 terrain product. IEEE Trans. Geosci. Remote Sens. 59 (10), 8195–8209. http://dx.doi.org/10.1109/TGRS.2021.3051086.

Toumelin, L.L., Gouttevin, I., Helbig, N., Galiez, C., Roux, M., Karbou, F., 2023. Emulating the adaptation of wind fields to complex terrain with deep learning. Artif. Intell. the Earth Syst. 2 (1), e220034. http://dx.doi.org/10.1175/AIES-D-22-0034.1, URL https://journals.ametsoc.org/view/journals/aies/2/1/AIES-D-22-0034.1.xml.

Treichler, D., Kääb, A., 2017. Snow depth from ICESat laser altimetry — A test study in Southern Norway. Remote Sens. Environ. 191, 389–401. http://dx.doi.org/10.1016/j.rse.2017.01.022.

Trujillo, E., Ramírez, J.A., Elder, K.J., 2007. Topographic, meteorologic, and canopy controls on the scaling characteristics of the spatial distribution of snow depth fields. Water Resour. Res. 43 (7), http://dx.doi.org/10.1029/2006WR005317.

Tryhorn, L., DeGaetano, A., 2013. A methodology for statistically downscaling seasonal snow cover characteristics over the Northeastern United States. Int. J. Climatol. 33 (12), 2728–2743. http://dx.doi.org/10.1002/joc.3626.

Tsang, L., Durand, M., Derksen, C., Barros, A.P., Kang, D.-H., Lievens, H., Marshall, H.-P., Zhu, J., Johnson, J., King, J., Lemmetyinen, J., Sandells, M., Rutter, N., Siqueira, P., Nolin, A., Osmanoglu, B., Vuyovich, C., Kim, E., Taylor, D., Merkouriadi, I., Brucker, L., Navari, M., Dumont, M., Kelly, R., Kim, R.S., Liao, T.-H., Borah, F., Xu, X., 2022. Review article: Global monitoring of snow water equivalent using high-frequency radar remote sensing. Cryosphere 16 (9), 3531–3573. http://dx.doi.org/10.5194/tc-16-3531-2022.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., Willemet, J.-M., 2012. The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2. Geosci. Model. Dev. 5 (3), 773–791. http://dx.doi.org/10.5194/gmd-5-773-2012, URL https://gmd.copernicus.org/articles/5/773/2012/.

Weiss, A., 2001. Topographic position and landforms analysis. in ESRI User Conference, San Diego, CA.

Xdem contributors, 2021. Xdem v0.0.11 [software]. http://dx.doi.org/10.5281/zenodo.4809698.

Zevenbergen, L.W., Thorne, C.R., 1987. Quantitative analysis of land surface topography. Earth Surf. Process. Landf. 12 (1), 47–56. http://dx.doi.org/10.1002/esp.3290120107.

Zhang, W., Quan, H., Srinivasan, D., 2018. Parallel and reliable probabilistic load forecasting via quantile regression forest and quantile determination. Energy 160, 810–819. http://dx.doi.org/10.1016/j.energy.2018.07.019.

Zhu, H., Liu, H., Zhou, Q., Cui, A., 2023. A XGBoost-based downscaling-calibration scheme for extreme precipitation events. IEEE Trans. Geosci. Remote Sens. 61, 1–12. http://dx.doi.org/10.1109/TGRS.2023.3294266.