# A seasonal algorithm of the snow-covered area fraction for mountainous terrain

Nora Helbig, Michael Schirmer, Jan Magnusson, Flavia Mäder, Alec van Herwijnen, Louis Quéno, Yves Bühler, Jeff Deems, Simon Gascoin

► To cite this version:

Nora Helbig, Michael Schirmer, Jan Magnusson, Flavia Mäder, Alec van Herwijnen, et al. A seasonal algorithm of the snow-covered area fraction for mountainous terrain. The Cryosphere, 2021, 15 (9), pp.4607-4624. 10.5194/tc-15-4607-2021. hal-03370375

HAL Id: hal-03370375

https://hal.science/hal-03370375

Submitted on 8 Oct 2021

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

# A seasonal algorithm of the snow-covered area fraction for mountainous terrain

Nora Helbig<fnref n="1" />, Michael Schirmer<fnref n="1" />, Jan Magnusson<fnref n="2" />, Flavia Mäder<fnref n="1" />,<fnref n="3" />, Alec van Herwijnen<fnref n="1" />, Louis Quéno<fnref n="1" />, Yves Bühler<fnref n="1" />, Jeff S. Deems<fnref n="4" />, and Simon Gascoin<fnref n="5" />

**Correspondence:** Nora Helbig (norahelbig@gmail.com)

Received: 22 December 2020 – Discussion started: 5 January 2021

Revised: 13 August 2021 – Accepted: 2 September 2021 – Published: 29 September 2021

**Abstract.** The snow cover spatial variability in mountainous terrain changes considerably over the course of a snow season. In this context, fractional snow-covered area (fSCA) is an essential model parameter characterizing how much ground surface in a grid cell is currently covered by snow. We present a seasonal fSCA algorithm using a recent scale-independent fSCA parameterization. For the seasonal implementation, we track snow depth (HS) and snow water equivalent (SWE) and account for several alternating accumulation–ablation phases. Besides tracking HS and SWE, the seasonal fSCA algorithm only requires subgrid terrain parameters from a fine-scale summer digital elevation model. We implemented the new algorithm in a multi-layer energy balance snow cover model. To evaluate the spatiotemporal changes in modeled fSCA, we compiled three independent fSCA data sets derived from airborne-acquired fine-scale HS data and from satellite and terrestrial imagery. Overall, modeled daily 1 km fSCA values had normalized root mean square errors of 7%, 12% and 21% for the three data sets, and some seasonal trends were identified. Comparing our algorithm performances to the performances of the CLM5.0 fSCA algorithm implemented in the multilayer snow cover model demonstrated that our full seasonal fSCA algorithm better represented seasonal trends. Overall, the results suggest that our seasonal fSCA algorithm can be applied in other geographic regions by any snow model application.

## 1 Introduction

In mountainous terrain, the large spatial variability in the snow cover is driven by the interaction of meteorological variables with the underlying topography. Over the course of a winter season, the dominating topographic interactions with wind, precipitation and radiation vary considerably, generating characteristic seasonal dynamics of spatial snow depth variability (e.g., Luce et al., 1999). This spatial variability, or how much of the ground is actually covered by snow, is typically characterized by the fractional snow-covered area (fSCA). The fSCA is a crucial parameter in model applications such as weather forecasts (e.g., Douville et al., 1995; Doms et al., 2011), hydrological modeling (e.g., Luce et al., 1999; Thirel et al., 2013; Magnusson et al., 2014; Griessinger et al., 2016, 2019) and avalanche forecasting (Bellaire and Jamieson, 2013; Horton and Jamieson, 2016; Vionnet et al., 2014), and it is also used for climate scenarios (e.g., Roesch et al., 2001; Mudryk et al., 2020).

The fSCA can be retrieved from various satellite sensor images, including Moderate Resolution Imaging Spectroradiometer (MODIS) or Sentinel-2 (e.g., Hall et al., 1995; Painter et al., 2009; Drusch et al., 2012; Masson et al., 2018; Gascoin et al., 2019). Nevertheless, solutions are required to correct for temporally and spatially inconsistent coverage due to time gaps between satellite revisits, data delivery and the frequent presence of clouds (Parajka and Blöschl, 2006; Gascoin et al., 2015). Though fine-scale spatial snow cover

<footnote n="1">¹WSL Institute for Snow and Avalanche Research SLF, Davos, Switzerland</footnote>

<footnote n="2">²Statkraft AS, Oslo, Norway</footnote>

<footnote n="3">³Institute of Geography, University of Bern, Bern, Switzerland</footnote>

<footnote n="4">⁴National Snow and Ice Data Center, University of Colorado, Boulder, CO, USA</footnote>

<footnote n="5">⁵Centre d'Etudes Spatiales de la Biosphère, CESBIO, Univ. Toulouse, CNES/CNRS/INRAE/IRD/UPS, 31401 Toulouse, France</footnote>
models provide spatial snow depth distributions that could be used to derive coarse-scale fSCA, applying such models to larger regions is generally not feasible. This is in part due to computational cost, a lack of detailed input data and limitations in model structure or parameters. While some of these limitations can be overcome using current snow cover model advances applying data assimilation routines (e.g., Andreadis and Lettenmaier, 2006; Nagler et al., 2008; Thirel et al., 2013; Griessinger et al., 2016; Huang et al., 2017; Baba et al., 2018; Griessinger et al., 2019; Cluzet et al., 2020), the inherent uncertainties in input or assimilation data still remain. Computationally efficient subgrid fSCA parameterizations, accounting for unresolved snow depth variability, are therefore still the method of choice for coarse-scale model systems, such as weather forecast, land surface and earth system models. Furthermore, fSCA parameterizations are essential when assimilating satellite snow-covered area data in model systems (e.g., Zaitchik and Rodell, 2009).

Several compact, closed-form fSCA parameterizations were suggested for coarse-scale model applications (e.g., Douville et al., 1995; Roesch et al., 2001; Yang et al., 1997; Niu and Yang, 2007; Su et al., 2008; Zaitchik and Rodell, 2009; Swenson and Lawrence, 2012). Some parameterizations introduced subgrid terrain parameters (e.g., Douville et al., 1995; Roesch et al., 2001; Swenson and Lawrence, 2012). The heuristic tanh form, suggested by Yang et al. (1997), was later confirmed by integrating theoretical log-normal snow distributions and fitting the resulting parametric depletion curves using the spatial snow depth distribution ($\sigma_{HS}$) in the denominator of fitted fSCA curves (Essery and Pomeroy, 2004). Through advances in remote sensing techniques, fine-scale spatial snow depth (HS) data became more readily available, allowing for the empirical parameterization of $\sigma_{HS}$ in complex topography at peak of winter (PoW) or during accumulation (Helbig et al., 2015b; Skaugen and Melvold, 2019). By parameterizing $\sigma_{HS}$ using subgrid terrain parameters, Helbig et al. (2015b) expanded the tanh fSCA parameterization of Essery and Pomeroy (2004) to account for topographic influence. Recently, Helbig et al. (2021a) re-evaluated this empirically derived fSCA parameterization with high-resolution spatial HS sets from seven different geographic regions at PoW and made it applicable across spatial scales $\ge$ 200 m by introducing a scale-dependency in the dominant model descriptors.

Many studies highlighted that the same mean HS in early winter or in late spring can lead to substantially different fSCA (Luce et al., 1999; Niu and Yang, 2007; Magand et al., 2014). This has led to the introduction of hysteresis in some fSCA parameterizations (e.g., Luce et al., 1999; Swenson and Lawrence, 2012). Previously found interannual time-persistent correlations between topographic parameters and snow depth distributions (e.g., Schirmer et al., 2011; Schirmer and Lehning, 2011; Revuelto et al., 2014; López-Moreno et al., 2017) suggest indeed that a time-dependent fSCA implementation might be feasible. However, a seasonal

model implementation of a closed form fSCA parameterization also needs to account for alternating snow accumulation and melt events during the season. Especially at lower elevations and increasingly so with climate change, the separation of the depletion curve in only one accumulation period followed by a melting period is no longer applicable (e.g., Egli and Jonas, 2009). A seasonal fSCA implementation in mountainous regions that accounts for these alternating periods is challenging. While some seasonal fSCA implementations of varying complexities were previously proposed (e.g., Niu and Yang, 2007; Su et al., 2008; Egli and Jonas, 2009; Swenson and Lawrence, 2012; Nitta et al., 2014; Magnuson et al., 2014; Riboust et al., 2019), a detailed evaluation of seasonally parameterized fSCA with fSCA derived from high-resolution spatial and temporal HS data or snow products is currently still missing.

Here, we present a seasonal fSCA implementation and evaluate it with high-resolution observation data in various geographic regions throughout Switzerland. The algorithm is based on the fSCA parameterization for complex topography proposed by Helbig et al. (2015b, 2021a). We apply two different empirical parameterizations for the spatial snow depth distribution, from Egli and Jonas (2009) and Helbig et al. (2021a), with seasonal and current HS values to describe the hysteresis. Snow accumulation and melt events during the season are accounted for by tracking the history of HS and snow water equivalent (SWE) values throughout the snow season. We implemented the algorithm in a multilayer energy balance snow cover model (modified JIM, the JULES investigation model by Essery et al., 2013) which we ran with COSMO-1 (operated by MeteoSwiss) reanalysis data, measured HS and RharesD precipitation data (MeteoSwiss). The seasonal performance of this algorithm was evaluated using fSCA data sets from terrestrial cameras, airborne surveys and satellite imagery. This allowed us to assess modeled fSCA using independent HS data sets with high spatial resolution and snow products with high temporal resolution. We further implemented the Community Land Model (CLM5.0) fSCA algorithm accounting for hysteresis in accumulation and ablation (Lawrence et al., 2018), which is based on the work of Swenson and Lawrence (2012), in the multilayer energy balance snow cover model. Modeled fSCA from the CLM5.0 fSCA algorithm was also assessed with the measured fSCA data sets and the performances compared to those of our seasonal fSCA algorithm.

## 2 Fractional snow-covered area algorithm

In the following, we introduce the seasonal fSCA algorithm in two parts. First we present the closed-form fSCA parameterization derived by Helbig et al. (2015b). This formulation uses the spatial subgrid variability in snow depth ($\sigma_{HS}$) and snow depth HS of a grid cell. To derive $\sigma_{HS}$, we used two different statistical parameterizations. Second, we describe our
seasonal fSCA algorithm, i.e., how we handle the distinctly different paths between $\sigma_{HS}$ and HS during accumulation and melt periods, i.e., the hysteresis.

## 2.1 The fSCA parameterization

The core of our seasonal algorithm is the PoW parameterization of Helbig et al. (2015b) relating fSCA to HS and $\sigma_{HS}$:

$$fSCA = \tanh \left( 1.3 \frac{HS}{\sigma_{HS}} \right). \qquad (1)$$

By including both HS and $\sigma_{HS}$, this formulation accounts for the close link between spatial subgrid snow depth variability and topography in fSCA. Although Eq. (1) was derived for PoW, in our seasonal fSCA algorithm we apply it throughout the entire snow season by using two different parameterizations for $\sigma_{HS}$, one accounting for subgrid topography (Helbig et al., 2021a), while the second only depends on HS (Egli and Jonas, 2009).

### 2.1.1 The $\sigma_{HS}$ parameterization accounting for topography

We use the PoW subgrid parameterization for $\sigma_{HS}$ in mountainous terrain originally developed by Helbig et al. (2015b) and later extended by Helbig et al. (2021a). This parameterization accounts for the impact of topography on the spatial snow depth distribution at PoW:

$$\sigma_{HS}^{\text{Helbig}} = HS^c \mu^d \exp[-(\xi/L)^2]. \qquad (2)$$

The parameterization contains two scale-dependent parameters c and d:

$$\begin{aligned} c &= 0.5330 L^{0.0389}, \\ d &= 0.3193 L^{0.1034}. \end{aligned} \qquad (3)$$

This $\sigma_{HS}$ subgrid parameterization is generally valid for domain sizes (i.e., the coarse grid cell size) $L \ge 200$ m. Besides domain size $L$, Eq. (3) requires snow depth HS and subgrid summer terrain parameters $\mu$ and $\xi$. The mean-squared-slope-related parameter $\mu = \left\{ \left[ (\partial_x z)^2 + (\partial_y z)^2 \right] / 2 \right\}^{1/2}$ is derived using partial derivatives of subgrid terrain elevations $z$, i.e., from a summer digital elevation model (DEM). The correlation length $\xi = \sqrt{2\sigma_z/\mu}$ is derived for each $L$ using the standard deviation $\sigma_z$ of terrain elevations $z$. The $L/\xi$ ratio in Eq. (3) describes the frequency of topographic features of length scale $\xi$ in a domain $L$. All terrain parameters are derived on linearly detrended summer DEMs (Helbig et al., 2015b). More details on Eqs. (2) and (3) can be found in Helbig et al. (2015b, 2021a).

### 2.1.2 The $\sigma_{HS}$ parameterization not accounting for topography

The second $\sigma_{HS}$ parameterization was developed by Egli and Jonas (2009) by fitting daily spatial HS means and standard deviation of HS from 77 weather stations distributed

throughout the Swiss Alps over six consecutive winter seasons during accumulation season. The resulting parameterization uses HS and a constant fit parameter:

$$\sigma_{HS}^{\text{Egli}} = HS^{0.839}. \qquad (4)$$

This parameterization does not account for the impact of topography on $\sigma_{HS}$.

## 2.2 Seasonal fSCA algorithm

To use the above fSCA formulation (Eq. 1) throughout an entire snow season, we track changes in HS with time. This is done to account for the fact that after a snowfall, fSCA can dramatically increase. Once the new snow has settled or started to melt, fSCA values then generally return to similar values as before. We account for this by computing two fSCA values in parallel, namely a seasonal fSCA (fSCA season) and a new snow fSCA (fSCA nsnow). The fSCA season accounts for the entire history of the snow season up to the current time step and thus all processes shaping the spatial snow depth distribution. It is therefore computed using $\sigma_{HS}^{\text{Helbig}}$ (Eq. 3), which accounts for subgrid topography. The fSCA nsnow only accounts for contributions from recent snowfall. As a snowfall generally covers most of the topography within a grid cell (i.e., all surfaces are initially covered by snow), we use $\sigma_{HS}^{\text{Egli}}$ (Eq. 4), which does not account for subgrid topography.

### 2.2.1 fSCA season

To compute fSCA season, we use extreme HS values at each time step per grid cell (Figure 1a). It is important to note that we identify these extremes using SWE rather than HS as, due to snow settlement, HS values can peak even before a precipitation event has ended. However, as our fSCA algorithm requires HS as input, we search for extreme SWE values in time and use the corresponding HS values. In the following we will not specify this anymore, and we only refer to extreme values of HS. To compute fSCA season we use $\sigma_{HS}^{\text{Helbig}}$ (Eq. 3) in the fSCA formulation (Eq. 1) as follows:

$$fSCA_{\text{season}} = \tanh \left( 1.3 \frac{HS_{\text{pseudo-min}}}{\sigma_{HS_{\text{max}}}^{\text{Helbig}}} \right). \qquad (5)$$

Here, $HS_{\text{pseudo-min}}$ is the current HS value or a recent minimum (pink dots in Figure 1a), and $\sigma_{HS_{\text{max}}}^{\text{Helbig}}$ is computed using the current seasonal maximum snow depth HS$_{max}$, i.e., the maximum in HS from the start of the season up to the current time step (green dots in Fig. 1a). We call $HS_{\text{pseudo-min}}$ a pseudo-minimum as it is not the absolute seasonal minimum. At each time step, $HS_{\text{pseudo-min}}$ and $HS_{\text{max}}$ are updated to compute fSCA. Note that after the PoW, $HS_{\text{max}}$ and $\sigma_{HS_{\text{max}}}^{\text{Helbig}}$ remain constant.

For the rare, completely flat grid cells, i.e., a subgrid mean slope angle of zero, Eq. (2) would always result in fSCA = 1.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic representation of snow depth HS extreme values used to compute fSCA for a grid cell. (a) To determine fSCA<sub>season</sub>, extremes in HS (black line) are tracked over the entire season. When HS decreases, the seasonal maximum snow depth HS<sub>max</sub> (green dots) remains constant until a new maximum is reached with subsequent snowfalls. The pseudo-minimum HS<sub>pseudo-min</sub> (pink dots) decreases when HS decreases until the next snowfall. It then remains constant until HS either exceeds HS<sub>max</sub> or decreases below the previous minimum. (b) To determine fSCA<sub>nsnow</sub>, several extremes in HS (black line) are tracked within the last 14 d (dashed black lines in a): the current value HS<sub>current</sub> (blue dot), the minimum within the last 14 d HS<sub>min</sub><sup>14d</sup> (pink dot), the maximum within the last 14 d HS<sub>max</sub><sup>14d</sup> (green dot) and the minimum prior to the most recent snowfall HS<sub>min</sub><sup>recent</sup> (yellow dot).</figcaption></figure>

In those cases, we therefore use Eq. (4) instead of Eq. (2) to compute fSCA<sub>season</sub>.

### 2.2.2 fSCA<sub>nsnow</sub>

To account for possible increases in fSCA after recent snow-falls, we evaluate fSCA (Eq. 1) using $\sigma_{HS}^{Egli}$ (Eq. 4) computed with differences in snow depth dHS (only positive changes) within the last 14 d (Fig. 1b). We use dHS rather than HS to only account for the contribution of new snow on changes in fSCA, thus as if the new snow fell on bare ground. A time window of 14 d provided reliable fSCA results after intensive testing, but the length of this period may require further investigation once more is known about changes in snow depth distributions in mountainous terrain after snowfall.

Within the 14 d time window, we compute two different fSCA values and then retain the maximum value. First, we evaluate fSCA<sub>nsnow</sub><sup>14d</sup> using the largest positive change in snow depth within the last 14 d:

$$fSCA_{nsnow}^{14d} = \tanh \left( 1.3 \frac{(HS_{current} - HS_{min}^{14d})}{\sigma_{dHS}^{Egli} HS_{min}^{14d}} \right). \quad (6)$$

Here, HS<sub>current</sub> is the snow depth at the current time step (blue dot in Figure 1b), HS<sub>min</sub><sup>14d</sup> is the minimum snow depth in the last 14 d (pink dot in Fig. 1b), and $\sigma_{dHS}^{Egli}$ is computed using the maximum difference in snow depth dHS<sup>14d</sup> = HS<sub>max</sub><sup>14d</sup> - HS<sub>min</sub><sup>14d</sup> in the last 14 d, with HS<sub>max</sub><sup>14d</sup> the maximum snow depth in the last 14 d (green dot in Fig. 1b).

Second, we evaluate fSCA<sub>nsnow</sub><sup>recent</sup> using only the most recent change in snow depth within the last 14 d:

$$fSCA_{nsnow}^{recent} = \tanh \left( 1.3 \frac{dHS_{min}^{recent}}{\sigma_{dHS}^{Egli} HS_{min}^{14d}} \right). \quad (7)$$

Here, dHS<sub>min</sub><sup>recent</sup> = HS<sub>current</sub> - HS<sub>min</sub><sup>recent</sup> is the change in snow since the most recent snowfall, where HS<sub>min</sub><sup>recent</sup> is the minimum snow depth prior to the snowfall (yellow dot in Fig. 1b). The fSCA<sub>nsnow</sub><sup>recent</sup> avoids spatial discontinuities: without this implementation, grid cells with HS > 0 m prior to a recent snowfall may have a lower fSCA value than grid cells where the same amount of new snow has fallen on the bare ground.

Finally, the maximum of fSCA<sub>nsnow</sub><sup>14d</sup> and fSCA<sub>nsnow</sub><sup>recent</sup> gives fSCA<sub>nsnow</sub> for the current time step and a grid cell.

### 2.2.3 Seasonal algorithm

Over the course of the snow season, we derive fSCA<sub>nsnow</sub> and fSCA<sub>season</sub> for each time step and grid cell (Fig. 2). The final fSCA was then obtained by taking the maximum of both values. This full seasonal fSCA algorithm, i.e., including the tracking of HS and SWE, was implemented in a distributed snow cover model. The code is publicly available on the WSL/SLF GitLab repository (see Code availability section). The data sets used to evaluate the performance of this algorithm are described in the next section.

## 3 Data

### 3.1 Modeled fSCA and HS maps

We model the snow cover evolution using the JULES investigation model (JIM). JIM is a multi-model framework of physically based energy-balance models solving the mass and energy balance for a maximum of three snow layers (Essery, 2013). While the multi-model framework JIM offers 1701 combinations of various process parameterizations, Magnusson et al. (2015) selected a specific combination that performed best for snowmelt modeling for Switzerland. The
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Illustration of modeled $f\text{SCA}_{\text{nsnow}}^{\text{recent}}$, $f\text{SCA}_{\text{nsnow}}^{14\,\text{d}}$ and $f\text{SCA}_{\text{season}}$ for one grid cell over a season. The fSCA is the maximum for each time step from $f\text{SCA}_{\text{nsnow}} = \max(f\text{SCA}_{\text{nsnow}}^{\text{recent}}, f\text{SCA}_{\text{nsnow}}^{14\,\text{d}})$ and $f\text{SCA}_{\text{season}}$. All terms are described in Sect. 2.2.</figcaption></figure>

latter model combination is used to predict daily snow mass and snowpack runoff for the operational snow hydrology service (OSHD) at WSL Institute of Snow and Avalanche Research SLF. We ran $JIM_{OSHD}$ at 1 km resolution with hourly meteorological data from the COSMO-1 model (operated by MeteoSwiss) for Switzerland. We used a reanalysis product of daily observed precipitation (RhiresD) from MeteoSwiss, as well as COSMO-1 data. Daily HS measurements from manual observers, as well as from a dense network of automatic weather stations (AWSs), were used to correct precipitation data via optimal interpolation (OI) (Magnusson et al., 2014), which is a computationally efficient data assimilation approach. Using OI in $JIM_{OSHD}$, Griessinger et al. (2019) obtained improved discharge simulations in 25 catchments over 4 hydrological years.

To describe the subgrid snow cover evolution in mountainous terrain, our seasonal fSCA algorithm was implemented in $JIM_{OSHD}$. As daily values, we used model output generated at 06:00 (UTC). In the following, modeled fSCA and HS maps refer to daily fSCA and HS from $JIM_{OSHD}$ model output.

We also computed the snow cover evolution using $JIM_{OSHD}$ with various simplifications in the seasonal fSCA algorithm, as well as with the fSCA parameterizations implemented in CLM5.0 (Lawrence et al., 2018), which are based on Swenson and Lawrence (2012) (see Table 1 for more details). This latter fSCA algorithm also accounts for heteresis in accumulation and ablation by using two different fSCA parameterizations and by tracking the seasonal maximum SWE. While subgrid topography is accounted for in the fSCA parameterization during ablation via $\sigma_z$, topography is not accounted for during snowfall events. The algorithm of Swenson and Lawrence (2012) was derived by linking daily satellite-retrieved fSCA to snow data. We implemented this algorithm in JIM using our snow tracking algorithm, i.e., the corresponding HS values such as $HS_{\text{pseudo-min}}$ (see Sect. 2.2). This was done to solely evaluate the differences in the fSCA parameterizations. In total, we performed

four additional snow cover simulations: $JIM_{\text{OSHD}}^{\text{season}}$, $JIM_{\text{OSHD}}^{\text{curr}}$, $JIM_{\text{OSHD}}^{\text{allHelbig}}$ and $JIM_{\text{OSHD}}^{\text{Swenson*}}$ (see Table 1).

## 3.2 Evaluation data

### 3.2.1 ADS fine-scale HS maps

We used fine-scale spatial HS maps gathered by airborne digital scanning (ADS) with an optoelectronic line scanner on an airplane. Data were acquired over the Wannengrat and Dischma area near Davos in the eastern Swiss Alps during winter and summer (Bühler et al., 2015). We used ADS-derived HS maps at three points in time during one winter season, namely during accumulation on 26 January (acc), at approximate peak of winter on 9 March (PoW) and during ablation season on 20 April 2016 (abl) (Marty et al., 2019). We used a summer DEM from ADS to derive the snow-free terrain parameters.

Each ADS data set covers about 150 km² at 2 m spatial resolution. Compared to TLS-derived (terrestrial laser scan) HS data, the 2 m ADS-derived HS maps had a root mean square error (RMSE) of 33 cm and a normalized median absolute deviation (NMAD) of 24 cm (Bühler et al., 2015).

### 3.2.2 ALS fine-scale HS maps

We used fine-scale spatial HS maps gathered by airborne laser scanning (ALS). The ALS campaign was a Swiss partner mission of the Airborne Snow Observatory (ASO) (Painter et al., 2016). Lidar setup and processing standards were similar to those in the ASO campaigns in California. Data were acquired over the Dischma area near Davos in the eastern Swiss Alps (see Fig. 3a in Helbig et al., 2021a). We used ALS-derived HS maps at three points in time during one winter season, namely at the approximate peak of winter on 20 March (PoW) and during the early and late ablation season on 31 March and 17 May 2017 (abl), respectively. We used a summer DEM from ALS from 29 August 2017 to derive the snow-free terrain parameters.

Each ALS data set covered about 260 km². The original 3 m resolution was aggregated to 5 m horizontal resolution. Comparing the ALS-derived HS data to manual snow probing within forest but outside canopy (i.e., not below a tree), Mazzotti et al. (2019) reported a RMSE of 13 cm and a bias of -5 cm for 20 March 2017.

### 3.2.3 Terrestrial camera images

We used camera images from terrestrial time-lapse photography in the visible band. The camera (Nikon Coolpix 5900 from 2016 to 2018, Canon EOS 400D from 2019 to 2020) was installed at the SLF/WSL in Davos Dorf in the eastern Swiss Alps (van Herwijnen and Schweizer, 2011; van Herwijnen et al., 2013). Photographs were taken of the Dorfberg in Davos, which is a large southeast-facing slope with a mean slope angle of about 30° (see Fig. 1 in Helbig et al., 2015a).
<table><caption>**Table 1.** Details of the different fSCA algorithms that are compared to the full fSCA algorithm in JIM<sub>OSHD</sub>.</caption><thead><tr><th>Algorithm name</th><th>fSCA<sub>season</sub></th><th>fSCA<sub>nsnow</sub></th><th>Tracking HS and SWE (Sect. 2.2)</th></tr></thead><tbody><tr><td>JIM<sub>OSHD</sub></td><td>Eq. (5)</td><td>Eqs. (6) and (7)</td><td>Season and 14 d</td></tr><tr><td>JIM<sub>season</sub><sub>OSHD</sub></td><td>Eq. (5)</td><td>—</td><td>Season</td></tr><tr><td>JIM<sub>current</sub><sub>OSHD</sub></td><td>tanh (1.3 <sup>HS<sub>current</sub></sup>/<sub>Helbig</sub> σ<sub>HS<sub>current</sub></sub>)</td><td>—</td><td>—</td></tr><tr><td>JIM<sub>allHelbig</sub><sub>OSHD</sub></td><td>Eq. (5)</td><td>Eqs. (6) and (7) with σ<sub>HS</sub><sup>Helbig</sup></td><td>Season and 14 d</td></tr><tr><td>JIM<sub>Swenson*</sub><sub>OSHD</sub></td><td>Eq. (8.2) in<br>Lawrence et al. (2018)</td><td>Eq. (8.1) in<br>Lawrence et al. (2018)</td><td>Season and 14 d</td></tr></tbody></table>

To obtain fSCA values from the camera images, we followed the workflow described by Portenier et al. (2020). We used the algorithm of Salvatori et al. (2011) to classify pixels in the images as snow or snow-free. Though images are taken at regular intervals (between 2 and 15 min, depending on the year), we used the image at noon to derive fSCA for that day. We evaluated images from five winter seasons (2016, 2017, 2018, 2019 and 2020) every time from 1 November to 30 June.

The resulting snow/no-snow map of the camera images had a horizontal resolution of 2 m. The field of view (FOV) overlaps with four 1 × 1 km JIM<sub>OSHD</sub> grid cells with projected visible fractions between 9 % to 40 % in each grid cell. The camera FOV covers about 0.76 km².

### 3.2.4 Sentinel-2 snow products

We used fine-scale snow-covered area maps obtained from the Theia snow collection (Gascoin et al., 2019). The satellite snow products were generated from Sentinel-2 L2A and L2B images. We used Sentinel-2 snow-covered area maps over one winter season from 20 December 2017 to 31 August 2018 for Switzerland. We further used Sentinel-2 snow maps over the Dischma area near Davos close to or on the date of the three ALS scans (18 and 28 March and 17 May 2017) and over the Dorberg area in Davos Dorf from 1 November 2017 to 30 June 2018.

The horizontal resolution of the snow product is 20 m. While the spatial coverage of the Sentinel-2 snow-covered area maps in Switzerland varies every time step, Sentinel-2 may cover several thousand square kilometers. A validation of the Theia snow product with snow depth from AWSs, through a comparison to snow maps with higher spatial resolution, as well as by visual inspection, indicated that snow is well detected, although there is a tendency to underdetect snow (Gascoin et al., 2019). The main difficulty of satellite snow products is to avoid false snow detection within clouds. Furthermore, snow omission errors may occur on steep, shaded slopes when the solar elevation is typically below 20°.

### 3.3 Derivation of 1 km fSCA evaluation data

For pre-processing, we masked out forest, rivers, glaciers or buildings in all fine-scale measurement data sets. Optical snow products that were obscured by clouds were also omitted. In all fine-scale HS data sets, we neglected HS values that were lower than 0 or above 15 m. We used a HS threshold of 0 m to decide whether or not a 2 or 5 m grid cell was snow-covered. This threshold could not be better adjusted due to a lack of independent observations.

We then aggregated all fine-scale snow data, as well as the snow products from optical imagery, in squared domain sizes L in regular grids of 1 km aligned with the OSHD model domain. For the spatial averages, we required at least 70 % valid data for the fine-scale snow data and at least 50 % valid for the satellite-derived fSCA data in each 1 km grid cell. We excluded 1 km grid cells with spatial mean slope angles larger than 60° and spatial mean measured or modeled HS < 5 cm. We further neglected 1 km grid cells with forest fractions larger than 10 %, derived from 25 m forest cover data. Overall, this led to a variable number of 1 km valid grid cells for the different data sets (Table 2). For the fine-scale snow data sets, this number ranged from 69 to 157 with a total of 668 valid 1 km grid cells. After cloud and forest removal, on average, every second day we had some valid Sentinel-2 data in Switzerland (153 valid days from the 255 calendar days). For the time period from 20 December 2017 to 31 August 2018, this resulted in 216 896 valid 1 km grid cells from a total of 2 274 991 valid OSHD grid cells in Switzerland, i.e., about 9.5 %.

These valid 1 km grid cells covered terrain elevations from 174 to 4278 m, subgrid mean slope angles from 0 to 60° and all terrain aspects. We used three of the four grid cells covered by the FOV of the terrestrial camera since one grid cell had a forest fraction larger than 10 %. On average, every fourth day we had valid camera data (337 valid days from the 1212 calendar days). Valid camera-derived fSCA for five seasons and the three grid cells covered by the FOV resulted in 931 valid 1 km grid cells from a total of 3018 valid OSHD grid cells, i.e., 31 %. The three grid cells have terrain elevations of 2077, 2168 and 2367 m and slope angles of 27, 34 and 39°. The diversity in each of the evaluation data sets af-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Probability density functions after preprocessing for all valid 1 km (a) fSCA, (b) snow depth and (c) elevation per measurement data set. All densities were normalized with the maximum in each data set. Colors represent the different measurement platforms as detailed in Sect. 3.2.</figcaption></figure>

<table><caption>**Table 2.** Details of the valid 1 km fSCA evaluation data sets after pre-processing as described in Sect. 3.3.</caption><thead><tr><th rowspan="2">Geographical region</th><th rowspan="2">Remote sensing method</th><th>Spatial coverage</th><th>Temporal coverage</th><th rowspan="2">σ<sub>fSCA</sub></th><th rowspan="2">Mean fSCA</th></tr><tr><th>[km<sup>2</sup>]</th><th>[days]</th></tr></thead><tbody><tr><td>Wannengrat and Dischma area (eastern CH)</td><td>ADS</td><td>232</td><td>3</td><td>0.05</td><td>0.98</td></tr><tr><td>Dischma and Engadin area (eastern CH)</td><td>ALS</td><td>436</td><td>3</td><td>0.08</td><td>0.96</td></tr><tr><td>Davos Dorfberg (eastern CH)</td><td>Terrestrial camera</td><td>931</td><td>337</td><td>0.23</td><td>0.81</td></tr><tr><td>Switzerland</td><td>Sentinel-2</td><td>216 896</td><td>153</td><td>0.18</td><td>0.93</td></tr></tbody></table>

ter pre-processing is indicated in Table 2 and is also shown for valid 1 km domains by means of the probability density function (pdf) for fSCA, HS and terrain elevation z in Fig. 3.

## 3.4 Performance measures

To evaluate the performance of modeled fSCA compared to the measurements, we used three measures: the root mean square error (RMSE), the normalized root mean square error (NRMSE; normalized by the mean of the measurements) and the mean percentage error (MPE; defined as measured minus modeled, normalized with the mean of the measurements).

## 4 Results

We present the evaluation of our seasonal fSCA algorithm in three sections: evaluation with fSCA derived from fine-scale HS maps near Davos, evaluation with fSCA from time-lapse photography in Davos Dorf and evaluation with fSCA from Sentinel-2 snow products over Switzerland. We further present some additional comparisons with Sentinel-2 snow products in the first two sections when Sentinel-2 data were available in the Davos area (see Sect. 3.2.4).

### 4.1 Evaluation with fSCA from fine-scale HS maps

Modeled fSCA compared well with fSCA derived from all six fine-scale HS data sets. Overall, we obtained a NRMSE of 7%, a RMSE of 0.07 and a MPE of 0.7% (Table 3). The best performance was for the two dates at the approximate PoW (NRMSE of 2%, a RMSE of 0.02 and a MPE of 0.3%), while the performance was somewhat lower during the ablation and accumulation periods.

To investigate the influence of elevation, we binned the data in 200 m elevation bands for the ADS and ALS data sets separately (Figs. 4 and 5). For ADS data, elevation-dependent modeled fSCA values were comparable to the measurements at PoW and early ablation, while the differences during accumulation were more pronounced (compare red and black dots in Fig. 4). There was also no consistent elevation trend, as during accumulation differences between modeled and measured fSCA increased with elevation, while during early ablation the opposite was true. For the ALS data, measurements were only available at PoW and during ablation. Overall, modeled fSCA values were again in line with the measurements (compare red and black dots in Fig. 5). The largest difference was observed for the lowest-elevation bin (0.15 at PoW at 1800 m; Fig. 5a), and for the late ablation data, modeled fSCA was consistently lower than ALS-derived fSCA, in particular for the lower-elevation bins (Fig. 5c).
<table><caption>**Table 3.** Performance measures for modeled fSCA with (I) fSCA derived from all fine-scale HS maps (combined ADS- and ALS-derived fSCA) and (II) Sentinel-derived fSCA (only available for ALS dates). Additionally, performance measures are shown for ALS-derived fSCA with Sentinel-derived fSCA (III) and for modeled fSCA using JIM<sup>Swenson*</sup> (IV). Given statistics are NRMSE, RMSE and MPE. For all differences we computed measured minus modeled values respectively Sentinel-derived fSCA minus ALS-derived fSCA for III. The different points in time of the season are specified in Sect. 3.2.</caption><thead><tr><th>fSCA</th><th>NRMSE</th><th>RMSE</th><th>MPE</th></tr><tr><th></th><th>[%]</th><th></th><th>[%]</th></tr></thead><tbody><tr><td colspan="4"><strong>I JIM<sub>OSHD</sub> vs. ADS&amp;ALS</strong></td></tr><tr><td>All dates</td><td>7</td><td>0.07</td><td>0.7</td></tr><tr><td>Accumulation date</td><td>8</td><td>0.08</td><td>-3.8</td></tr><tr><td>PoW dates</td><td>2</td><td>0.02</td><td>0.3</td></tr><tr><td>Ablation dates</td><td>8</td><td>0.08</td><td>1.8</td></tr><tr><td colspan="4"><strong>II JIM<sub>OSHD</sub> vs. Sentinel-2 (at ALS dates)</strong></td></tr><tr><td>All dates</td><td>9</td><td>0.08</td><td>-1.4</td></tr><tr><td>PoW dates</td><td>3</td><td>0.03</td><td>2.5</td></tr><tr><td>Ablation dates</td><td>9</td><td>0.08</td><td>-1.5</td></tr><tr><td colspan="4"><strong>III Sentinel-2 vs. ALS</strong></td></tr><tr><td>All dates</td><td>11</td><td>0.10</td><td>3.1</td></tr><tr><td>PoW date</td><td>9</td><td>0.08</td><td>-5.9</td></tr><tr><td>Ablation dates</td><td>11</td><td>0.10</td><td>3.4</td></tr><tr><td colspan="4"><strong>IV JIM<sup>Swenson*</sup><sub>OSHD</sub> vs. ADS&amp;ALS</strong></td></tr><tr><td>All dates</td><td>14</td><td>0.14</td><td>-1.2</td></tr><tr><td>Accumulation date</td><td>9</td><td>0.09</td><td>-6.1</td></tr><tr><td>PoW dates</td><td>6</td><td>0.06</td><td>-0.6</td></tr><tr><td>Ablation dates</td><td>18</td><td>0.18</td><td>-0.7</td></tr></tbody></table>

Valid Sentinel-2 data were only available on dates close to the ALS measurements (green dots in Fig. 5), not to the ADS measurement dates. Overall, modeled and Sentinel-derived fSCA values were in good agreement for the three ALS dates (II in Table 3), there was no clear elevation dependence (compare green and red dots in Fig. 5), and differences were at most 0.05 (for elevations between 2300 and 2500 m in Fig. 5c). The Sentinel-derived fSCA values can also be compared to those from the ALS scans. In this case, the performance measures were somewhat lower (compare II and III in Table 3), and Sentinel-derived fSCA values were especially lower than the ALS data in late ablation (compare green and black dots in Fig. 5c).

Our seasonal fSCA algorithm is implemented in a complex operational snow cover model framework (Sect. 3.1). Uncertainties related to input or model structure therefore impact modeled HS and fSCA values. To investigate the influence of these uncertainties more closely, we also derived two benchmark fSCA models based on Eq. (1) using measured rather than modeled HS data. The first benchmark fSCA<sup>measured</sup><sub>curr</sub> (light blue stars in Figs. 4 and 5) uses measured HS and σ<sub>HS</sub> from the current scan. The second benchmark

fSCA<sup>measured</sup><sub>PoW</sub> (orange stars in Figs. 4 and 5) combines current HS measurements with σ<sub>HS</sub> values measured at PoW. At PoW, fSCA<sup>measured</sup><sub>PoW</sub> and fSCA<sup>measured</sup><sub>curr</sub> are the same, and fSCA<sup>measured</sup><sub>PoW</sub> can only be derived at or after PoW. Results obtained with both benchmark models were similar except for the lowest-elevation bin in the ALS data set (Fig. 5b and c). Overall, the values of fSCA<sup>measured</sup><sub>curr</sub> were somewhat closer to the measured fSCA values (e.g., Figs. 4c or 5b). Both benchmark models were closest to the measured fSCA values during the ablation season (Figs. 4c and 5c), and overall the agreement was better for higher-elevation bins. Our seasonal fSCA implementation (red dots in Figs. 4 and 5) was also similar to both benchmark models. The largest differences were during the accumulation period (Fig. 4a).

As a final benchmark, we also compared our seasonal fSCA implementation with the parameterizations implemented in CLM5.0 (see Table 1). Modeled fSCA using JIM<sub>OSHD</sub> performed better than that modeled with JIM<sup>Swenson*</sup><sub>OSHD</sub> (compare I and IV in Table 3). During most of the season, fSCA values from JIM<sup>Swenson*</sup><sub>OSHD</sub> were close to 1 and showed little elevation dependence (blue stars in Figs. 4 and 5). The only exception was during the late-ablation season, when fSCA values from JIM<sub>OSHD</sub> and from JIM<sup>Swenson*</sup><sub>OSHD</sub> were very similar (red dots and dark blue stars in Fig. 5c).

To investigate the origin of the discrepancies between modeled and observed fSCA values more closely, we compared modeled and measured HS in 200 m elevation bins for the ADS and ALS data sets separately (Figs. 6 and 7). For both data sets, modeled HS was substantially lower than measured HS at higher elevations. The only exception was for the accumulation date, when modeled and measured HS values were in good agreement for all elevations (Fig. 6a). For all dates and data sets, the NRMSE between modeled and measured HS was 12%, and the MPE was 14%. Note that seasonal variations in ALS HS across all elevations were generally much lower than those in the ADS HS data. This was in part because the time intervals between the three ALS scans (20 March, 31 March and 17 May 2017) were shorter than for the ADS scans (26 January, 9 March and 20 April 2016), and there were also some snowfall events during the ALS ablation period (spring 2017).

## 4.2 Evaluation with fSCA from camera images

The high temporal resolution of camera-derived fSCA allowed us to evaluate the seasonal model performance. The seasonal trend in modeled fSCA using JIM<sub>OSHD</sub> was generally in line with that from camera-derived fSCA (compare red and black dots in Fig. 8). For the grid cell at 2168 m, however, the agreement was somewhat poorer as there was a delay in the modeled start of the ablation season, and modeled fSCA values were too high during accumulation (Fig. 8b, e).

For all winter seasons (2016 to 2020) and for the three grid cells, we obtained a NRMSE of 21%, a RMSE of 0.17 and a MPE of -7% (I in Table 4). Note that the inter-annual per-
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Modeled and ADS-derived fSCA in 200 m elevation bins for three dates: (a) during accumulation, (b) at approximate peak of winter (PoW) and (c) during ablation. Two benchmarks based on Eq. (1) are shown where applicable: $fSCA_{PoW}^{measured}$ (orange stars) uses HS form from the current ADS scan and $\sigma_{HS}$ from the ADS scan at PoW, while $fSCA_{curr}^{measured}$ (light blue stars) uses HS and $\sigma_{HS}$ from the current ADS scan. The bars show the valid data percentage per bin.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Modeled and ALS-derived, as well as Sentinel-derived, fSCA in 200 m elevation bins for three dates: (a) at approximate PoW, (b) during early ablation and (c) during late ablation. The same two benchmarks based on Eq. (1) as in Fig. 4 are also shown where applicable. Sentinel-derived fSCA (green dots) was available 2 d before the PoW scan, 3 d before the early ablation scan and on the same day as the late ablation scan. The bars show the valid data percentage per bin.</figcaption></figure>

performance varied substantially, as did the performance among the three grid cells. For instance, for all three grid cells, the overall best performance was for the season 2018 (NRMSE = 14%, RMSE = 0.11, MPE = -4%), while the worst performance was for the season 2019 (NRMSE = 25%, RMSE = 0.2, MPE = -12%).

For winter season 2018, we used Sentinel-derived fSCA to evaluate modeled and camera-derived fSCA values. While overall the agreement between modeled and Sentinel-derived fSCA was good (NRMSE 2% and MPE of 1%), the agreement between camera- and Sentinel-derived fSCA was poorer (NRMSE = 12%, MPE = 5%). The latter performance values were, however, comparable to the agreement between modeled and camera-derived fSCA for days with valid Sentinel-derived data (NRMSE = 12%, MPE = -4%).

The camera-derived fSCA was also used to evaluate the relevance of applying our full seasonal fSCA algorithm as opposed to simplifications and $JIM_{OSHD}^{Swenson*}$ (see Table 1 for details). While overall fSCA from $JIM_{OSHD}^{season}$ and $JIM_{OSHD}$

agreed well, there were substantial differences after snow-fall events on partly snow-free ground (compare orange stars and red dots in Fig. 8). Specifically, after such a snowfall event, modeled fSCA using $JIM_{OSHD}$ generally increased, while $JIM_{OSHD}^{season}$ remained constant. Using $JIM_{OSHD}^{curr}$, modeled fSCA values were less in line with those from $JIM_{OSHD}$ (compare light blue stars and red dots in Fig. 8). While discrepancies were again large after snowfall events, they were also pronounced during the ablation periods. In general, with $JIM_{OSHD}^{curr}$ the ablation season started later and was followed by a much steeper melt-out period. Using $JIM_{OSHD}^{curr}$ can result in a substantially shorter snow season compared to $JIM_{OSHD}$, with a maximum difference of 21 d at 2168 m in the season 2017. Overall, compared to camera-derived fSCA, both simplified models performed less well than $JIM_{OSHD}$ (Table 4). The performance using $JIM_{OSHD}^{allHelbig}$ was very similar to fSCA from $JIM_{OSHD}$; i.e., applying $\sigma_{HS}^{Helbig}$ instead of $\sigma_{HS}^{Egli}$ for fSCA$_{nsnow}$ did not substantially affect model performance. On the contrary, fSCA from $JIM_{OSHD}^{Swenson*}$ had the worst overall performances when compared to camera-derived fSCA
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Modeled and ADS-derived HS in 200 m elevation bins for three dates: (a) during accumulation, (b) at approximate PoW and (c) during ablation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Modeled and ALS-derived HS in 200 m elevation bins for three dates: (a) at approximate PoW, (b) during early ablation and (c) during ablation.</figcaption></figure>

(VII in Table 4). Similar to $JIM_{OSHD}^{cur}$, using $JIM_{OSHD}^{Swenson*}$ considerably delayed the ablation season, followed by a much steeper melt out. The snow season was substantially shortened again by at most 32 d in the 2017 season at 2077 m. Modeled fSCA using $JIM_{OSHD}^{Swenson*}$ also largely overestimates fSCA during the accumulation period (blue dots in Fig. 8). Overall, using $JIM_{OSHD}^{Swenson*}$ led to much steeper increases and decreases in fSCA, i.e., an almost binary seasonal fSCA trend that was not in line with camera-derived fSCA.

## 4.3 Evaluation with fSCA from Sentinel-2 snow products

Overall, modeled fSCA using $JIM_{OSHD}$ compared well with Sentinel-derived fSCA throughout the season (I in Table 5). To investigate the elevation-dependent differences between modeled and Sentinel-derived fSCA in more detail, we binned the data in 250 m elevation bands for each day throughout the entire season (Fig. 9). To estimate the end of the accumulation (1 April 2018) and ablation season (30 June 2018), we used the spatial mean HS (solid black line at bottom of Fig. 9). Overall, differences in performance between the accumulation and the ablation period were small (I in Table 5). However, there were marked differences with el-

evation throughout the season. Up to the end of the accumulation period, the largest differences between modeled and Sentinel-derived fSCA were at elevations lower than 1500 m, whereas at elevations above around 3000 m, the agreement was good (Fig. 9a). During the ablation period, most of the snow at lower elevations was gone, and modeled fSCA was generally larger than Sentinel-derived fSCA at higher elevations (> 2500 m), in particular towards the end of the ablation season. During the summer (30 June to 30 August 2018), i.e., after the end of the ablation season, modeled fSCA was larger than Sentinel-derived fSCA at the highest elevations (> 3500 m), whereas between the snow line and these highest elevations, modeled fSCA was generally lower.

Given the high temporal resolution of the Sentinel-derived fSCA data set, we again evaluated the fSCA algorithm simplifications and $JIM_{OSHD}^{Swenson*}$ (see Table 1). Compared to our seasonal implementation, the overall performance values of the fSCA algorithm simplifications were similar except for $JIM_{OSHD}^{cur}$ and $JIM_{OSHD}^{Swenson*}$ (Table 5). Modeled fSCA values with $JIM_{OSHD}^{cur}$ and $JIM_{OSHD}^{Swenson*}$ were generally larger than Sentinel-derived fSCA, resulting in larger MPE values with the largest ones for $JIM_{OSHD}^{Swenson*}$ (compare I, III and V in Table 5). This is also clearly reflected in the elevation-
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Modeled and camera- and Sentinel-derived fSCA for the three 1 km grid cells within the field of view of the camera for two seasons: (a–c) winter 2017 and (d–f) winter 2018.</figcaption></figure>

dependent differences between fSCA using $JIM_{OSHD}^{Swenson*}$ and
Sentinel-derived fSCA throughout the season (Fig. 9b).

# 5 Discussion

## 5.1 Fractional snow-covered area fSCA algorithm

Our seasonal fSCA algorithm is based on the closed-form
fSCA parameterization of Helbig et al. (2015a) (Eq. 1) and
combines two statistical parameterizations for $\sigma_{HS}$, together
with a tracking method, to account for changes in maxi-
mum snow depth and precipitation events. The algorithm is
modular, meaning that individual parts can easily be com-
plemented or replaced with new parameterizations, e.g., for
fSCA<sub>nsnow</sub>. Overall, our algorithm only requires subgrid cell
summer terrain parameters, which are a slope-related param-
eter and the terrain correlation length, and tracking snow in-
formation.

We evaluated the performance of our seasonal fSCA im-
plementation in Switzerland. We could not explicitly eval-
uate the performance for completely flat grid cells, i.e., grid
cells with a subgrid mean slope angle of zero. After removing
rivers/lakes, we only had five 1 km grid cells for Switzerland
with a subgrid mean slope angle of zero, i.e., 0.01 % of all
grid cells. For these grid cells, using $\sigma_{HS}^{Helbig}$ (Eq. 2) always
results in a fSCA of 1. As a first approach, we therefore pro-
posed to use $\sigma_{HS}^{Egli}$ (Eq. 4). Although we see no reason why
our fSCA algorithm could not be used in other geographic
region, it remains unclear at this point if our seasonal fSCA
implementation can also be used in flat regions.

We used $\sigma_{HS}^{Egli}$ (Eq. 4), which does not account for subgrid
topography, to derive fSCA<sub>nsnow</sub>. We did this to account for
uniform blanketing after a snowfall, i.e., to account for pos-
sible increases in fSCA after a recent snowfall. When sub-
stituting $\sigma_{dHS}^{Egli}$ with $\sigma_{dHS}^{Helbig}$ in Eqs. (6) and (7) ($JIM_{OSHD}^{allHelbig}$;
see Table 1), the overall performance was very similar (Ta-
bles 4 and 5). Thus, while applying $\sigma_{dHS}^{Egli}$ might not describe
the true spatial new snow distribution in mountainous ter-
rain, the formulation is simple and is therefore used here as
a first approach. Based on the modular algorithm setup, dif-
ferent closed-form fSCA parameterizations can be applied
in our seasonal algorithm, e.g., for a flat grid cell or for
fSCA<sub>nsnow</sub> (for some empirical examples, see Essery and
Pomeroy, 2004).

## 5.2 Evaluation

### 5.2.1 Evaluation with fSCA from fine-scale HS maps

The evaluation of the seasonal fSCA algorithm with fSCA
from fine-scale HS maps showed that overall the model per-
formed well, especially at PoW (I in Table 3). Modeled
fSCA using $JIM_{OSHD}^{Swenson*}$, on the other hand, generally over-
estimated fSCA (MPE< 0). This algorithm intercomparison
shows that the seasonal fSCA evolution is better captured by
$JIM_{OSHD}$ most likely because the $JIM_{OSHD}^{Swenson*}$ model does
not sufficiently account for the high spatial variability in
snow distribution in complex terrain.

During accumulation at higher elevations, modeled fSCA
using $JIM_{OSHD}$ overestimated ADS-derived fSCA even
though modeled HS agreed reasonably well with the mea-
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Difference between Sentinel-derived and modeled fSCA for Switzerland as a function of date and elevation z (in 250 m elevation bins) for available satellite dates for (a) JIM<sub>OSHD</sub> and (b) JIM<sub>OSHD</sub><sup>Swenson*</sup>. Daily spatial mean snow depth HS is also shown (solid black line). The vertical lines indicate the dates for the end of accumulation (dashed) and ablation (line with stars) seasons.</figcaption></figure>

<table><caption>**Table 4.** Performance measures for (I) modeled fSCA using JIM<sub>OSHD</sub> and camera-retrieved fSCA for the winter seasons 2016 to 2020, (II) modeled fSCA using JIM<sub>OSHD</sub> and Sentinel-derived fSCA for the three grid cells for the winter season 2018, (III) camera-derived fSCA with Sentinel-derived fSCA for the three grid cells, and (IV to VII) all JIM-modeled fSCA versions (for details see Table 1), namely for JIM<sub>OSHD</sub><sup>season</sup>, JIM<sub>OSHD</sub><sup>curr</sup>, JIM<sub>OSHD</sub><sup>allHelbig</sup> and JIM<sub>OSHD</sub><sup>Swenson*</sup>, with camera-derived fSCA.</caption><thead><tr><th>fSCA</th><th>NRMSE</th><th>RMSE</th><th>MPE</th></tr><tr><td></td><td>[%]</td><td></td><td>[%]</td></tr></thead><tbody><tr><td>I JIM<sub>OSHD</sub> vs. camera</td><td></td><td></td><td></td></tr><tr><td></td><td>21</td><td>0.17</td><td>-7.1</td></tr><tr><td>II JIM<sub>OSHD</sub> vs. Sentinel-2</td><td>2</td><td>0.02</td><td>0.8</td></tr><tr><td>III Camera vs. Sentinel-2</td><td>12</td><td>0.11</td><td>5.0</td></tr><tr><td>IV JIM<sub>OSHD</sub><sup>season</sup> vs. camera</td><td>22</td><td>0.18</td><td>-6.1</td></tr><tr><td>V JIM<sub>OSHD</sub><sup>curr</sup> vs. camera</td><td>26</td><td>0.21</td><td>-9.2</td></tr><tr><td>VI JIM<sub>OSHD</sub><sup>allHelbig</sup> vs. camera</td><td>21</td><td>0.17</td><td>-7.6</td></tr><tr><td>VII JIM<sub>OSHD</sub><sup>Swenson*</sup> vs. camera</td><td>30</td><td>0.25</td><td>-10.6</td></tr></tbody></table>

surements (Figs. 4a and 6a). We also used a different model configuration (JIM<sub>OSHD</sub><sup>allHelbig</sup> in Table 1), yet fSCA values did not substantially change for the accumulation date (not shown). Based on this we assume that both σ<sub>HS</sub> parameterizations cannot sufficiently describe snow redistribution during accumulation likely due to periods with strong winds following snowfall. The description of σ<sub>HS</sub> during the accumulation period thus needs to be improved. This will, however, require more than one spatial HS data set during accumulation.

At PoW and during the ablation season, JIM<sub>OSHD</sub> mostly underestimated fSCA compared to fSCA from fine-scale HS maps, without a clear elevation trend (Figs. 4 and 5). Discrepancies between modeled and measured HS, on the other hand, generally increased with elevation (Figs. 6 and 7). Obviously for larger snow depth, correctly modeling HS has little effect on fSCA. The overall underestimated modeled fSCA values were likely a consequence of the HS threshold of 0 m we used to decide whether a 2 or 5 m grid cell was snow-covered or not. In reality, due to measurement uncertainties, both small positive or negative measured HS values can still be associated with snow-free areas. When arbitrarily increasing the HS threshold to ±10 cm for the ALS data, modeled 1 km fSCA values were rather larger than the measurements (not shown). This is not contradictory but emphasizes the need to accurately model HS along snow lines, in which small inaccuracies in HS can have large impacts on fSCA. For instance, during early ablation modeled and measured fSCAs are larger in the lowest-elevation bin than at higher elevations (see Fig. 4c). Unfortunately, we currently do not have detailed snow observations available to define robust HS threshold values which take into account the different points in time of the season, as well as the influence of terrain and ground cover. However, the overall good agreement between Sentinel- and ALS-derived fSCA (Fig. 5 and
<table><caption>**Table 5.** Performance measures (I) for modeled fSCA using JIM<sub>OSHD</sub> and Sentinel-retrieved fSCA for the winter season 2018 for all valid 1 km grid cells of Switzerland and for all dates (20 December 2017 to 30 June 2018), for the accumulation period (20 December to 1 April), and for the ablation period (1 April to 30 June), as well as (II to V) for all JIM-modeled fSCA versions (for details, see Table 1), namely for JIM<sub>OSHD</sub>, JIM<sub>OSHD</sub><sup>season</sup>, JIM<sub>OSHD</sub><sup>curr</sup>, JIM<sub>OSHD</sub><sup>allHelbig</sup> and JIM<sub>OSHD</sub><sup>Swenson*</sup>.</caption><thead><tr><th>fSCA vs. Sentinel-2</th><th>NRMSE</th><th>RMSE</th><th>MPE</th></tr><tr><th></th><th>[%]</th><th></th><th>[%]</th></tr></thead><tbody><tr><td colspan="4"><strong>I JIM<sub>OSHD</sub></strong></td></tr><tr><td>All dates</td><td>12</td><td>0.11</td><td>0.4</td></tr><tr><td>Accumulation period</td><td>11</td><td>0.11</td><td>0.3</td></tr><tr><td>Ablation period</td><td>14</td><td>0.12</td><td>0.5</td></tr><tr><td colspan="4"><strong>II JIM<sub>OSHD</sub><sup>season</sup></strong></td></tr><tr><td>All dates</td><td>12</td><td>0.12</td><td>0.4</td></tr><tr><td>Accumulation period</td><td>11</td><td>0.11</td><td>0.3</td></tr><tr><td>Ablation period</td><td>14</td><td>0.12</td><td>0.5</td></tr><tr><td colspan="4"><strong>III JIM<sub>OSHD</sub><sup>curr</sup></strong></td></tr><tr><td>All dates</td><td>14</td><td>0.13</td><td>-0.8</td></tr><tr><td>Accumulation period</td><td>11</td><td>0.11</td><td>0.1</td></tr><tr><td>Ablation period</td><td>18</td><td>0.16</td><td>-2.4</td></tr><tr><td colspan="4"><strong>IV JIM<sub>OSHD</sub><sup>allHelbig</sup></strong></td></tr><tr><td>All dates</td><td>12</td><td>0.11</td><td>0.3</td></tr><tr><td>Accumulation period</td><td>11</td><td>0.11</td><td>0.2</td></tr><tr><td>Ablation period</td><td>14</td><td>0.12</td><td>0.5</td></tr><tr><td colspan="4"><strong>V JIM<sub>OSHD</sub><sup>Swenson*</sup></strong></td></tr><tr><td>All dates</td><td>18</td><td>0.17</td><td>-1.8</td></tr><tr><td>Accumulation period</td><td>17</td><td>0.16</td><td>-0.7</td></tr><tr><td>Ablation period</td><td>21</td><td>0.19</td><td>-3.6</td></tr></tbody></table>

III in Table 3) provides some confidence in the fine-scale HS data-derived fSCA used here to evaluate modeled fSCA.

The two benchmark fSCA models based on Eq. (1) using measured rather than modeled HS data (fSCA<sup>measured</sup><sub>curr</sub> and fSCA<sup>measured</sup><sub>PoW</sub>) generally showed similar trends as HS-derived and modeled fSCA (Figs. 4 and 5). At PoW, fSCA<sup>measured</sup><sub>curr</sub> agreed less well with measured fSCA than our seasonal implementation (see Figs. 4b and 5a). This may indicate uncertainties in the empirical fSCA parameterization (Eq. 1), which requires the further investigation of spatial HS data sets during accumulation. During ablation, we expected that fSCA<sup>measured</sup><sub>PoW</sub> would be closer to measured fSCA than fSCA<sup>measured</sup><sub>curr</sub>, which was, however, not the case (see Figs. 4c and 5b). Since the true PoW date is elevation and aspect dependent, we cannot assume that one date for PoW is representative for the entire catchment, covering several hundred square kilometers and large elevation gradients. Thus, mea-

sured σ<sub>HS</sub> at the date we defined as PoW might not have been representative for the true σ<sub>HSmax</sub> in each grid cell as required by Eq. (5). Besides possible uncertainties in the empirical fSCA parameterization (Eq. 1), we assume this is the main reason why these two benchmark models using measured HS data did not outperform our seasonal implementation. Overall, these comparisons emphasize the need for tracking snow information per grid cell, as is done by our seasonal fSCA algorithm.

### 5.2.2 Evaluation with camera-derived fSCA

The evaluation with fine-scale HS maps revealed overall good model performance at six points in time. It was, however, not possible to comprehensively evaluate the performance over the season. For this, we used daily camera-derived fSCA, showing that the modeled seasonal fSCA trend was mostly in line with observations (Fig. 8).

Model performance compared to the camera-derived fSCA values was overall worse than when comparing to HS-derived fSCA (e.g., NRMSE of 21% for I in Table 4 compared to NRMSE of 7% for I in Table 3). Since the higher temporal resolution of the camera data set leads to the largest spread in fSCA values compared to the other two data sets (see Table 2 and Fig. 3), a larger portion of intermediate fSCA values (e.g., close to the snow line) are included which are generally more difficult to model correctly than fSCA values close to 1. The poorer model performance is, however, likely also due to the overall lower accuracy of camera-derived fSCA. For instance, the projection of the 2D camera image to a 3D DEM may introduce errors and distortions. Furthermore, when deriving fSCA from camera images, clouds/fog and uneven illumination, for instance due to shading or partial cloud cover, may deteriorate the accuracy (e.g., Farinotti et al., 2010; Fedorov et al., 2016; Härer et al., 2016; Portenier et al., 2020). Another factor affecting the performance measures was the threshold for the number of valid fine-scale data per 1 km grid cell. When aggregating to 1 km fSCA maps for the Sentinel-derived values, we required at least 50% valid fine-scale data. This requirement could not be met for camera-derived fSCA as the projected fractions of the camera FOV on the 1 km model grid cells were only 9%, 13% and 14%. This is reflected in the better agreement between modeled and Sentinel-derived fSCA than between camera- and Sentinel-derived fSCA (NRMSE of 2% vs. 12% in Table 4). Finally, as the camera was installed at valley bottom, steep slope sections cover larger areas of the FOV, while flatter slope parts remain invisible. This likely led to underestimated fSCA values. On the other hand, valid Sentinel-derived fSCA has a much lower temporal resolution and did not cover the entire ablation period. Instead, Sentinel-derived fSCA was often available throughout the period when fSCA was rather close to 1 (see Fig. 8d, e). Thus, while there is likely more uncertainty in camera-derived fSCA, the high temporal resolution of this product
still provides valuable information on model performance
throughout the season.

We used the camera-derived fSCA to also evaluate sim-
plifications of our seasonal fSCA algorithm, as well as
JIM<sup>Swenson*</sup><sub>OSHD</sub> (Table 1). Compared to our seasonal fSCA im-
plementation, the more simple implementations did not cap-
ture the seasonal variation as well (Fig. 8). With JIM<sub>OSHD</sub><sup>curr</sup>,
the start of the ablation season was delayed, and the abla-
tion season was also considerably shortened by up to 21 d.
In this respect, the results for JIM<sub>OSHD</sub><sup>Swenson*</sup> were very simi-
lar as overall the increases and decreases in fSCA were very
steep, leading to shortened snow seasons and poorer perfor-
mances (see Table 4). In principle, JIM<sub>OSHD</sub><sup>curr</sup> considers each
day as PoW, leading to rapid changes in fSCA, in particu-
lar when HS values are low (i.e., early accumulation or ab-
lation season). In JIM<sub>OSHD</sub><sup>season</sup>, the seasonal maximum value of
HS was additionally tracked, substantially improving the sea-
sonal fSCA trend, in particular during the ablation season.
However, changes in fSCA due to snowfall events were still
not captured well with this implementation, showing that our
new snow tracking algorithm further improves the overall
model performance. Since the impact of using JIM<sub>OSHD</sub><sup>allHelbig</sup>
on modeled fSCA is mainly restricted to snowfall follow-
ing melt periods, overall performances were very similar to
JIM<sub>OSHD</sub> (see Tables 4 and 5). This again indicates that the
description of σ<sub>HS</sub> following snowfall events requires further
investigation.

5.2.3 Evaluation with Sentinel-derived fSCA

By including Sentinel-derived fSCA in our evaluation, we
added a data set with both a high temporal resolution and
a much larger spatial coverage (see Table 2). The Sentinel-
derived fSCA data set comprised about 217 000 1 km grid
cells covering a wide range of terrain elevations, slope angles
and terrain aspects.

For the investigated winter season, results showed an
overall good seasonal agreement across Switzerland, though
there was some elevation-dependent scatter (Fig. 9a). Dis-
crepancies during accumulation occurred mostly along the
snowline at lower elevations, where lower spatial HS values,
as well as more cloudy weather, prevail during accumulation.
Both can lead to inaccurate modeled and Sentinel-derived
fSCA. Furthermore, we assume that some of the overestima-
tions in modeled fSCA at higher elevations during accumu-
lation could also stem from underestimated σ<sub>HS</sub> during peri-
ods when strong winds follow snowfall events, as was also
observed in the HS data sets (Fig. 4a and Sect. 5.2.1). The
scatter at high elevations during ablation and summer likely
originates from lower modeled fSCA due to underestimated
precipitation as there are fewer AWSs at high elevations for
data assimilation in our model.

Performance measures were somewhat poorer than those
from fine-scale HS maps (e.g., NRMSE of 12% for Sentinel
vs. 7% for fSCA for HS data). Uncertainties introduced by

reduced visibility in the snow products of Sentinel-2 are the
most likely reason for this. Both our camera and the Sentinel-
2 data sets cover long time periods at higher temporal res-
olution; i.e., they include also periods under unfavorable
weather conditions. On the contrary, clear sky dates were
carefully selected for the on-demand high-quality data acqui-
sitions from the air for our fSCA data sets derived from fine-
scale HS maps. Nevertheless, the camera and the Sentinel-2
data sets enabled us to evaluate seasonal fSCA model trends
which would not have been possible from only six fSCA data
sets derived from HS data.

When evaluating the simplified fSCA algorithms and
JIM<sup>Swenson*</sup><sub>OSHD</sub>, model performance measures were compar-
ble to our seasonal implementation except for JIM<sub>OSHD</sub><sup>curr</sup> and
JIM<sup>Swenson*</sup><sub>OSHD</sub> (Table 5), as was also the case for the compari-
son with camera-derived fSCA (Table 4). For Sentinel- and
camera-derived fSCA, the main reason is likely the limited
availability of fSCA data during or shortly after snowfall due
to bad visibility and clouds. Additionally, for the Sentinel-
derived fSCA, local performance differences across Switzer-
land are likely averaged out. Nevertheless, fSCA values when
using JIM<sup>Swenson*</sup><sub>OSHD</sub> were overestimated compared to Sentinel-
derived values (Fig. 9b, and negative MPE for V in Table 5).
Similar results were also observed when using JIM<sub>OSHD</sub><sup>curr</sup> (see
negative MPE for III in Table 5). These biases are most likely
related to the rather steep increases and decreases in mod-
eled fSCA over the season, as we also observed with the
camera-derived fSCA (Fig. 8). We further assume that over-
estimated fSCA using JIM<sup>Swenson*</sup><sub>OSHD</sub> at higher elevations due
to underestimating spatial snow depth variability in complex
terrain may have compensated for other modeled fSCA er-
ror sources (e.g., from underestimated precipitation input at
these elevations), leading to an overall lower bias at higher el-
evations during accumulation compared to our fSCA imple-
mentation. Finally, note that the scatter above zero between
Sentinel-derived and JIM<sup>Swenson*</sup><sub>OSHD</sub> fSCA (Fig. 9b) almost dis-
appears when we neglect all 1 km domains with modeled
HS < 5 cm using JIM<sub>OSHD</sub><sup>Swenson*</sup> (not shown). While the over-
all NRMSE values for JIM<sup>Swenson*</sup><sub>OSHD</sub> are then comparable to
our seasonal implementation (e.g., NRMSE of 12% for all
dates instead of 18%; see V in Table 5), it reveals the overall
overestimation of JIM<sup>Swenson*</sup><sub>OSHD</sub> (e.g., increased negative MPE
of -4.1% for all dates instead of -1.8%). Clearly, our sea-
sonal fSCA implementation is better suited to more realisti-
cally represent seasonal changes in mountainous terrain, in
particular following snowfall and during the ablation period.

6 Conclusions

We presented a seasonal fractional snow-covered area
(fSCA) algorithm based on the fSCA parameterization of
Helbig et al. (2015b, 2021a). The seasonal algorithm is based
on tracking HS and SWE values accounting for alternating
snow accumulation and melt events. Two empirical parame-
terizations were used to describe the spatial snow depth dis-
tribution: one for mountainous terrain and one not accounting
for subgrid topography. An implementation in a multilayer
energy balance snow cover model system (JIMOSHD; JIM,
JULES investigation model; Essery et al., 2013) allowed us
to evaluate seasonally modeled fSCA for Switzerland.

Compiling independent fSCA data sets with different spa-
tiotemporal characteristics enabled a thorough analysis of
the seasonal fSCA algorithm in mountainous terrain of daily
1 km fSCA values. While the evaluation with the three data
sets showed overall good seasonal performance, each of
the evaluation data sets allowed specific conclusions to be
drawn. The evaluation with fine-scale spatial HS-derived
fSCA showed that HS uncertainties along the snow line
likely contributed most to the underestimation of fSCA dur-
ing ablation and PoW, emphasizing the need to accurately
model HS along snow lines. The camera-derived fSCA data
set with the highest temporal resolution confirmed the need
for tracking HS over the season, as well as accounting for in-
termediate snowfalls, to avoid a delayed melt start and a dras-
tic shortening of the ablation season. The Sentinel-derived
fSCA data set, with the largest spatial coverage, together
with a rather high temporal resolution, demonstrated that the
seasonal fSCA algorithm performs well across a range of
elevations, slope angles, terrain aspects and snow regimes.
This comparison showed that there were some differences
at low elevation or along the snowline coinciding with low
HS, while discrepancies occurred mostly at high elevations
towards the end of the season during summer.

Overall, NRMSEs for seasonally modeled fSCA increased
from 7% for HS data-derived fSCA to 12% for Sentinel-
derived fSCA and to 21% for camera-derived fSCA. While
the large variation in performance measures is likely tied to
the various temporal and spatial resolutions of the data sets
and measurement uncertainties, it also demonstrates the dif-
ficulties in drawing conclusions when evaluating a model al-
gorithm with evaluation data from different acquisition plat-
forms. Nevertheless, this comparison with data covering a
wide range of spatiotemporal scales allowed us to obtain a
comprehensive overview of the strengths and weaknesses of
our seasonal fSCA implementation. We are not aware of any
seasonal fSCA implementation that has been evaluated in
such detail by exploiting independent HS and snow product
data sets at high spatial and temporal resolution.

By implementing the fSCA parameterizations applied in CLM5.0 (Lawrence et al., 2018) in JIM<sub>OSHD</sub>, we also evaluated modeled fSCA using <i>JIM</i><sup>Swenson*</sup>. This showed that our seasonal fSCA algorithm captures the seasonal variation best and that seasonal variation in <i>JIM</i><sup>Swenson*</sup> was limited. <i>JIM</i><sub>OSHD</sub><sup>Swenson*</sup> resulted in often overestimated fSCA values likely because the high spatial variability in snow depth distribution in complex terrain is not sufficiently described.

The implementation of the seasonal fSCA algorithm in a
model only requires subgrid terrain parameters from a fine-
scale summer DEM in combination with tracking HS and

SWE for coarse grid cells. The algorithm is set up such that improvements or adaptations of individual algorithm parts can easily be implemented. The PoW fSCA parameterization of Helbig et al. (2015b) forms the centerpiece of the presented seasonal fSCA algorithm. The recent re-evaluation with various spatial PoW snow depth data sets from seven geographic regions showed an overall NRMSE of only 2% (Helbig et al., 2021a). This detailed evaluation at PoW in different geographic regions, together with the seasonal assessment with the three fSCA data pools presented here, suggests that the seasonal fSCA algorithm may also be used in other geographic regions. However, further investigations, once more spatial HS data sets before and after snowfalls in complex topography become available, would be advantageous for improvements of our seasonal fSCA algorithm, especially during the accumulation period.

***Code availability.*** The code of the full algorithm is made avail-
able on the WSL/SLF GitLab repository with the link and meta-
data available on the Envidat data repository of the Swiss Federal
Institute WSL (https://doi.org/10.16904/envidat.244, Helbig et al.,
2021b).

***Data availability.*** All data used in this study are described in the data section. The data can be downloaded from the referenced repositories, or data availability is described in the referenced publications. Theia snow maps are freely distributed via the Theia portal (https://doi.org/10.24400/329360/F7Q52MNK, Gascoin et al., 2019).

***Author contributions.*** The seasonal algorithm was developed by JM, NH and MS. LQ set up and ran the seasonal JIM simulations. The evaluation data sets were designed, gathered and pre-processed by YB, SG, FM, AvH and JSD. The design and performance of the evaluation of the algorithm with measurement data were conducted by NH with discussions and input from MS, AvH, SG and YB. NH had the lead in writing the manuscript with contributions from all co-authors.

***Competing interests.*** The authors declare that they have no conflict of interest.

***Disclaimer.*** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

***Acknowledgements.*** We thank Andreas Stoffel at SLF for his help with GIS processing of the satellite images.
*Financial support.* Nora Helbig was funded by a grant of the Swiss National Science Foundation (SNF) (grant no. IZSEZ_186887), as well as being partly funded by the Federal Office of the Environment FOEN. Simon Gascoin acknowledges the support of the Centre national d'études spatiales (CNES).

*Review statement.* This paper was edited by Carrie Vuyovich and reviewed by two anonymous referees.

## References

Andreadis, K. M. and Lettenmaier, D. P.: Assimilating remotely sensed snow observations into a macroscale hydrology model, Adv. Water Resour., 29, 872–886, 2006.

Baba, M. W., Gascoin, S., and Hanich, L.: Assimilation of Sentinel-2 Data into a Snowpack Model in the High Atlas of Morocco, Remote Sens., 10, 1982, https://doi.org/10.3390/rs10121982, 2018.

Bellaire, S. and Jamieson, B.: Forecasting the formation of critical snow layers using a coupled snow cover and weather model, Cold. Reg. Sci. Technol., 94, 37–44, 2013.

Bühler, Y., Marty, M., Egli, L., Veitinger, J., Jonas, T., Thee, P., and Ginzler, C.: Snow depth mapping in high-alpine catchments using digital photogrammetry, The Cryosphere, 9, 229–243, https://doi.org/10.5194/tc-9-229-2015, 2015.

Cluzet, B., Revuelto, J., Lafaysse, M., Tuzet, F., Cosme, E., Pocard, G., Arnaud, L., and Dumont, M.: Towards the assimilation of satellite reflectance into semi-distributed ensemble snowpack simulations, Cold Reg. Sci. Technol., 170, 102918, https://doi.org/10.1016/j.coldregions.2019.102918, 2020.

Doms, G., Förstner, J., Heise, E., Herzog, H. J., Mironov, D., Raschendorfer, M., Reinhardt, T., Ritter, B., Schrodin, R., Schulz, J. P., and Vogel, G.: A Description of the Nonhydrostatic Regional COSMO Model, Part II: Physical Parameterization, LM F90 4.20 38, Consortium for Small-Scale Modelling, Printed at Deutscher Wetterdienst, 63004 Offenbach, Germany, 2011.

Douville, H., Royer, J.-F., and Mahfouf, J.-F.: A new snow parameterization for the Météo-France climate model Part II: validation in a 3-D GCM experiment, Clim. Dynam., 1, 37–52, 1995.

Drusch, M., Del Bello, U., Carlier, S., Colin, O., Fernandez, V., Gascon, F., Hoersch, B., Isola, C., Laberinti, P., Martimort, P., et al.: Sentinel-2: ESA's optical high-resolution mission for GMES operational services, Remote Sens. Environ., 120, 25–36, 2012.

Egli, L. and Jonas, T.: Hysteretic dynamics of seasonal snow depth distribution in the Swiss Alps, Geophys. Res. Lett., 36, L02501, https://doi.org/10.1029/2008GL035545, 2009.

Essery, R.: Large-scale simulations of snow albedo masking by forests, Geophys. Res. Lett., 40, 5521–5525, https://doi.org/10.1002/grl.51008, 2013.

Essery, R. and Pomeroy, J.: Implications of spatial distributions of snow mass and melt rate for snow-cover depletion: theoretical considerations, Ann. Glaciol., 38, 261–265, https://doi.org/10.3189/172756404781815275, 2004.

Essery, R., Morin, S., Lejeune, Y., and Ménard, C. B.: A comparison of 1701 snow models using observations from an alpine site, Adv. Water Resour., 55, 131–148, 2013.

Farinotti, D., Magnusson, J., Huss, M., and Bauder, A.: Snow accumulation distribution inferred from time-lapse photography and simple modelling, Hydrol. Process., 24, 2087–2097, https://doi.org/10.1002/hyp.7629, 2010.

Fedorov, R., Camerada, A., Fraternali, P., and Tagliasacchi, M.: Estimating Snow Cover From Publicly Available Images, IEEE T. Multimedia, 18, 1187–1200, https://doi.org/10.1109/TMM.2016.2535356, 2016.

Gascoin, S., Hagolle, O., Huc, M., Jarlan, L., Dejoux, J.-F., Szczypta, C., Marti, R., and Sánchez, R.: A snow cover climatology for the Pyrenees from MODIS snow products, Hydrol. Earth Syst. Sci., 19, 2337–2351, https://doi.org/10.5194/hess-19-2337-2015, 2015.

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

Griessinger, N., Seibert, J., Magnusson, J., and Jonas, T.: Assessing the benefit of snow data assimilation for runoff modeling in Alpine catchments, Hydrol. Earth Syst. Sci., 20, 3895–3905, https://doi.org/10.5194/hess-20-3895-2016, 2016.

Griessinger, N., Schirmer, M., Helbig, N., Winstral, A., Michel, A., and Jonas, T.: Implications of observation-enhanced energy-balance snowmelt simulations for runoff modeling of Alpine catchments, Adv. Water Resour., 133, 103410, https://doi.org/10.1016/j.advwatres.2019.103410, 2019.

Hall, D. K., Riggs, G. A., and Salomonson, V. V.: Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer data, Remote Sens. Environ., 54, 127–140, https://doi.org/10.1016/0034-4257(95)00137-P, 1995.

Härer, S., Bernhardt, M., and Schulz, K.: PRACTISE – Photo Rectification And ClassificaTion SoftwarE (V.2.1), Geosci. Model Dev., 9, 307–321, https://doi.org/10.5194/gmd-9-307-2016, 2016.

Helbig, N., van Herwijnen, A., and Jonas, T.: Forecasting wet-snow avalanche probability in mountainous terrain, Cold Reg. Sci. Technol., 120, 219–226, https://doi.org/10.1016/j.coldregions.2015.07.001, 2015a.

Helbig, N., van Herwijnen, A., Magnusson, J., and Jonas, T.: Fractional snow-covered area parameterization over complex topography, Hydrol. Earth Syst. Sci., 19, 1339–1351, https://doi.org/10.5194/hess-19-1339-2015, 2015b.

Helbig, N., Bühler, Y., Eberhard, L., Deschamps-Berger, C., Gascoin, S., Dumont, M., Revuelto, J., Deems, J. S., and Jonas, T.: Fractional snow-covered area: scale-independent peak of winter parameterization, The Cryosphere, 15, 615–632, https://doi.org/10.5194/tc-15-615-2021, 2021a.

Helbig, N., Schirmer, M., and Magnusson, J.: Seasonal fractional snow-covered area algorithm, https://doi.org/10.16904/envidat.244, EnviDat [data set], https://doi.org/10.16904/envidat.244, 2021b.

Horton, S. and Jamieson, B.: Modelling hazardous surface hoar layers across western Canada with a coupled weather and snow cover model, Cold. Reg. Sci. Technol., 128, 22–31, 2016.

Huang, C., Newman, A. J., Clark, M. P., Wood, A. W., and Zheng, X.: Evaluation of snow data assimilation using the ensemble Kalman filter for seasonal streamflow prediction in the
western United States, Hydrol. Earth Syst. Sci., 21, 635–650,
https://doi.org/10.5194/hess-21-635-2017, 2017.

Lawrence, D., Fisher, R., Koven, C., Oleson, K., Swenson, S., and Vertenstein, M.: Technical Description of version 5.0 of the Community Land Model (CLM), available at: https://www.cesm.ucar.edu/models/cesm2/land/CLM50_Tech_Note.pdf, last access: 31 March 2018.

López-Moreno, J. I., Revuelto, J., Alonso-Gonzáles, E., Sanmiguel-Vallelado, A., Fassnacht, S. R., Deems, J., and Morán-Tejeda, E.: Using very long-range Terrestrial Laser Scanning to Analyze the Temporal Consistency of the Snowpack Distribution in a High Mountain Environment, J. Mt. Sci., 14, 823–842, 2017.

Luce, C. H., Tarboton, D. G., and Cooley, K. R.: Sub-grid param-
eterization of snow distribution for an energy and mass balance
snow cover model, Hydrol. Process., 13, 1921–1933, 1999.

Magand, C., Ducharne, A., Moine, N. L., and Gascoin, S.: Introduc-
ing Hysteresis in Snow Depletion Curves to Improve the Water
Budget of a Land Surface Model in an Alpine Catchment, J. Hy-
drometeorol., 15, 631–649, https://doi.org/10.1175/JHM-D-13-
091.1, 2014.

Magnusson, J., Gustafsson, D., Hüsler, F., and Jonas, T.: Assimila-
tion of point SWE data into a distributed snow cover model com-
paring two contrasting methods, Water Resour. Res., 50, 7816–
7835, 2014.

Magnusson, J., Wever, N., Essery, R., Helbig, N., Winstral, A., and Jonas, T.: Evaluating snow models with varying process representations for hydrological applications, Water Resour. Res., 51, 2707–2723, https://doi.org/10.1002/2014WR016498, 2015.

Marty, M., Bühler, Y., and Ginzler, C.: Snow Depth Mapping,
https://doi.org/10.16904/envidat.62, 2019.

Masson, T., Dumont, M., Mura, M., Sirguey, P., Gas-
coin, S., Dedieu, J.-P., and Chanussot, J.: An Assess-
ment of Existing Methodologies to Retrieve Snow Cover
Fraction from MODIS Data, Remote Sensing, 10, 619,
https://doi.org/10.3390/rs10040619, 2018.

Mazzotti, G., Currier, W. R., Deems, J. S., Pflug, J. M., Lundquist, J. D., and Jonas, T.: Revisiting Snow Cover Variability and Canopy Structure Within Forest Stands: Insights From Airborne Lidar Data, Water Resour. Res., 55, 6198–6216, 2019.

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

Nagler, T., Rott, H., Malcher, P., and Müller, F.: Assimilation of meteorological and remote sensing data for snowmelt runoffforecasting, Remote Sens. Environ., 112, 1408–1420, 2008.

Nitta, T., Yoshimura, K., Takata, K., O'ishi, R., Sueyoshi, T., Kanae, S., Oki, T., Abe-Ouchi, A., and Liston, G. E.: Representing Variability in Subgrid Snow Cover and Snow Depth in a Global Land Model: Offline Validation, J. Climate, 27, 3318–3330, https://doi.org/10.1175/JCLI-D-13-00310.1, 2014.

Niu, G. Y. and Yang, Z. L.: An observation-based formula-
tion of snow cover fraction and its evaluation over large
North American river basins, J. Geophys. Res., 112, D21101,
https://doi.org/10.1029/2007JD008674, 2007.

Painter, T., Berisford, D., Boardman, J., Bormann, K., Deems, J., Gehrke, F., Hedrick, A., Joyce, M., Laidlaw, R., Marks, D., Mattmann, C., Mcgurk, B., Ramirez, P., Richardson, M.,

Skiles, S. M., Seidel, F., and Winstral, A.: 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, https://doi.org/10.1016/j.rse.2016.06.018, 2016.

Painter, T. H., Rittger, K., McKenzie, C., Slaughter, P., Davis, R. E., and Dozier, J.: Retrieval of subpixel snow covered area, grain size, and albedo from MODIS, Remote Sens. Environ., 113, 868–879, https://doi.org/10.1016/j.rse.2009.01.001, 2009.

Parajka, J. and Blöschl, G.: Validation of MODIS snow cover images over Austria, Hydrol. Earth Syst. Sci., 10, 679–689, https://doi.org/10.5194/hess-10-679-2006, 2006.

Portenier, C., Hüsler, F., Härer, S., and Wunderle, S.: To-
wards a webcam-based snow cover monitoring network:
methodology and evaluation, The Cryosphere, 14, 1409–1423,
https://doi.org/10.5194/tc-14-1409-2020, 2020.

Revuelto, J., López-Moreno, J. I., Azorin-Molina, C., and Vicente-Serrano, S. M.: Topographic control of snowpack distribution in a small catchment in the central Spanish Pyrenees: intra- and inter-annual persistence, The Cryosphere, 8, 1989–2006, https://doi.org/10.5194/tc-8-1989-2014, 2014.

Riboust, P., Thirel, G., Le Moine, N., and Ribstein, P.: Revisiting a simple degree-day model for integrating satellite data: implementation of SWE-SCA hystereses, J. Hydrol. Hydromech., 67, 70–81, 2019.

Roesch, A., Wild, M., Gilgen, H., and Ohmura, A.: A new snow cover fraction parameterization for the ECHAM4 GCM, Clim. Dynam., 17, 933–946, 2001.

Salvatori, R., Plini, P., Giusto, M., Valt, M., Salzano, R., Montag-
noli, M., Cagnati, A., Crepaz, G., and Sigismondi, D.: Snow
cover monitoring with images from digital camera systems, Ital.
J. Remote. Sens., 43, 2, https://doi.org/10.5721/ltJRS201143211,
2011.

Schirmer, M. and Lehning, M.: Persistence in intra-annual
snow depth distribution: 2. Fractal analysis of snow
depth development, Water Resour. Res., 47, W09517,
https://doi.org/10.1029/2010WR009429, 2011.

Schirmer, M., Wirz, V., Clifton, A., and Lehning, M.: Persis-
tence in intra-annual snow depth distribution: 1. Measurements
and topographic control, Water Resour. Res., 47, W09516,
https://doi.org/10.1029/2010WR009426, 2011.

Skaugen, T. and Melvold, K.: Modeling the snow depth vari-
ability with a high-resolution lidar data set and nonlin-
ear terrain dependency, Water Resour. Res., 55, 9689–9704,
https://doi.org/10.1029/2019WR025030, 2019.

Su, H., Yang, Z. L., Niu, G. Y., and Dickinson, R. E.: Enhancing the estimation of continental-scale snow water equivalent by assimilating MODIS snow cover with the ensemble Kalman filter, J. Geophys. Res., 113, D08120, https://doi.org/10.1029/2007JD009232, 2008.

Swenson, S. C. and Lawrence, D. M.: A new fractional snow-
covered area parameterization for the Community Land Model
and its effect on the surface energy balance, J. Geophys. Res.-
Atmos., 117, D21107, https://doi.org/10.1029/2012JD018178,
2012.

Thirel, G., Salamon, P., Burek, P., and Kalas, M.: Assimilation of MODIS snow cover area data in a distributed hydrological model using the particle filter, Remote Sensing, 5, 5825–5850, 2013.
van Herwijnen, A. and Schweizer, J.: Seismic sensor array for monitoring an avalanche start zone: design, deployment and preliminary results, J. Glaciol., 57, 257–264, 2011.

van Herwijnen, A., Berthod, N., Simenhois, R., and Mitterer, C.: Using time-lapse photography in avalanche research, in: Proceedings of the International Snow Science Workshop, Grenoble, France, 950–954, 2013.

Vionnet, V., Martin, E., Masson, V., Guyomarc'h, G., Naaim-Bouvet, F., Prokop, A., Durand, Y., and Lac, C.: Simulation of wind-induced snow transport and sublimation in alpine terrain using a fully coupled snowpack/atmosphere model, The Cryosphere, 8, 395–415, https://doi.org/10.5194/tc-8-395-2014, 2014.

Yang, Z. L., Dickinson, R. E., Robock, A., and Vinnikov, K. Y.: On validation of the snow sub-model of the biosphere atmosphere transfer scheme with Russian snow cover and meteorological observational data, J. Climate, 10, 353–373, 1997.

Zaitchik, B. F. and Rodell, M.: Forward-Looking Assimilation of MODIS-Derived Snow-Covered Area into a Land Surface Model, J. Hydrometeorol., 10, 130–148, https://doi.org/10.1175/2008JHM1042.1, 2009.