# Processing of VENµS Images of High Mountains: A Case Study for Cryospheric and Hydro-Climatic Applications in the Everest Region (Nepal)

Zoé Bessin, Jean-Pierre Dedieu, Yves Arnaud, Patrick Wagnon, Fanny Brun, Michel Esteves, Baker Perry, Tom Matthews

► To cite this version:

Zoé Bessin, Jean-Pierre Dedieu, Yves Arnaud, Patrick Wagnon, Fanny Brun, et al. Processing of VENµS Images of High Mountains: A Case Study for Cryospheric and Hydro-Climatic Applications in the Everest Region (Nepal). Remote Sensing, 2022, 24 (5), pp.1098. 10.3390/rs14051098. hal-03590789

HAL Id: hal-03590789

https://hal.science/hal-03590789

Submitted on 28 Feb 2022

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>

Article

# Processing of VENµS Images of High Mountains: A Case Study for Cryospheric and Hydro-Climatic Applications in the Everest Region (Nepal)

Zoé Bessin <fnref n="1" />,<fnref n="2" />,<fnref n="3" />,† ID, Jean-Pierre Dedieu <fnref n="3" />,*,†, Yves Arnaud <fnref n="3" /> ID, Patrick Wagnon <fnref n="3" />, Fanny Brun <fnref n="3" /> ID, Michel Esteves <fnref n="3" />, Baker Perry <fnref n="4" /> ID and Tom Matthews <fnref n="5" />

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

check for
updates

**Citation:** Bessin, Z.; Dedieu, J.-P.; Arnaud, Y.; Wagnon, P.; Brun, F.; Esteves, M.; Perry, B.; Matthews, T. Processing of VENµS Images of High Mountains: A Case Study for Cryospheric and Hydro-Climatic Applications in the Everest Region (Nepal). *Remote Sens.* **2022**, *14*, 1098. https://doi.org/10.3390/rs14051098

Academic Editors: Arnon Karnieli, Gerard Dedieu and Olivier Hagolle

Received: 23 December 2021

Accepted: 19 February 2022

Published: 23 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Copyright:** © 2022 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In the Central Himalayas, glaciers and snowmelt play an important hydrological role, as they ensure the availability of surface water outside the monsoon period. To compensate for the lack of field measurements in glaciology and hydrology, high temporal and spatial resolution optical remotely sensed data are necessary. The French–Israeli VENµS Earth observation mission has been able to complement field measurements since 2017. The aim of this paper is to evaluate the performance of different reflectance products over the Everest region for constraining the energy balance of glaciers and for cloud and snow cover mapping applied to hydrology. Firstly, the results indicate that a complete radiometric correction of slope effects such as the Gamma one (direct and diffuse illumination) provides better temporal and statistical metrics ($R^2$ = 0.73 and RMSE = 0.11) versus ground albedo datasets than a single cosine correction, even processed under a fine-resolution digital elevation model (DEM). Secondly, a mixed spectral-textural approach on the VENµS images strongly improves the cloud mapping by 15% compared with a spectral mask thresholding process. These findings will improve the accuracy of snow cover mapping over the watershed areas downstream of the Everest region.

**Keywords:** optical remote sensing; VENµS; glaciers; snow; mountains; cloud mapping

# 1. Introduction

## 1.1. Context

Hosting the highest mountains in the world and the largest glacierized area outside of the polar regions, the Himalayas are a unique hydro-climatic system. The numerous Himalayan watersheds supply more than 800 million people downstream, providing water resources from 760,000 km<fnref n="2" /> of snow-covered areas and 40,800 km<fnref n="2" /> of glaciers (including Karakoram and the Himalayas). Therefore these high mountains are often called “water towers” for the buffering role they play in the Asian water cycle [1–3] with snow and ice melt maintaining streamflow during periods of meteorological drought [3–5]. However, this region is a challenge for high-elevation glacio-hydrologic research. Distributed hydro-meteorological information is limited by the difficulty of collecting in situ data [6], translating into limited knowledge of the hydrological regime and its evolution. Gridded

<footnote n="1">¹ Geo-Ocean Laboratory (LGO), University Brest, CNRS, Ifremer, UMR 6538, 29280 Plouzané, France;
zoe.bessin@univ-brest.fr</footnote>

<footnote n="2">² LETG-Brest, University Brest, CNRS, UMR 6554, 29280 Plouzané, France</footnote>

<footnote n="3">³ Institute for Geosciences and Environmental Research (IGE), University Grenoble-Alpes/CNRS/IRD/Grenoble-INP, 38058 Grenoble, France; yves.arnaud@ird.fr (Y.A.); patrick.wagnon@univ-grenoble-alpes.fr (P.W.); fanny.brun@univ-grenoble-alpes.fr (F.B.); michel.esteves@ird.fr (M.E.)</footnote>

<footnote n="4">⁴ Department of Geography and Planning, Appalachian State University, Boone, NC 28608, USA;
perrylb@appstate.edu</footnote>

<footnote n="5">⁵ Department of Geography, King's College London, London WC2B 4BG, UK; tom.matthews@kcl.ac.uk</footnote>

<footnote>* Correspondence: jean-pierre.dedieu@univ-grenoble-alpes.fr; Tel.: +33-456-520-977</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
meteorological data are available in mountain areas (e.g., ERA5 [7] or HAR [8] climate reanalysis) but they remain too coarse to directly resolve processes in mountain catchments. Their performance on a larger scale in mountainous regions is also often unknown, again due to the lack of sufficiently dense observational networks. These limitations create considerable uncertainty about the future evolution of the Himalayan Water Towers [9,10].

The increasing resolution of spaceborne sensors and the lengthening of available time series means that optical remote sensing is offering new opportunities to fill the data gap, by up-scaling in situ measurements at the basin scale, and/or by monitoring unsurveyed regions. This is especially true in the field of snow monitoring [11]. For example, in the Everest region, MODIS (Moderate-Resolution Imaging Spectroradiometer) imagery has been used to help constrain the contributions of snow and ice to streamflow [12–14] and to track seasonal migration of the snow line [15]. A new technical step was reached in terms of spatio-temporal resolution with the Copernicus Sentinel platforms from ESA (European Space Agency), and the recent VENμS (Vegetation and Environment New Micro Satellite) program.

## 1.2. Advantages of VENμS Reflectance Products

The VENμS micro satellite mission is jointly developed and managed by the French Space Agency (CNES) and the Israeli Space Agency (ISA). Since August 2017, VENμS sun-synchronous sensor has provided optical images every 2 days from 12 narrow spectral bands from blue visible (0.424 µm) to near-infrared (0.910 µm), with a 5 m ground resolution. For our study area, the registration time is 11:00 local time. The VENμS program aims to highlight the benefits of repetitive high-resolution image acquisitions for monitoring land surface dynamics. With its unique high spatio-temporal resolution characteristics, the VENμS sensor is well suited to snow monitoring, particularly in high mountain areas that are difficult to access. A total of 130 experimental sites worldwide were selected for the VM1 (VENμS Mission 1) step program (2017–2020), including Mount Everest and the Khumbu valley in Nepal for this study (KHUMBU tile, “Pointing” project). L1C (Level 1C, single date ortho-rectified image expressed in top-of-atmosphere reflectance) and L2A (Level 2A, single date ortho-rectified image expressed in surface reflectance) images are processed by CNES (MAJA chain [16]) providing surface reflectance (SRE) and flat reflectance (FRE) output products, distributed by the THEIA site [17]. The 90 m SRTM (Shuttle Radar Topography Mission) DEM is used by CNES to compute the FRE product with topographic and radiometric corrections using the Gamma process [18,19]. In this way, CNES offers a very high-quality reflectance product with a high revisit frequency and spatial resolution. However, the 8 m HMA (High Mountain Asia) DEM is available in the Khumbu region and could be used to apply the slope effect correction in place of the 90 m SRTM DEM. Thus, we could analyze the influence of the DEM resolution on the topographic correction.

## 1.3. Snow and Cloud Discrimination in Optical Remote Sensing

Concerning the question of snow cover changes at a regional scale adapted to mountainous watersheds, long-term databases offer regular snow map collection (albedo, Snow Cover Area (SCA)); i.e., the National Snow and Ice Data Center (NSIDC) [20] or the Copernicus Snow and Ice Monitoring Service [21]. For this purpose, different methods have been developed to compute changes in terms of the snow cover duration (SCD in days) and snow cover area (in km²). A complete review is given by [22]. The most efficient method to separate clouds and snow is given by the NDSI (Normalized-Difference Snow Index) [23] based on the normalized ratio between Green and SWIR (Short-Wave Infrared) spectral bands:

$$NDSI = \frac{R_{green} - R_{swir}}{R_{green} + R_{swir}} \qquad (1)$$

where $R_{green}$ is the surface reflectance in the green spectral band and $R_{swir}$ is the surface reflectance in the SWIR (1.6 µm) spectral band. Indeed, snow has a very low reflectance in the SWIR band [24]. Mainly based on the NDSI method, MODIS datasets at regional scale
have been widely used [25–29], as well as Sentinel-2 and -3 products [30,31] to compute SCD and SCA. Unlike Sentinel-2 images, the lack of 1.6 µm SWIR spectral band on the VENµS sensor prevents direct snow and ice mapping through the NDSI calculation [28]; thus, alternative methods are needed. Only a few studies have already addressed this issue in other sensors. A method has been developed on Worldview-2 sensor by Bühler et al. [32] using the visible and the highest NIR (Near Infra-Red) band in a normalized-difference band ratios index. It should be noted that no atmospheric or radiometric correction has been applied to these satellite images. In Marchane et al. [33], they apply a supervised classification using the parallelepiped approach with the green, red and NIR bands of the FORMOSAT-2 sensor to retrieve the snow cover. Recently, to overcome this constraint with VENµS, Baba & Gascoin [34] suggested an approach using VENµS database (named the NDSIµ index), based on the Red-edge and NIR spectral bands with classification method for SCA mapping in mountainous environments. However, the 15 FORMOSAT-2 images used in Marchane et al. had almost no clouds, which avoids the difficulty of differentiating snow and cloud. Baba & Gascoin [34] used the VENµS cloud mask to exclude clouds in the process. Therefore, the discrimination between snow and clouds remains unaddressed.

Cloud cover is a major constraint in optical remote sensing, impacting the retrieval of surface properties. Different masking algorithms exist, based on multi-temporal or elevation methods [35,36] but spectral similarity between cloud and snow can be challenging [37,38]. Most solutions developed recently to detect clouds in remote-sensed images involve machine learning or deep learning methods [39–42] using convolutional neural network architectures. Although these methods offer attractive results, the implementation of such a process requires considerable computing resources and is often time consuming. We intend to propose here a simple and parsimonious alternative. Generally, in “land-use” thematic applications of optical remote sensing, images with strong cloud cover are not processed or rejected because they prevent any analysis of the ground surfaces. On the contrary, we decide in this study to use all available VENµS images, including those with considerable cloud cover. Since CNES already provides a robust cloud mask for each VENµS date, the purpose is to enhance it and to discriminate between clouds and snow. This will allow us to study the seasonal and spatio-temporal evolution of clouds versus orography and monsoon/non-monsoon impacts related to local precipitation, and to complement climatological studies in the Himalayas. Moreover, such analysis at high spatial resolutions might also be of interest for surface radiative forcing.

## 1.4. Aims

The aims of this paper are therefore: (i) to evaluate the performances of CNES reflectance products (SRE and FRE) relative to ground-based observation and to explore the possibility of improving them by including a high resolution DEM in the processing chain; (ii) to combine spectral and textural parameters to improve cloud mapping without the help of manual intervention; and (iii) to use these improved cloud masks to provide better SCA maps at a 5 m ground resolution and at a 2-day frequency.

We anticipate that more accurate remote-sensed surface reflectance products will offer benefits across disciplines, including:

* in glaciology: to provide spatially distributed albedo maps, ideally every 2 days, from optical remote-sensing data, as already documented by [12–14], and useful for cryosphere modeling at a regional scale [43,44];

* in climatology: to analyze the occurrence of cloud cover and its spatio-temporal variability at 5 m resolution for radiative impact on surface energy as given by [45,46];

* in hydrology: to provide snow cover area (SCA) maps at 2-day temporal and 5 m spatial resolution, useful to optimize hydrological modeling in the Everest/Khumbu basins [5,12,14].

Firstly, we describe the study area, materials and methods in Section 2. Results are detailed and discussed in Sections 3 and 4 respectively. Finally, the conclusions of this study are presented in Section 5.
the complete application area (1485 km²). From 27 November 2017 to 2 June 2020, out of 460 dates, 238 VENµS images were available and preprocessed by CNES over the Everest/Khumbu area. Furthermore, 104 dates were not acquired for technical reasons and 118 images were discarded by CNES due to high cloud cover extent (over 70%), mainly during the monsoon.

<table><caption>**Table 1.** Summary of VENµS characteristics and CNES products (SRE: surface reflectance; FRE: flat reflectance).</caption><thead><tr><th>Sensor</th><th>Wavelengths (µm)</th><th>Temporal Resolution</th><th>Spatial Resolution</th><th>Tile Size (km²)</th><th>Number of SRE and FRE Processed Images</th></tr></thead><tbody><tr><td>VENµS</td><td>12 bands: 0.420 to 0.910</td><td>2-day</td><td>5-m</td><td>27.5 × 54</td><td>238 × 2</td></tr></tbody></table>

### 3.1.2. DEM Products

Concerning our Nepalese site, the CNES preprocessing chain in the VM1 phase used two different DEM products: (i) the geometric orthorectification using the High Mountain Asia (HMA) DEM at 8 m ground resolution, but (ii) the radiometric correction of slope effects was only applied based on the SRTM DEM at 90 m resolution. In this study, we use the HMA DEM, aiming to evaluate if a higher-resolution DEM could enhance the reflectance product relative to CNES. The HMA DEM was resampled using the bicubic method at 5 m resolution by S. Gascoin (CESBIO-CNES) to be phased with VENµS images.

### 3.1.3. In-Situ Measurements

The two watersheds of Pheriche and Dingboche have been instrumented with a water pressure level sensor since 2010 for runoff monitoring [5]. Data from three automatic weather stations (AWS) are included in our area: Pyramid (5035 m a.s.l.) on the edge of Khumbu glacier, West Changri Nup glacier (5387 m a.s.l.) and the Everest South Col site (7945 m a.s.l.) [6] (Figure 1). Table 2 provides the surface types for each AWS as well as the sensor types of the 3 data used for this study: albedo values are expressed here using Equation (2). The AWSs also provide air temperature (°C), snow height (m) (only available for West Changri Nup site) and precipitation (mm water equivalent (w.e.), only for Pyramid site).

$$ \text{Albedo} = \frac{\text{mean}_{24\text{h}}(S_{\text{wout}})}{\text{mean}_{24\text{h}}(S_{\text{win}})}, \qquad (2) $$

where $S_{wout}$ and $S_{win}$ are respectively reflected and incident short-wave radiation in W.m⁻².

<table><caption>**Table 2.** Specifications (location, elevation, sensor types, surface types) of each AWS used for the VENµS application.</caption><thead><tr><th>Site</th><th>Lat/long</th><th>Elevation (m a.s.l.)</th><th>Albedo</th><th>Air Temp.</th><th>Snow Height</th><th>Precip.</th><th>Surface Type</th></tr></thead><tbody><tr><td>Pyramid</td><td>N 86°48'47"/E 27°57'32"</td><td>5035</td><td>Kipp and Zonen CNR4</td><td>Vaisaila HMP45C</td><td>NA</td><td>Geonor T-200</td><td>grassy surface</td></tr><tr><td>Changri Nup</td><td>N 86°46'40"/E 27°58'57"</td><td>5387</td><td>Kipp and Zonen CNR4</td><td>Vaisaila HMP155</td><td>Campbell SR50A</td><td>NA</td><td>Debris-covered glacier</td></tr><tr><td>South Col</td><td>N 27°58'18"/E 86° 55'46"</td><td>7945</td><td>Hukseflux NR01</td><td>Vaisaila HMP155</td><td>NA</td><td>NA</td><td>rocky bare surface</td></tr></tbody></table>

### 3.1.4. VENµS Cloud Mask "CLM_XS"

The CNES cloud mask (named CLM_XS) is an 8-bit image, and each bit refers to a cloud type and a detection method (Table 3). To detect clouds on VENµS images, the MAJA chain applies three tests using: (i) the 1.3 µm band available on LANDSAT-8 and SENTINEL-2; (ii) a multi-temporal test on the blue surface reflectance increase and (iii) the
correlation of the pixel neighborhood with previous images. For our study, we extracted the 7th bit containing the high clouds detected by stereoscopy (value 128).

<table><caption>**Table 3.** CLM_XS product description.</caption><thead><tr><td>Bit</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>0</td><td>1</td><td>All clouds except the thinnest and all shadows</td></tr><tr><td>1</td><td>2</td><td>All clouds (except the thinnest)</td></tr><tr><td>2</td><td>4</td><td>Cloud shadows cast by a detected cloud</td></tr><tr><td>3</td><td>8</td><td>Cloud shadows cast by a cloud outside image</td></tr><tr><td>4</td><td>16</td><td>Clouds detected via mono-temporal thresholds</td></tr><tr><td>5</td><td>32</td><td>Clouds detected via multi-temporal thresholds</td></tr><tr><td>6</td><td>64</td><td>Thinnest clouds</td></tr><tr><td>7</td><td>128</td><td>High clouds detected by stereoscopy</td></tr></tbody></table>

## 3.2. Methods

### 3.2.1. Slope Effect Correction

Atmospheric and radiometric slope effect corrections over mountainous areas are necessary for state-of-the-art preprocessing in optical remote sensing [53,54]. This post-processing is essential for surface albedo retrieval in rugged terrain, particularly for glaciers and snow-covered areas [55–58]. We analyze the performance of two approaches to correct the differential solar illumination effects: (1) the cosine correction taking into account only the direct illumination (lambertian case) (Figure 2) [59,60]; and (2) the Gamma correction, more complex, considering both direct and diffuse illuminations (non-Lambertian case) [18]. The Gamma process, logically more efficient for snow cover, is applied by CNES to generate the VENµS FRE product (flat reflectance) from the SRE product (surface reflectance), but is based on the coarse 90 m SRTM DEM.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Principles of relative geometry of sun, sensor and target. $\theta_s$ and $\phi_s$ are (respectively) the sun zenith and azimuthal angles (yellow); $\theta_n$ and $\phi_n$ are the slope value and orientation (green), respectively. $\theta_v$ is the sensor zenith and $\phi_v$ is the sensor azimuth (blue). The illumination angle $\gamma$ (magenta) is the relative angle between the normal to the surface (slope) and the solar ray. All these angles are linked by Equations (3) and (4) using Euler's rotation theorem.</figcaption></figure>

Our objective here is to test whether a simple correction of direct illumination (cosine correction), if based on a higher-resolution DEM (8-m HMA DEM in this case), would be more accurate than the Gamma reflectance correction applied by CNES. In order to apply the cosine correction to the VENµS SRE product, we developed a Python script whose major steps are described below and in Figure 3.
<figure><img src="image_6.png" /><figcaption>Figure 3. TopCos HMA products (B11 and B7) processing flow chart. This workflow describes the method used to correct VENµS surface reflectance (SRE) product with the cosine algorithm and the HMA DEM.</figcaption></figure>

As illustrated in Figure 2, the use of a DEM makes it possible to determine the respective illumination cos(γ) of each specific slope applying the following formula:

$$ \cos(\gamma) = \cos(\theta_s) \cos(\theta_n) + \sin(\theta_s) \sin(\theta_n) \cos(\varphi_s - \varphi_n), \quad (3) $$

where θs and φs are the sun zenithal and azimuthal angles; θn and φn are the slope value and orientation, respectively. Equation (3) creates the illumination angle and the normalized horizontal reflectance is then calculated by Equation (4):

$$ R_{hor} = \begin{cases} \frac{R \times \cos(\theta_s)}{\cos(\gamma)}, & \cos(\gamma) > 0.2 \\ \text{No Date Value}, & \cos(\gamma) \le 0.2 \end{cases} \quad (4) $$

where R is the input pixel sloping reflectance and R<sub>hor</sub> is the output pixel corrected as in horizontal condition. Due to the low reflectance values in shaded areas, artifacts are created. We use a threshold fixed at 0.2 as in [53], below which R<sub>hor</sub> data are discarded, as the correlation between the ground reflectance data and the corrected reflectance of the satellite images is too low. The cosine correction can be simply applied to any reflectance band, but some low illumination angles are over-corrected [61]. The reason is that these areas receive more diffuse radiation not accounted for in the cos(γ) direct calculation. We apply Equations (3) and (4) using the HMA DEM. Our cosine output product is then named TopCos HMA (referred to as such here-after). This model is run on the VENµS SRE product for (i) the B11 near-infrared spectral band (865 µm) concerning albedo application, because it is least impacted by atmospheric effects, and (ii) the B7 red spectral band (662 µm) for cloud and snow mapping applications (Figure 3).

Additionally, we compute a fourth reflectance product, named TopCos SRTM (Table 4), which is based on the SRE product corrected with the cosine correction using the SRTM DEM (90 m). The aim is to analyze the influence of different DEM spatial resolutions with the same radiometric correction method.
<table><caption>**Table 4.** Summary of the DEM and corrections used for each reflectance product in this study.</caption><thead><tr><th>Product Name</th><th>DEM</th><th>Correction</th></tr></thead><tbody><tr><td>SRE</td><td>None</td><td>None</td></tr><tr><td>FRE</td><td>SRTM 90-m</td><td>Gamma</td></tr><tr><td>TopCos HMA</td><td>HMA 8-m</td><td>Cosine</td></tr><tr><td>TopCos SRTM</td><td>SRTM 90-m</td><td>Cosine</td></tr></tbody></table>

### 3.2.2. Reflectance Product Statistical Evaluation Strategy

A statistical comparison is conducted among the four reflectance products; thus, the contribution of the bias and STD to the RMSE is calculated for each reflectance product. Appendix A exhibit the detailed metrics results for each site. As a reminder, the root-mean-square error (RMSE) is retrieved by the bias and dispersion (STD) calculations. Their respective contribution (%) is simply given by:

$$ \text{Bias contribution} = \frac{100 \times \text{Bias}^2}{\text{RMSE}^2}, \qquad (5) $$

and,

$$ \text{STD contribution} = \frac{100 \times \text{STD}^2}{\text{RMSE}^2}, \qquad (6) $$

with:

$$ \text{RMSE}^2 = \text{Bias}^2 + \text{STD}^2, \qquad (7) $$

where RMSE is the root-mean-square error and STD is the standard deviation. The Pearson coefficient R² is also calculated with:

$$ R^2 = 1 - \frac{SS_{res}}{SS_{tot}}, \qquad (8) $$

where $SS_{res}$ is the residual sum of squares and $SS_{tot}$ is the total sum of squares.

Finally, in order to highlight the best product and the best site in terms of performance, we use a normalized skill score. The transformation to a generic skill enables us to avoid the problem of having different units for the different metrics. The normalized skill score (NSS) for each kind of metric is calculated with:

$$ NSS_{\text{RMSE}}(\text{site, product}) = 1 - \frac{\text{RMSE (site, product)}}{\max(\text{RMSE})}, \qquad (9) $$

where $NSS_{\text{RMSE}}$ represents the normalized RMSE skill score for a specific site and product, RMSE is the root-mean-square error for a specific site and product, and $max(RMSE)$ is the highest RMSE among all sites and products.

$$ NSS_{\text{Bias}}(\text{site, product}) = 1 - \frac{|\text{Bias (site, product)}|}{\max(|\text{Bias}|)}, \qquad (10) $$

where $NSS_{\text{Bias}}$ represents the normalized bias skill score for a specific site and product, $|\text{Bias}|$ is the bias absolute value for a specific site and product, and $max(Bias)$ is the highest bias among all sites and products.

$$ NSS_{\text{STD}}(\text{site, product}) = 1 - \frac{\text{STD (site, product)}}{\max(\text{STD})}, \qquad (11) $$
where $NSS_{STD}$ represents the normalized standard deviation skill score for a specific site and product, STD is the standard deviation for a specific site and product, and $max(STD)$ is the highest STD among all sites and products.

$$ NSS_{R^2}(\text{site, product}) = \frac{R^2(\text{site, product})}{\max(R^2)}, \qquad (12) $$

where $NSS_{R^2}$ represents the normalized Pearson coefficient skill score for a specific site and product, $R^2$ is the Pearson coefficient for a specific site and product, and $max(R^2)$ is the highest $R^2$ among all sites and products. Once all skill scores are calculated for each value at each site, we calculate the average normalized skill score (ANSS) for each type of metric, by product and by site (ANSS by row and column). Higher values correspond to better performance.

### 3.2.3. Cloud Mask Processing

As a reminder, in this study, our goal is not only to apply a cloud mask to the VENµS images to retrieve a satisfactory snow mapping chain, but also to use all cloudy images in the aim of mapping the cloud cover evolution over the measuring period for a climate application (monsoon regime). Based on the Hagolle algorithm [51], we observe that the CLM_XS product applied to the SRE product obscures several snow-covered areas in our study site, which are clearly visible on the Level-1C images of VENµS. This is partly due to the 2-pixel-wide dilation of the cloud mask at reduced resolution, applied to include thin clouds at the edge of thicker clouds [16]. Due to the lack of a SWIR band on the VENµS sensor, we have developed another approach dedicated to cloud mapping improvement. We aim to leverage and refine the CNES cloud masks to retrieve more information about the snow cover on the ground and create cloud and SCA maps databases with a 2-day frequency, at a 5 m spatial resolution. We therefore propose a complementary process based on a hybrid textural/radiometric approach in order to improve and optimize the CLM_XS product. The Python algorithm to retrieve cloud masks for VENµS is part of the VENUS_coseine.correction algorithm.

Mainly applied during the 1970s for geographical application, textural features extraction from optical data are based on mathematical morphology applied to grey-tone images. The Haralick textural decomposition is applied on co-occurrence matrices [62,63] and offers attractive results using entropy and energy as parameters. Energy indicates the angular second-moment feature (ASM), which is a measurement of local grey-level uniformity. Taking into account the textural behavior of the snow and clouds according to uniformity, the energy parameter is more appropriate for this study than entropy. Energy is defined by the following relationship:

$$ \text{Energy} = \sqrt{\text{ASM}} \qquad (13) $$

where:

$$ \text{ASM} = \sum_{i} \sum_{j} p_{d}^{2}(i, j) \qquad (14) $$

with $p_d(i, j)$ as (i,j)-th element of the grey-level co-occurrence matrix [63]. We also use the NDVI (Normalized-Difference Vegetation Index) between red and near-infrared bands (respectively B7 and B11) to help separate snow, ice and cloud features from other land uses. NDVI is defined as:

$$ \text{NDVI} = \frac{R_{\text{NIR}} - R_{\text{red}}}{R_{\text{NIR}} + R_{\text{red}}} \qquad (15) $$

where $R_{NIR}$ is the surface reflectance in the NIR (861.1 nm) spectral band and $R_{red}$ (666.2 nm) is the surface reflectance in the red spectral band.

Our cloud mask product, named CloudCov, is the intersection of masks based on the NDVI and Energy values (Table 5), and the CLM_XS product. Thresholds were determined by radiometric tests on several VENµS dates in winter, with lots of snow. The CLM_XS
product is a conservative mask, which overestimates the cloud cover. This is why we intersect it (∩ denotes the intersection function) with the combination of NDVI and Energy, which tends to miss some cloud-covered areas. The CloudCov product is defined according to Equation (16):

$$ \text{CloudCov} = \text{NDVI} \geq -0.06 \cap \text{NDVI} \leq 0.05 \cap \text{CLM\_XS} = 128 \cap \text{Energy} > 0.8 \quad (16) $$

<table><caption>**Table 5.** NDVI and Energy thresholding window values.</caption><thead><tr><th rowspan="2"></th><th colspan="2">NDVI</th><th>ENERGY</th></tr><tr><th>Low</th><th>High</th><th></th></tr></thead><tbody><tr><td>Clouds</td><td>-0.06</td><td>0.05</td><td>&gt;0.8</td></tr><tr><td>Snow</td><td>-0.16</td><td>-0.02</td><td>&gt;0.8</td></tr></tbody></table>

Clouds over a snow-covered area that are not detected in CLM_XS are also not detected as clouds in our cloud mask.

Once the CloudCov mask has been computed for each VENµS date using TopCos HMA products (NDVI from B7 and B11), all TopCos HMA (B11 only) products (138 dates) are re-computed, this time considering the CloudCov mask. This second process enables the provision of TopCos HMA reflectance products where cloud pixels have the value 128, and thus can be easily discriminated for further applications (Figure 4 describes the complete workflow).

<figure><img src="image_7.png" /><figcaption>**Figure 4.** CloudCov product processing flow chart. This workflow describes the method used to enhance the CNES cloud mask (CLM_XS) with NDVI and Energy thresholds.</figcaption></figure>

### 3.2.4. Cloud Mask Statistical Evaluation Strategy

In order to assess the concordance between a control sample and the cloud mask product, we use the Kappa coefficient [64,65], equivalent to the Heidke Skill Score (HSS) [27,34,66]. By calculating the improvement in the algorithm's behavior relative to a reference behavior such as chance, the HSS measures relative skill. This is expressed as the proportion of pixels correctly classified relative to the number of pixels correct by chance in the complete lack of skill [66]. It allows us to evaluate the classification performances of the CloudCov mask at the pixel scale. The CloudCov mask, manually corrected and improved, is considered to be the actual situation. Table 6, based on Equations (17) and (18), summarizes the metrics used to assess the CloudCov mask.
<table><caption>**Table 6.** Confusion matrix statistics used to assess CLM and CloudCov products.</caption><thead><tr><th>Recall</th><th>Accuracy</th><th>Precision</th><th>Kappa (HSS)</th></tr></thead><tbody><tr><td><sup>TP</sup>/<sub>TP+FN</sub></td><td><sup>TP+TN</sup>/<sub>Total</sub></td><td><sup>TP</sup>/<sub>TP+FP</sub></td><td><sup>p<sub>0</sub> - p<sub>e</sub></sup>/<sub>1 - p<sub>e</sub></sub></td></tr></tbody></table>

The Kappa coefficient is retrieved using:

$$p_0 = \frac{TP + TN}{TP + FP + FN + TN} \qquad (17)$$

and:

$$p_e = \frac{(TN + FP) \times (TN + FN) + (FN + TP) \times (TP + FP)}{(TP + FP + FN + TN)^4} \qquad (18)$$

where $p_0$ is the observed proportionate agreement, $p_e$ is the overall random agreement probability, $TN$ is the true negative pixels; $TP$ is the true positive pixels; $FN$ is the false negative pixels and $FP$ is the false positive pixels. $TN$ means the algorithm correctly predicted that the pixel has no cloud. $TP$ means the algorithm correctly predicted that the pixel has clouds.

### 3.2.5. SCA Outputs

The SCA product is generated using the TopCos and CloudCov chains. Following a similar reasoning to the CloudCov mask, snow pixels are extracted by calculating the following intersection (Equation (19)). Note that a CloudCov value of 0 means no cloud.

$$\text{Snow pixels} = \text{NDVI} \geq -0.16 \cap \text{NDVI} \leq -0.02 \cap \text{CloudCov} = 0 \cap \text{Energy} > 0.8 \quad (19)$$

If these conditions are fulfilled, a snow pixel will get the value 1. Then we code cloud pixels with the value 128 (using the CloudCov mask computed formerly) and other land uses with a value of 0. Figure 5 describes the SCA product workflow and example of maps are given in Figure A1. Under cloudy conditions, our classification algorithm clearly distinguishes snow from clouds. Note that bare ice areas of glaciers are not classified as “snow”, while accumulation areas are. Snow-covered lakes are included, and if not snow-covered, are classified as 0 (other land use).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** SCA product processing flow chart. This workflow presents the method used to retrieve snow cover area maps using the CloudCov mask and thresholds on the NDVI and Energy.</figcaption></figure>
# 4. Results

## 4.1. Producing Reflectance Maps from VENµS Images as a Proxy of Albedo Maps

### 4.1.1. Comparison of TopCos Products (HMA, SRTM) and CNES Products (SRE, FRE)

To compare those four products, we select a subset (Figure 1) in the NE part of the VENµS tile, representative of fully snow-covered areas without cloud coverage (Figure A2). Figure 6 presents the frequency of reflectance occurrences for SRE, FRE and the two TopCos products, HMA and SRTM. Note that for Figure 6, normalization of all data was performed using Equation (20):

$$z_i = \frac{x_i - \min(x)}{\max(x) - \min(x)}, \qquad (20)$$

where $z_i$ is the ith normalized data and $x = (x_1, ..., x_n)$ is the reflectance for each pixel.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Histograms for SRE, FRE and TopCos (HMA and SRTM) normalized products for a subset of the 2018-03-27 VENµS image.</figcaption></figure>

The plots in Figure 6 are representative of the general pattern observed in the complete database. On the one hand, for SRE and FRE products, over the entire range of values, the frequency of occurrence of each reflectance value is homogeneous. On the other hand, TopCos products (HMA and SRTM) indicate reflectance values with a peak around 0.25 (Figure 6). As a reminder, when charting the normalized reflectance occurrence frequency, peak values are left-shifted, to fit the 0–1 interval. Low reflectance values with limited dispersion mean that the image has limited radiometric contrast. The TopCos products have less contrasts than the SRE and FRE. In TopCos, there is indeed a reduction in the difference between illuminated and shaded areas induced by the cosine correction [67].

Figure 7 shows the probability density of the reflectances for the four reflectance products on the same subset as Figure 6. All products have their mode close to 1, characteristic of the snow response. On the one hand, SRE and FRE reflectances are constrained to 1.3 due to a LUTs (Look Up Tables) limit in the atmospheric correction process of VENµS images. On the other hand, the cosine method is known for its overcorrection effect, especially on less illuminated regions: the lower the value of $cos(γ)$, the higher the overcorrection [53]. This phenomenon may be accentuated by the darker pixels in neighboring shadow areas (no correction applies since $cos(γ)$ is negative or zero) [59]. Here, the cosine method is applied only for $cos(γ) > 0.2$ ($γ < 78°$) to avoid artifacts. The larger reflectance amplitude observed for TopCos products is logically induced by the cosine correction, whereas FRE computed via the Gamma method (direct and diffuse illumination) is closer to broad band albedo from ground data. It should be noted that the fine HMA DEM (5 m) reduces this overcorrection in regard to the larger SRTM one (90 m), suggesting that a finer DEM enables better consideration of local topography, and thus more realistic reflectances.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Statistical probability for SRE, FRE and TopCos products (HMA and SRTM) for a subset of the 27 March 2018 VENμS image.</figcaption></figure>

### 4.1.2. Difference between the Four Reflectance Products

In order to understand the high-resolution DEM influence on the radiometric correction, we analyzed the differences between the various reflectance products (Figure A3). The analysis is carried out here on the 7 December 2017 VENμS image, with no cloud and average snow cover extent. As a reminder, a summary of all the products computed is available in Table 4. Comparison between the two CNES products (SRE and FRE) shows that surfaces with steep slopes exposed to the sun are corrected negatively, with the reflectance of these areas decreasing when the gamma correction is applied. In contrast, flat areas less exposed to the sun are corrected positively (i.e., reflectance increases in the FRE product). The difference between the SRE and TopCos HMA indicates that the DEM resolution increases the level of detail and thus the TopCos HMA product should be preferred. Comparison between the two corrected products FRE and TopCos HMA reveals that the TopCos HMA product further reduces the high reflectance values in sun-exposed areas.

Finally, to analyze only the influence of DEM resolution, we compare the two TopCos products (HMA and SRTM); the difference between them shows the same mean bias correction all over the tile (around 0.01 difference in their mean values). Use of an initial DEM with an 8 m resolution, close to the VENμS spatial resolution (5-m) offers a sharp relief description and more variation in the reflectance within a neighborhood which is a benefit in snow detection.

### 4.1.3. Comparison of TopCos and CNES Products with Shortwave Albedo Measurements

* Temporal evolution

Here we compare satellite data registered in a spectral band of near-infrared (B11, 865 µm), with a narrow bandwidth of 40 µm, with broadband albedo sensors on the ground (0.305 < λ < 2.8 µm) for incident and reflected short-wave radiation. Although not a like-for-like comparison, we pursue this comparison to explore the possibility of using the remotely sensed data to infer broadband albedo relevant for the surface energy balance. Figure 8 illustrates the results for the three ground measurement sites. It is important to note that during the monsoon (JJAS), Pyramid site is not covered with snow due to its lower elevation (5 035 m a.s.l.), leading to low reflectance values. Apart from the monsoon, the four reflectance products follow the ground albedo pattern with each instance of snowfall characterized by a sharp albedo increase, followed by a decay. At Pyramid (Figure 8a), when there is snow on the ground, TopCos HMA and FRE products are both closer to the terrain albedo values, whereas the SRE and TopCos SRTM products are respectively impacted by (i) the lack of slope radiometric correction, and (ii) the low DEM resolution effect. Changri Nup site (Figure 8b) is regularly covered by snow in winter (DJF), and
unlike Pyramid site, is consistent with more frequent high reflectance values. During the monsoon, there is no snow on the ground, the reflectance values are low and similar to the ground albedo values as for the Pyramid site. When snowfall occurs in the winter period, the four reflectance images follow the temporal variability of the in situ albedo values; the reflectance increases, then slowly decreases (gradual melting of the snow), and the four products offer a similar response considering alternatively each pair of products as for Pyramid site. At South Col site (Figure 8c) we can observe several high reflectance values during the monsoon, suggesting significant snow cover, induced by solid precipitation at such elevation (7945 m a.s.l.).

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Temporal comparison between four reflectance products and terrain albedo measurements; (a) Pyramid site (May 2019 to June 2020); (b) Changri Nup site (Dec 2017 to Nov 2019); (c) South Col site (May 2019 to June 2020).</figcaption></figure>
The three sites of this study have distinct seasonal behavior due to their specific locations and surface states (Table 2), but similar temporal variabilities occur between the four reflectance products and ground measurements: (i) when there is little or no snow on the ground, all four reflectance products exhibit low reflectances close to the in situ albedo measurements. The three corrected products (FRE, TopCos HMA and TopCos SRTM) are close to each other; (ii) during snowfalls the four products follow the ground albedo rise. More precisely, TopCos products (HMA and SRTM) are closer to the terrain albedo when the site is snow-covered while CNES products (SRE and FRE) tend to underestimate the in situ albedo values. The raw SRE products are the furthest from the observations.

* Statistical comparison

Table 7 offers a synthesis of all metrics for the four reflectance products compared to albedo measured in situ. The details of the normalized skill scores used to calculate the ANSS are described in Table A1. We can observe that, globally, all metrics offer satisfactory results. In particular, we can see that (i) for the FRE product, RMSE values are low on the three sites, and (ii) from the three corrected products (FRE, TopCos HMA and TopCos SRTM), the FRE product has the best mean RMSE on the three sites. The dispersion is always the most impactful on the RMSE. Using the average normalized skill score (ANSS) reveals that the best product according to statistical performances is the FRE from CNES, and the site for which this is most effective is Pyramid (Figure A4). However, since Pyramid site is rarely snow-covered, observing the correlation between reflectance products and in situ albedo is more difficult than the two other sites. Changri Nup (Figure A5) is the site with the lowest performance but it has more temporal variability. On the one hand, the CNES FRE and SRE products are similar and FRE introduces great metrics (STD contributing to the RMSE for 84% while the bias contributes only for 11%). On the other hand, the TopCos HMA has a dispersion and bias that similarly impact the RMSE (60% and 40% respectively) and TopCos products present significant biases (around 0.1). South Col site (Figure A6) is the second-best site of application according to Table 7. The two CNES products are similar, likewise at the Changri Nup site. The FRE product offers great metrics; however, its RMSE is larger than TopCos products (dispersion and bias contribute equally to the RMSE of the FRE product, respectively 51% and 41%). TopCos HMA product shows a low RMSE and the lowest bias, so it would therefore be less affected by overcorrection effects (dispersion and bias contributions to RMSE are 100% and 2% respectively). We can observe that Pyramid and South Col sites have limited temporal variability, with respectively little or much snow cover throughout the year, which impacts the statistical results.

<table><caption>**Table 7.** Average normalized skill score derived from metrics comparison between SRE, FRE, TopCos HMA, and TopCos SRTM products and in situ albedo AWS data for the three studied sites.</caption><thead><tr><th>Site</th><th>Metrics</th><th>SRE</th><th>FRE</th><th>TopCos HMA</th><th>TopCos SRTM</th><th>Average Normalized<br>Skill Score (ANSS) by<br>Site</th></tr></thead><tbody><tr><td rowspan="4">Pyramid</td><td>RMSE</td><td>0.10</td><td>0.08</td><td>0.08</td><td>0.11</td><td>0.58</td></tr><tr><td>Bias</td><td>-0.05</td><td>-0.03</td><td>-0.03</td><td>-0.05</td><td>0.71</td></tr><tr><td>STD</td><td>0.09</td><td>0.08</td><td>0.08</td><td>0.09</td><td>0.50</td></tr><tr><td>R<sup>2</sup></td><td>0.72</td><td>0.71</td><td>0.71</td><td>0.72</td><td>0.94</td></tr><tr><td rowspan="4">Changri Nup</td><td>RMSE</td><td>0.12</td><td>0.12</td><td>0.22</td><td>0.18</td><td>0.27</td></tr><tr><td>Bias</td><td>-0.04</td><td>-0.04</td><td>-0.14</td><td>-0.10</td><td>0.43</td></tr><tr><td>STD</td><td>0.12</td><td>0.11</td><td>0.17</td><td>0.14</td><td>0.21</td></tr><tr><td>R<sup>2</sup></td><td>0.76</td><td>0.76</td><td>0.73</td><td>0.74</td><td>0.98</td></tr><tr><td rowspan="4">South Col</td><td>RMSE</td><td>0.14</td><td>0.14</td><td>0.13</td><td>0.12</td><td>0.40</td></tr><tr><td>Bias</td><td>0.08</td><td>0.09</td><td>0.02</td><td>0.04</td><td>0.59</td></tr><tr><td>STD</td><td>0.11</td><td>0.10</td><td>0.13</td><td>0.11</td><td>0.34</td></tr><tr><td>R<sup>2</sup></td><td>0.71</td><td>0.73</td><td>0.71</td><td>0.73</td><td>0.95</td></tr><tr><td rowspan="4">Average Normalized<br>Skill Score (ANSS) by<br>product</td><td>ANSS RMSE</td><td>0.45</td><td>0.48</td><td>0.35</td><td>0.38</td><td></td></tr><tr><td>ANSS Bias</td><td>0.60</td><td>0.62</td><td>0.55</td><td>0.55</td><td></td></tr><tr><td>ANSS STD</td><td>0.37</td><td>0.43</td><td>0.25</td><td>0.33</td><td></td></tr><tr><td>ANSS R<sup>2</sup></td><td>0.96</td><td>0.96</td><td>0.94</td><td>0.96</td><td></td></tr></tbody></table>
#### 4.1.4. Distribution of Reflectance Values according to the Correction Used

By presenting the reflectance values under classes, an analysis of the reflectance distribution according to each reflectance product, and thus according to the correction method, can be performed (Figure 9). Pyramid site (low snow cover) (Figure 9a) reveals a similar behavior for the four reflectance products in low reflectances (<0.4). Above 0.4, SRE and TopCos SRTM products follow the same trend, as does the FRE-TopCos HMA pair. For reflectance higher than 0.6, Figure 9a suggests that the cosine correction algorithm used to compute TopCos HMA has the same effect as the MAJA Gamma algorithm (CNES). However, most reflectance values are below 0.7, which means Pyramid site is rarely covered by snow. Results on the Changri Nup site (more often snow-covered) (Figure 9b) indicate that for reflectances below 0.6, all the studied products are very similar to each other. However, for higher reflectances, two product pairs can be distinguished: the TopCos products show analogous behavior between each other, as do CNES products. At South Col site (Figure 9c) the behavior of all the products remains the same for low reflectances (<0.6). Above this threshold, the same two pairings as in the Changri Nup site are observed.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Boxplot comparison between the four reflectance products; (a) Pyramid site; (b) Changri Nup site; (c) South Col site.</figcaption></figure>
## 4.2. Cloud Mapping

### 4.2.1. Statistical Analysis of the CloudCov Product

Our algorithm developed to retrieve cloud cover is presented in Section 3.2.3. (Table 5 and Figure 4) and an example is presented in Figure 10. The cloud mask produced is named CloudCov. Our approach enables significant optimization of the CNES cloud mask CLM_XS (Figure 10b). A representative example (8 August 2018), describing more precisely the cloud outlines versus snow-covered and non-snow-covered areas on the ground, is presented in Figure 10c. On average, the CLM_XS mask classifies 27% of the pixel as clouds, while CloudCov classifies only 12% as clouds on the same images. This is a significant gain in usable imagery. However, as the CloudCov product extends near the edge of the clouds, it also contains fewer cloud-cast shadows (Figure 10). Depending on the objective, one can still favor the CLM_XS mask as it is more conservative.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Cloud products and masks footprint comparison. (a) VENµS image of 8 August 2018, (b) CNES Cloud mask (CLM_XS), (c) CloudCov mask.</figcaption></figure>

The CloudCov product estimated from our algorithm is then manually corrected (manual delineation) and used as a reference in the confusion matrix (Figure 11), computed on the 238 VENµS dates, to quantify the gain obtained with the CloudCov mask. Our cloud mask contains 96% well-classified pixels, while CLM_XS cloud mask only has 81% (Figure 11 and Table 8). Thus, the CloudCov mask is 15% more accurate than CLM_XS mask, but the precision improvement is even more significant, with a 35% gain compared to the CLM_XS mask.

<table><caption>**Table 8.** Confusion matrix statistics between CLM and CloudCov products.</caption><thead><tr><td></td><td>Recall</td><td>Accuracy</td><td>Precision</td><td>Kappa (HSS)</td></tr></thead><tbody><tr><td>CLM_XS</td><td>99.9%</td><td>80.9%</td><td>37.6%</td><td>45.6%</td></tr><tr><td>CloudCov</td><td>99.7%</td><td>95.5%</td><td>72.1%</td><td>81.2%</td></tr></tbody></table>

We compare both the CNES CLM mask and the CloudCov mask retrieved from our algorithm to the actual situation. Table 8 summarizes the results for the CNES CLM_XS mask and the CloudCov mask pixel study. While the CLM_XS mask has an HSS score under 50%, the CloudCov mask computed by our algorithm seems to perform very well with a Kappa score above 80%.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Mean confusion matrix. Left: on the CLM_XS (CNES) cloud mask, right: on the CloudCov mask.</figcaption></figure>

### 4.2.2. Temporal Analysis

Figure 12 shows the synthesis of the CloudCov masks computed over the complete VENµS acquisition period and seasonally aggregated. This figure enables the monitoring of cloud behavior through the different seasons of the year. Post-monsoon and winter are the least cloudy seasons, with a mean cloud occurrence of 10%, while monsoon is the cloudiest one with a mean occurrence of 40% (note that the monsoon cloud occurrence is likely underestimated, as fully cloud covered images are not available—118 VENµS images were discarded).

## 4.3. Snow Cover Area Mapping

### 4.3.1. Comparison of Snow Cover with Meteorological Data

For unavailable VENµS data due to high cloud cover, we selected the MAJA chain cloud percentage and filled the cloud cover percentage database retrieved from available images. We compared ground data versus estimations of the extent of snow and cloud covers in Pheriche basin from VENµS images, in order to verify their consistency (Figure 13). We consider here the complete VENµS time period database (November 2018 to June 2020). Figure 13 (top) shows the precipitation at Pyramid site. Outside the rainy season, high cloud cover events are observed intermittently. During the monsoon, periods of high cloud cover persist for several days, whereas outside the rainy season they are more short-lived (Figure 13, bottom). Looking at Pyramid site's precipitation, two phenomena are apparent: (i) the high frequency of pre-monsoon precipitation coincides with the snow height in this period (and so with the cloud cover percentage increase), and (ii) the even greater precipitation frequency during the monsoon. However, the link between cloud cover and precipitation is not obvious and heavy cloud cover does not mean heavy precipitation, or precipitation at all. In addition, the analysis is partial, due to the unavailability of images during the monsoon (high cloud cover percentage), which induce a lack of information on the extent of snow cover. The comparison between cloud cover, snow cover and snow height does not show an obvious link (Figure 13, bottom). Some increases in snow height are related to high cloud cover, but not systematically. Surface snow cover is more persistent in the winter months, with coverage then increasing in the Feb-May pre-monsoon period, before accelerating and then declining during the monsoon. Note that this decline may be an artefact from the strong cloudiness in the monsoon, which prevents detection of the snow cover on the ground. As expected, at Changri Nup site, the snow is higher in the winter season and before the monsoon; the Pheriche basin snow cover percentage follows the same pattern. By joining those data, it seems that for the Pheriche watershed, there is a decrease in the extent of snow cover at the beginning of the monsoon, but likely an increase in the thickness of the snow cover at elevations higher than the glacier equilibrium line, i.e., approx. 5500 m a.s.l.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Cloud cover statistical and temporal occurrence comparison according to seasonal time periods. (a) post-monsoon and winter, (b) pre-monsoon, (c) monsoon, (d) location map. Since we ceased adding to our database on 2 June 2020, only one date of the 2020 monsoon is available; thus, it has been discarded from Figure 12 and only 237 VENµS images are considered here.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** VENµS clouds estimation in Pheriche watershed compared to Pyramid daily cumulative precipitation (**top**), Changri Nup snow height and Pheriche basin interpolated snow coverage. Time period is Nov 2018 to June 2020 (**bottom**). VENµS data availability online is also presented. Snow-covered area has been interpolated in between available Venus images.</figcaption></figure>

### 4.3.2. Seasonal Evolution of SCA versus Cloud Cover at the Watershed Scale

With snow cover area maps (SCA) retrieval, a comparison between the two watersheds of Dingboche and Pheriche is conducted. An example of VENµS SCA maps is presented in Figure A1. For each basin, Figure 14 presents the percentage of the surface area covered by snow or clouds. Although VENµS images with high cloud cover are not available online, the cloud cover percentage has been set to 100% for dates where the CNES' MAJA chain detects 100% cloudiness. The two watersheds offer a similar pattern, coherent with the season. In post monsoon and winter from October to February, there are few clouds and snow coverage is quite constant. In 2018 and 2019, during the pre-monsoon, from April to May, the snow cover progressively decreases and the cloud cover increases. Then, in the monsoon, from June to September, there is little snow cover detected because of the high cloud cover and of the lack of available VENµS images during massive cloud cover episodes.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Time evolution 2018–2020 of snow and cloud covers for (a) Dingboche and (b) Pheriche watersheds. Dashed lines express a linear interpolation.</figcaption></figure>
# 5. Discussion

## 5.1. TopCos HMA Product

As a reminder, we compare in situ albedo measurements (broad band albedo, 0.305–2.8 µm) to narrow reflectance observations (0.865 µm). Overall, none of the reflectance products stands out as an ideal proxy for albedo when compared with field measurements. While all the reflectance products follow the right temporal evolution of albedo, they appear to have biases and performances that are site dependent. This means that empirical calibration is required before using FRE or TopCos HMA reflectances as albedo proxies for applications, such as snow and ice surface energy balance computation. The FRE product seems to be the best product among all sites (Table 7); however, the TopCos HMA product has a better NSS for the Pyramid and South Col sites than for the Changri Nup site (Table A1). Additional reflectance bands in long-wave infrared are required to calculate albedo [68]. The lack of such measurements from the VENµS plateform is a limitation in snow and glaciological applications for energy balance. However, both FRE and TopCos HMA products have satisfactory results with a low RMSE, especially at South Col site: 0.14 for the FRE product and 0.13 for the TopCos HMA product. They show great potential in differentiating snow from rock or bare ice.

Figure 9 highlights the overcorrection in the TopCos HMA reflectance image. A threshold of 0.2 is set on the cos(γ), equivalent to a maximum illumination angle of 78.5°, whereas in Teillet et al. [59], a threshold of 55° is recommended. We made this choice to retain sufficient data, knowing some artifacts and overcorrection could appear. Indeed, an illumination angle close to 90° leads to a high dividing factor and thus to an overcorrection.

Looking at the statistical analysis, the high-resolution HMA DEM does not yield better results than the FRE product corrected with SRTM. In [67], radiometric slope correction using coarse spatial resolution of SRTM offers slightly better results than using fine spatial-resolution DEM from an optical sensor, due to the better vertical accuracy in SRTM retrieved from the active radar (SAR) mode. This means that the cosine correction is likely not the most efficient correction in the rough terrain of the study area. The MAJA processing chains that calculate indirect and ambient reflection (Gamma correction) seem to be better suited (FRE), despite being applied on a lower resolution DEM. It would be worthwhile to test the MAJA chain using HMA DEM as input DEM.

## 5.2. CloudCov Product

Different methods exist in the literature for cloud detection processes and a recent and complete review is given by Nepal et al. [69]. Cloud detection in remote sensing, based on deep learning and machine learning using dense neural networks, actually represents a way of interest to detect clouds (clustering) with better accuracy than existing mask generation methods, i.e., threshold-based methods [70,71]. However, these algorithms need strong calculation architecture that we could not apply for this brief project. This is the reason why we proposed a simple scheme based on a hybrid method which combines co-occurrence textural matrices and radiometric parameters to improve the CLM_XS mask. Another advantage of the CloudCov product is that it enables working at the VENµS ground resolution, since we used a pixel-wise study. Moreover, to our current knowledge, no other study has used VENµS or any other satellite sensor with high spatial and temporal resolution to monitor the spatio-temporal evolution of cloud cover.

In Figure 12, a striking topographic effect is observed in post-monsoon, winter and pre-monsoon at the level of the Mount Everest (Figure 12a,b). Westerly winds blow the snow, creating a snow plume and an orographic banner cloud on the leeward eastern side of Everest's summit, above the Kangshung glacier [72]. During the pre-monsoon, the early stages of the monsoon can be observed with clouds coming from South to North in the low-elevation valleys, and the Everest plume emphasizes this. Cloud intrusions from the South, in the pre-monsoon and monsoon, are constrained to elevation below 5000 m a.s.l. In the pre-monsoon and monsoon, the Tibetan plateau on the leeward side is less impacted by clouds. As a reminder, the quality and reliability of those maps is related to the amount of VENµS images available through the seasons. Note, however, that most of the monsoon
images are not available for analysis here because of the impact of cloud cover (24 images not acquired for technical reason and 65 images discarded in 2018 and 2019 monsoons). To enhance the CloudCov product, an altitudinal threshold could be added in the process for the mapping of low clouds. The step of combining with solar angle could give us the cloud shadow projection on the ground. Thus, cloud shadows would be integrated into the CloudCov product. In addition, we could experiment with this method on the 6th bit of the CLM_XS mask to try to detect thin clouds.

### 5.3. SCA Maps

Many applications using optical remote sensing to retrieve SCA maps in Himalayas are proposed in the literature [73–75], particularly for this Nepalese region using MODIS images [5,14]. The GR4J [73] or the J2000 [75] hydrological models could be improved by using the SCA 5 m maps retrieved from the VENµS dataset thanks to their high spatial resolution. Indeed, the great asset of this SCA database retrieved from VENµS images relies on the high spatial (5 m) and temporal (2-day) resolution of this sensor. These SCA maps will be dedicated to the validation of the snow module from a distributed hydrological model DHSVM-GDM currently applied at IGE [5,13], and actually only based on a daily MODIS 500 m images dataset. The very high temporal resolution of the VENµS satellite is an asset, especially in monsoon conditions, even if there are clouds. Furthermore, such a database is of interest to the community studying the hydro-glaciological context of the region.

In the recent study from Baba & Gascoin [34], an unsupervised classification method has been performed and compared to a supervised classification method to detect snow in VENµS images, showing that a supervised classification gives better results. Since one of the main limitations of this approach is the cloud detection (due to the absence of SWIR band on VENµS sensor) our CloudCov mask, with a 5 m resolution, could be applied ahead of the classification. Then, their method could be tested on this new product, to better discriminate clouds from snow and improve the computed SCA maps.

## 6. Conclusions

The aim of this study was to evaluate and quantify the performance of different products retrieved by the CNES processing of the VENµS database over the Everest region in Nepal for glaciological and hydro-climatic applications. In such a high-mountain context, the focus was set on the interest in adding a fine DEM to the processing chain to enhance the pre-existing quality of the high temporal and spatial resolution from remote data provided by the VENµS sensor. The major findings are the following:

(i) No consistent benefits for assessing the spatio-temporal evolution of surface albedo are retrieved using a cosine radiometric correction enhanced by a fine 5 m DEM regarding the Gamma one carried out with 90 m SRTM by CNES. The CNES FRE product offers efficient results versus in situ measurements for albedo retrieval, because Gamma correction takes into account both direct and diffuse illumination contributions. The cosine correction based only on the direct illumination can induce higher values versus broad band albedo from ground data, even using a precise DEM intended to better reduce the radiometric slope effects.

(ii) We obtained a significant improvement with the TopCos product as an enhancement for cloud cover mapping and satisfactory results for seasonal snow mapping. Our novelty is to compute a hybrid product by merging a radiometric index (NDVI) and a textural approach (Haralick energy matrix) to overtake the CNES initial cloud mask performance.

(iii) Furthermore, SCA maps are also improved since they are obtained using the Cloud-Cov mask.

It appears that no other study using VENµS or any Very High Resolution (VHR) sensor dedicated to monitoring the spatio-temporal evolution of cloud cover has been conducted. Using images usually discarded because of cloudiness to study the spatio-
temporal evolution of clouds in this region could be of interest, in particular to track
cloud cover over the year. In conclusion, the recommendation (for users of future VENµS
missions) is to apply a Gamma pre-processing chain with fine DEM at the same sensor
resolution (i.e., 5 m for VENµS) and a mixed spectral-textural approach on the CNES cloud
mask to improve the cloud masking/tracking in high mountain environments, including
snow and glacial areas.

**Author Contributions:** Conceptualization, J.-P.D. and Y.A.; Methodology, Z.B., J.-P.D. and Y.A.; Software, Z.B.; Validation and Formal Analysis; Z.B., J.-P.D., P.W., Y.A., F.B. and M.E.; Data Curation, Z.B., P.W., M.E., B.P. and T.M.; Writing-Original Draft preparation, Z.B. and J.-P.D.; Writing-Review & Editing, P.W., F.B., Y.A., M.E., B.P. and T.M.; Supervision, J.-P.D.; Project administration P. W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been funded and supported by the French Service of Observation GLACIO-CLIM, now called CRYOBS-CLIM, part of IR OZCAR, and by the French National Research Agency (ANR) through the PRESHINE project, Grant Number ANR-13-SENV-0005.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data collected, pre-processed, processed, or analyzed during this study are included in this work and algorithm codes are available on GitHub at [76].

**Acknowledgments:** This research was part of the selected "Pointing" project from the CNES-ISA VENµS Program, thanks to CNES for the VENµS images collection (KHUMBU site). The authors are grateful to G. Dedieu and O. Hagolle from CESBIO-CNES for their support and advice, S. Pelou from CNES for technical imaging help. Thanks to L. Davaze (IGE) for his help concerning the cosine correction processing and to S. Gascoin (CESBIO-CNES) for the HMA DEM providing. The in situ measurements at Changri Nup and Pyramid sites have been supported by the JEAI HIMALICE Program (PIs: D. Shrestha, and Y. Arnaud), within the framework of the Ev-K2-CNR Project in collaboration with the Nepal Academy of Science and Technology and Tribhuvan University. We are deeply grateful to the National Geographic Society concerning fieldwork measurements at the South Col site: the albedo data collection was conducted in partnership with the National Geographic Society, Rolex, and Tribhuvan University, with approval from all relevant agencies of the Government of Nepal. Special thanks to Aurora Elmore from the National Geographic Society for her review and advice.

**Conflicts of Interest:** The authors declare no conflict of interest.

Appendix A

<figure><img src="image_18.png" /><figcaption>**Figure A1.** Snow-covered area (SCA) mapping for Pheriche and Dingboche watersheds. Example of VENµS image of 2 March 2019.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure A2.** Subset of the SRE product on which the TopCos no data mask is applied. Blue pixels signify no data due to the cos(γ) threshold. VENμS image of 27 March 2018.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure A3.** Mapping comparison of radiometric differences between the four reflectance products. No data values are in white.</figcaption></figure>
<table><caption>**Table A1.** Normalized skill scores (NSS) between SRE, FRE, TopCos HMA, and TopCos SRTM products and in situ albedo AWS data for the three studied sites.</caption><thead><tr><th>Site</th><th>Metrics</th><th>SRE</th><th>FRE</th><th>TopCos HMA</th><th>TopCos SRTM</th></tr></thead><tbody><tr><td rowspan="4">Pyramid</td><td>NSS RMSE</td><td>0.55</td><td>0.64</td><td>0.64</td><td>0.50</td></tr><tr><td>NSS Bias</td><td>0.64</td><td>0.79</td><td>0.79</td><td>0.64</td></tr><tr><td>NSS STD</td><td>0.47</td><td>0.53</td><td>0.53</td><td>0.47</td></tr><tr><td>NSS R<sup>2</sup></td><td>0.95</td><td>0.93</td><td>0.93</td><td>0.95</td></tr><tr><td rowspan="4">Changri Nup</td><td>NSS RMSE</td><td>0.45</td><td>0.45</td><td>0.00</td><td>0.18</td></tr><tr><td>NSS Bias</td><td>0.71</td><td>0.71</td><td>0.00</td><td>0.29</td></tr><tr><td>NSS STD</td><td>0.29</td><td>0.35</td><td>0.00</td><td>0.18</td></tr><tr><td>NSS R<sup>2</sup></td><td>1.00</td><td>1.00</td><td>0.96</td><td>0.97</td></tr><tr><td rowspan="4">South Col</td><td>NSS RMSE</td><td>0.36</td><td>0.36</td><td>0.41</td><td>0.45</td></tr><tr><td>NSS Bias</td><td>0.43</td><td>0.36</td><td>0.86</td><td>0.71</td></tr><tr><td>NSS STD</td><td>0.35</td><td>0.41</td><td>0.24</td><td>0.35</td></tr><tr><td>NSS R<sup>2</sup></td><td>0.93</td><td>0.96</td><td>0.93</td><td>0.96</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Figure A4.** Comparison between daily field-measured albedo at Pyramid site and four reflectance products (85 days). **Top**: SRE vs. TopCos SRTM. **Bottom**: FRE vs. TopCos HMA.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure A5.** Comparison between daily field-measured albedo at Changri Nup site and four reflectance products (155 days). **Top:** SRE vs. TopCos SRTM. **Bottom:** FRE vs. TopCos HMA.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure A6.** Comparison between daily field-measured albedo at South Col site and four reflectance products (84 days). **Top:** SRE vs. TopCos SRTM. **Bottom:** FRE vs. TopCos HMA.</figcaption></figure>
## References

1. Immerzeel, W.W.; Lutz, A.F.; Andrade, M.; Bahl, A.; Biemans, H.; Bolch, T.; Hyde, S.; Brumby, S.; Davies, B.J.; Elmore, A.C.; et al. Importance and vulnerability of the world's water towers. *Nature* **2020**, *577*, 364–369. [CrossRef] [PubMed]

2. Azam, M.F.; Kargel, J.S.; Shea, J.M.; Nepal, S.; Haritashya, U.K.; Srivastava, S.; Maussion, F.; Qazi, N.; Chevallier, P.; Dimri, A.P.; et al. Glaciohydrology of the Himalaya-Karakoram. *Science* **2021**, *373*, eabf3668. [CrossRef] [PubMed]

3. Pritchard, H.D. Asia's shrinking glaciers protect large populations from drought stress. *Nature* **2019**, *569*, 649–654. [CrossRef]

4. Savéan, M.; Delclaux, F.; Chevallier, P.; Wagnon, P.; Gonga-Saholiariliva, N.; Sharma, R.; Neppel, L.; Arnaud, Y. Water budget on the Dudh Koshi River (Nepal): Uncertainties on precipitation. *J. Hydrol.* **2015**, *531*, 850–862. [CrossRef]

5. Mimeau, L.; Esteves, M.; Zin, I.; Jacobi, H.W.; Brun, F.; Wagnon, P.; Koirala, D.; Arnaud, Y. Quantification of different flow components in a high-altitude glacierized catchment (Dudh Koshi, Himalaya): Some cryospheric-related issues. *Hydrol. Earth Syst. Sci.* **2019**, *23*, 3969–3996. [CrossRef]

6. Matthews, T.; Perry, L.B.; Koch, I.; Aryal, D.; Khadka, A.; Shrestha, D.; Abernathy, K.; Elmore, A.C.; Seimon, A.; Tait, A.; et al. Going to Extremes: Installing the World's Highest Weather Stations on Mount Everest. *Bull. Am. Meteorol. Soc.* **2020**, *101*, E1870–E1890. [CrossRef]

7. Hersbach, H.; Bell, B.; Berrisford, P.; Hirahara, S.; Horányi, A.; Muñoz-Sabater, J.; Nicolas, J.; Peubey, C.; Radu, R.; Schepers, D.; et al. The ERA5 global reanalysis. *Q. J. R. Meteorol. Soc.* **2020**, *146*, 1999–2049. [CrossRef]

8. Maussion, F.; Scherer, D.; Mölg, T.; Collier, E.; Curio, J.; Finkelnburg, R. Precipitation seasonality and variability over the Tibetan Plateau as resolved by the high Asia reanalysis. *J. Clim.* **2014**, *27*, 1910–1927. [CrossRef]

9. Bolch, T.; Shea, J.M.; Liu, S.; Azam, F.M.; Gao, Y.; Gruber, S.; Immerzeel, W.W.; Kulkarni, A.; Li, H.; Tahir, A.A.; et al. Status and Change of the Cryosphere in the Extended Hindu Kush Himalaya Region. In *The Hindu Kush Himalaya Assessment*; Wester, P., Mishra, A., Mukherji, A., Shrestha, A., Eds.; Springer International Publishing: Cham, Switzerland, 2019; pp. 209–255.

10. Viviroli, D.; Archer, D.R.; Buytaert, W.; Fowler, H.J.; Greenwood, G.B.; Hamlet, A.F.; Huang, Y.; Koboltschnig, G.; Litaor, M.I.; López-Moreno, J.I.; et al. Climate change and mountain water resources: Overview and recommendations for research, management and policy. *Hydrol. Earth Syst. Sci.* **2011**, *15*, 471–504. [CrossRef]

11. Gascoin, S.; Grizonnet, M.; Bouchet, M.; Salgues, G.; Hagolle, O. Theia Snow collection: High-resolution operational snow cover maps from Sentinel-2 and Landsat-8 data. *Earth Syst. Sci. Data* **2019**, *11*, 493–514. [CrossRef]

12. Eeckman, J.; Chevallier, P.; Boone, A.; Neppel, L.; De Rouw, A.; Delclaux, F.; Koirala, D. Providing a non-deterministic representation of spatial variability of precipitation in the Everest region. *Hydrol. Earth Syst. Sci. Discuss.* **2017**, *21*, 4879–4893. [CrossRef]

13. Mimeau, L.; Esteves, M.; Jacobi, H.W.; Zin, I. Evaluation of gridded and in situ precipitation datasets on modeled glacio-hydrologic response of a small glacierized himalayan catchment. *J. Hydrometeorol.* **2019**, *20*, 1103–1121. [CrossRef]

14. Bouchard, B.; Eeckman, J.; Dedieu, J.P.; Delclaux, F.; Chevallier, P.; Gascoin, S.; Arnaud, Y. On the interest of optical remote sensing for seasonal snowmelt parameterization, applied to the Everest Region (Nepal). *Remote Sens.* **2019**, *11*, 2598. [CrossRef]

15. Pelto, M.; Panday, P.; Matthews, T.; Maurer, J.; Perry, L.B. Observations of Winter Ablation on Glaciers in the Mount Everest Region in 2020–2021. *Remote Sens.* **2021**, *13*, 2692. [CrossRef]

16. Hagolle, O.; Huc, M.; Desjardins, C.; Auer, S.; Richter, R. MAJA Algorithm Theoretical Basis Document. Available online: https://zenodo.org/record/1209633 (accessed on 14 March 2020).

17. Réflectance de surface Venus–THEIA-LAND. Available online: https://www.theia-land.fr/product/venus-2/ (accessed on 9 June 2021).

18. Sirguey, P.; Mathieu, R.; Arnaud, Y. Subpixel monitoring of the seasonal snow cover with MODIS at 250 m spatial resolution in the Southern Alps of New Zealand: Methodology and accuracy assessment. *Remote Sens. Environ.* **2009**, *113*, 160–181. [CrossRef]

19. Sirguey, P.; Still, H.; Cullen, N.J.; Dumont, M.; Arnaud, Y.; Conway, J.P. Reconstructing the mass balance of Brewster Glacier, New Zealand, using MODIS-derived glacier-wide albedo. *Cryosphere* **2016**, *10*, 2465–2484. [CrossRef]

20. Hall, D.K.; Riggs, G.A.; Salomonson, V.V.; DiGirolamo, N.E.; Bayr, K.J. MODIS snow-cover products. *Remote Sens. Environ.* **2002**, *83*, 181–194. [CrossRef]

21. Barrou Dumont, Z.; Gascoin, S.; Hagolle, O.; Ablain, M.; Jugier, R.; Salgues, G.; Marti, F.; Dupuis, A.; Dumont, M.; Morin, S. Brief communication: Evaluation of the snow cover detection in the Copernicus High Resolution Snow & Ice Monitoring Service. *Cryosphere* **2021**, *15*, 4975–4980. [CrossRef]

22. Dumont, M.; Gascoin, S. Optical Remote Sensing of Snow Cover. In *Land Surface Remote Sensing in Continental Hydrology*; Baghdadi, N., Zribi, M., Eds.; Elsevier: Amsterdam, The Netherlands, 2016; pp. 115–137. ISBN 9781785481048.

23. Dozier, J. Spectral signature of alpine snow cover from the landsat thematic mapper. *Remote Sens. Environ.* **1989**, *28*, 9–22. [CrossRef]

24. Warren, S.G. Optical properties of snow. *Rev. Geophys.* **1982**, *20*, 67. [CrossRef]

25. Parajka, J.; Blöschl, G. Validation of MODIS snow cover images over Austria. *Hydrol. Earth Syst. Sci.* **2006**, *10*, 679–689. [CrossRef]

26. Painter, T.H.; Rittger, K.; McKenzie, C.; Slaughter, P.; Davis, R.E.; Dozier, J. Retrieval of subpixel snow covered area, grain size, and albedo from MODIS. *Remote Sens. Environ.* **2009**, *113*, 868–879. [CrossRef]

27. Notarnicola, C.; Duguay, M.; Moelg, N.; Schellenberger, T.; Tetzlaff, A.; Monsorno, R.; Costa, A.; Steurer, C.; Zebisch, M. Snow cover maps from MODIS images at 250 m resolution, part 1: Algorithm description. *Remote Sens.* **2013**, *5*, 110–126. [CrossRef]
28. Dedieu, J.P.; Lessard-Fontaine, A.; Ravazzani, G.; Cremonese, E.; Shalpykova, G.; Beniston, M. Shifting mountain snow patterns in a changing climate from remote sensing retrieval. *Sci. Total Environ.* **2014**, *493*, 1267–1279. [CrossRef]

29. Dozier, J.; Painter, T.H. Multispectral and hyperspectral remote sensing of alpine snow properties. *Annu. Rev. Earth Planet. Sci.* **2004**, *32*, 465–494. [CrossRef]

30. Gascoin, S.; Dumont, Z.B.; Deschamps-Berger, C.; Marti, F.; Salgues, G.; López-Moreno, J.I.; Revuelto, J.; Michon, T.; Schattan, P.; Hagolle, O. Estimating fractional snow cover in open terrain from Sentinel-2 using the normalized difference snow index. *Remote Sens.* **2020**, *12*, 2904. [CrossRef]

31. Kokhanovsky, A.; Lamare, M.; Danne, O.; Brockmann, C.; Dumont, M.; Picard, G.; Arnaud, L.; Favier, V.; Jourdain, B.; Le Meur, E.; et al. Retrieval of Snow Properties from the Sentinel-3 Ocean and Land Colour Instrument. *Remote Sens.* **2019**, *11*, 2280. [CrossRef]

32. Buhler, Y.; Meier, L.; Meister, R. Continuous, High Resolution Snow Surface Type Mapping in High Alpine Terrain Using WorldView-2 Data. Available online: https://www.researchgate.net/publication/267859153_Continuous_highResolution_snowSurface_type_mapping_in_high_alpine_terrain_using_WorldView-2_data (accessed on 10 November 2021).

33. Marchane, A.; Jarlan, L.; Hanich, L.; Boudhar, A.; Gascoin, S.; Tavernier, A.; Filali, N.; Le Page, M.; Hagolle, O.; Berjamy, B. Assessment of daily MODIS snow cover products to monitor snow cover dynamics over the Moroccan Atlas mountain range. *Remote Sens. Environ.* **2015**, *160*, 72–86. [CrossRef]

34. Baba, M.W.; Gascoin, S.; Hagolle, O.; Bourgeois, E.; Desjardins, C.; Dedieu, G. Evaluation of Methods for Mapping the Snow Cover Area at High Spatio-Temporal Resolution with VEN µS. *Remote Sens.* **2020**, *12*, 3058. [CrossRef]

35. Zhu, Z.; Wang, S.; Woodcock, C.E. Improvement and expansion of the Fmask algorithm: Cloud, cloud shadow, and snow detection for Landsats 4–7, 8, and Sentinel 2 images. *Remote Sens. Environ.* **2015**, *159*, 269–277. [CrossRef]

36. Melchiorre, A.; Boschetti, L.; Roy, D.P. Global evaluation of the suitability of MODIS-Terra detected cloud cover as a proxy for Landsat 7 cloud conditions. *Remote Sens.* **2020**, *12*, 202. [CrossRef]

37. Stillinger, T.; Roberts, D.A.; Collar, N.M.; Dozier, J. Cloud Masking for Landsat 8 and MODIS Terra Over Snow-Covered Terrain: Error Analysis and Spectral Similarity Between Snow and Cloud. *Water Resour. Res.* **2019**, *55*, 6169–6184. [CrossRef]

38. Hall, D.K.; Hall, D.K.; Riggs, G.A.; Digirolamo, N.E.; Román, M.O. Evaluation of MODIS and VIIRS cloud-gap-filled snow-cover products for production of an Earth science data record. *Hydrol. Earth Syst. Sci.* **2019**, *23*, 5227–5241. [CrossRef]

39. Zhan, Y.; Wang, J.; Shi, J.; Cheng, G.; Yao, L.; Sun, W. Distinguishing Cloud and Snow in Satellite Images via Deep Convolutional Network. *IEEE Geosci. Remote Sens. Lett.* **2017**, *14*, 1785–1789. [CrossRef]

40. Marais, W.J.; Holz, R.E.; Reid, J.S.; Willett, R.M. Leveraging spatial textures, through machine learning, to identify aerosols and distinct cloud types from multispectral observations. *Atmos. Meas. Tech.* **2020**, *13*, 5459–5480. [CrossRef]

41. Xia, M.; Li, Y.; Zhang, Y.; Weng, L.; Liu, J. Cloud/snow recognition of satellite cloud images based on multiscale fusion attention network. *J. Appl. Remote Sens.* **2020**, *14*, 32609. [CrossRef]

42. Fang, Z.; Ji, W.; Wang, X.; Li, L.; Li, Y. Automatic cloud and snow detection for GF-1 and PRSS-1 remote sensing images. *J. Appl. Remote Sens.* **2021**, *15*, 24516. [CrossRef]

43. Huss, M.; Hock, R. Global-scale hydrological response to future glacier mass loss. *Nat. Clim. Chang.* **2018**, *8*, 135–140. [CrossRef]

44. Maussion, F.; Butenko, A.; Champollion, N.; Dusch, M.; Eis, J.; Fourteau, K.; Gregor, P.; Jarosch, A.H.; Landmann, J.; Oesterle, F.; et al. The Open Global Glacier Model (OGGM) v1.1. *Geosci. Model Dev.* **2019**, *12*, 909–931. [CrossRef]

45. Gardner, A.S.; Sharp, M.J. A review of snow and ice albedo and the development of a new physically based broadband albedo parameterization. *J. Geophys. Res.* **2010**, *115*, F01009. [CrossRef]

46. Izeboud, M.; Lhermitte, S.; Van Tricht, K.; Lenaerts, J.T.M.; Van Lipzig, N.P.M.; Wever, N. The Spatiotemporal Variability of Cloud Radiative Effects on the Greenland Ice Sheet Surface Mass Balance. *Geophys. Res. Lett.* **2020**, *47*. [CrossRef]

47. Wagnon, P.; Vincent, C.; Arnaud, Y.; Berthier, E.; Vuillermoz, E.; Gruber, S.; Ménégoz, M.; Gilbert, A.; Dumont, M.; Shea, J.M.; et al. Seasonal and annual mass balances of Mera and Pokalde glaciers (Nepal Himalaya) since 2007. *Cryosphere* **2013**, *7*, 1769–1786. [CrossRef]

48. Perry, L.B.; Matthews, T.; Guy, H.; Koch, I.; Khadka, A.; Elmore, A.C.; Shrestha, D.; Tuladhar, S.; Baidya, S.K.; Maharjan, S.; et al. Precipitation Characteristics and Moisture Source Regions on Mt. Everest in the Khumbu, Nepal. *One Earth* **2020**, *3*, 594–607. [CrossRef]

49. Khadka, A.; Matthews, T.; Perry, L.B.; Koch, I.; Wagnon, P.; Shrestha, D.; Sherpa, T.C.; Aryal, D.; Tait, A.; Sherpa, T.G.; et al. Weather on Mount Everest during the 2019 summer monsoon. *Weather* **2021**, *76*, 205–207. [CrossRef]

50. Bookhagen, B.; Burbank, D.W. Topography, relief, and TRMM-derived rainfall variations along the Himalaya. *Geophys. Res. Lett.* **2006**, *33*, L08405. [CrossRef]

51. Hagolle, O.; Huc, M.; Pascual, D.V.; Dedieu, G. A multi-temporal method for cloud detection, applied to FORMOSAT-2, VENµS, LANDSAT and SENTINEL-2 images. *Remote Sens. Environ.* **2010**, *114*, 1747–1755. [CrossRef]

52. Ferrier, P.; Crebassol, P.; Dedieu, G.; Hagolle, O.; Meygret, A.; Tinto, F.; Yaniv, Y.; Herscovitz, J. VENµS (Vegetation and environment monitoring on a new micro satellite). In Proceedings of the 2010 IEEE International Geoscience and Remote Sensing Symposium, Honolulu, HI, USA, 25–30 July 2010; IEEE: Piscataway, NJ, USA, 2010; pp. 3736–3739.

53. Meyer, P.; Itten, K.I.; Kellenberger, T.; Sandmeier, S.; Sandmeier, R. Radiometric corrections of topographically induced effects on Landsat TM data in an alpine environment. *ISPRS J. Photogramm. Remote Sens.* **1993**, *48*, 17–28. [CrossRef]

54. Chavez, P.S. Image-based atmospheric corrections—Revisited and improved. *Photogramm. Eng. Remote Sens.* **1996**, *62*, 1025–1036.
55. Naegeli, K.; Damm, A.; Huss, M.; Wulf, H.; Schaepman, M.; Hoelzle, M. Cross-Comparison of Albedo Products for Glacier Surfaces Derived from Airborne and Satellite (Sentinel-2 and Landsat 8) Optical Data. *Remote Sens.* **2017**, *9*, 110. [CrossRef]

56. Lamare, M.; Dumont, M.; Picard, G.; Larue, F.; Tuzet, F.; Delcourt, C.; Arnaud, L. Simulating optical top-of-atmosphere radiance satellite images over snow-covered rugged terrain. *Cryosphere* **2020**, *14*, 3995-4020. [CrossRef]

57. Picard, G.; Dumont, M.; Lamare, M.; Tuzet, F.; Larue, F.; Pirazzini, R.; Arnaud, L. Spectral albedo measurements over snow-covered slopes: Theory and slope effect corrections. *Cryosphere* **2020**, *14*, 1497-1517. [CrossRef]

58. Robledano, A.; Picard, G.; Arnaud, L.; Larue, F.; Ollivier, I. Modelling surface temperature and radiation budget of snow-covered complex terrain. *Cryosph.* **2022**, *16*, 559-579. [CrossRef]

59. Teillet, P.M.; Guindon, B.; Goodenough, D.G. On the Slope-Aspect Correction of Multispectral Scanner Data. *Can. J. Remote Sens.* **1982**, *8*, 84-106. [CrossRef]

60. Riano, D.; Chuvieco, E.; Salas, J.; Aguado, I. Assessment of different topographic corrections in landsat-TM data for mapping vegetation types. *IEEE Trans. Geosci. Remote Sens.* **2003**, *41*, 1056-1061. [CrossRef]

61. Davaze, L.; Rabatel, A.; Arnaud, Y.; Sirguey, P.; Six, D.; Letreguilly, A.; Dumont, M. Monitoring glacier albedo as a proxy to derive summer and annual surface mass balances from optical remote-sensing data. *Cryosph.* **2018**, *12*, 271-286. [CrossRef]

62. Haralick, R.M.; Shanmugam, K.; Dinstein, I. Textural Features for Image Classification. *IEEE Trans. Syst. Man. Cybern.* **1973**, SMC-3, 610-621. [CrossRef]

63. Haralick, R.M. Statistical and structural approaches to texture. *Proc. IEEE* **1979**, *67*, 786-804. [CrossRef]

64. Cohen, J. A Coefficient of Agreement for Nominal Scales. *Educ. Psychol. Meas.* **1960**, *20*, 37-46. [CrossRef]

65. Rosenfield, G.H.; Fitzpatrick-Lins, K. A coefficient of agreement as a measure of thematic classification accuracy. *Photogramm. Eng. Remote Sens.* **1986**, *52*, 223-227.

66. Hüsler, F.; Jonas, T.; Wunderle, S.; Albrecht, S. Validation of a modified snow cover retrieval algorithm from historical 1-km AVHRR data over the European Alps. *Remote Sens. Environ.* **2012**, *121*, 497-515. [CrossRef]

67. Wu, Q.; Jin, Y.; Fan, H. Evaluating and comparing performances of topographic correction methods based on multi-source DEMs and Landsat-8 OLI data. *Int. J. Remote Sens.* **2016**, *37*, 4712-4730. [CrossRef]

68. Dumont, M.; Gardelle, J.; Sirguey, P.; Guillot, A.; Six, D.; Rabatel, A.; Arnaud, Y. Linking glacier annual mass balance and glacier albedo retrieved from MODIS data. *Cryosphere* **2012**, *6*, 1527-1539. [CrossRef]

69. Mahajan, S.; Fataniya, B. Cloud detection methodologies: Variants and development—A review. *Complex Intell. Syst.* **2020**, *6*, 251-261. [CrossRef]

70. Shi, M.; Xie, F.; Zi, Y.; Yin, J. Cloud detection of remote sensing images by deep learning. In *Proceedings of the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)*, Beijing, China, 10-15 July 2016; IEEE: Piscataway, NJ, USA, 2016; pp. 701-704.

71. Li, Z.; Shen, H.; Cheng, Q.; Liu, Y.; You, S.; He, Z. Deep learning based cloud detection for medium and high resolution remote sensing images of different sensors. *ISPRS J. Photogramm. Remote Sens.* **2019**, *150*, 197-212. [CrossRef]

72. Moore, G.W.K. Mount Everest snow plume: A case study. *Geophys. Res. Lett.* **2004**, *31*, 1-4. [CrossRef]

73. Nepal, S.; Chen, J.; Penton, D.J.; Neumann, L.E.; Zheng, H.; Wahid, S. Spatial GR4J conceptualization of the Tamor glaciated alpine catchment in Eastern Nepal: Evaluation of GR4JSG against streamflow and MODIS snow extent. *Hydrol. Process.* **2017**, *31*, 51-68. [CrossRef]

74. Shukla, S.; Kansal, M.L.; Jain, S.K. Snow cover area variability assessment in the upper part of the Satluj River Basin in India. *Geocarto Int.* **2017**, *32*, 1285-1306. [CrossRef]

75. Nepal, S. Impacts of climate change on the hydrological regime of the Koshi river basin in the Himalayan region. *J. Hydro-Environment Res.* **2016**, *10*, 76-89. [CrossRef]

76. GitHub-Bessinz/VENuS_Cosine_Correction. Available online: https://github.com/bessinz/VENuS_cosine_correction (accessed on 21 February 2022).