Satellite-derived products of solar and longwave irradiances used for snowpack modelling in mountainous terrain

In mountainous terrain, the snowpack is strongly affected by incoming shortwave and longwave radiation. In this study, a thorough evaluation of the solar and longwave downwelling irradiance products (DSSF and DSLF) derived from the Meteosat Second Generation satellite was undertaken in the French Alps and the Pyrenees. The satellitederived products were compared with forecast fields from the meteorological model AROME and with analysis fields from the SAFRAN system. A new satellite-derived product (DSLFnew) was developed by combining satellite observations and AROME forecasts. An evaluation against in situ measurements showed lower errors for DSSF than AROME and SAFRAN in terms of solar irradiances. For longwave irradiances, we were not able to select the best product due to contrasted results falling in the range of uncertainty of the sensors. Spatial comparisons of the different datasets over the Alpine and Pyrenean domains highlighted a better representation of the spatial variability of solar fluxes by DSSF and AROME than SAFRAN. We also showed that the altitudinal gradient of longwave irradiance is too strong for DSLFnew and too weak for SAFRAN. These datasets were then used as radiative forcing together with AROME near-surface forecasts to drive distributed snowpack simulations by the model Crocus in the French Alps and the Pyrenees. An evaluation against in situ snow depth measurements showed higher biases when using satellite-derived products, despite their quality. This effect is attributed to some error compensations in the atmospheric forcing and the snowpack model. However, satellite-derived irradiance products are judged beneficial for snowpack modelling in mountains, when the error compensations are solved.

1 Introduction 20 Seasonal snowpacks are a key component of mountain hydrological systems. Snow accumulation and ablation processes set up the temporal evolution of the snow cover and its spatial distribution, controlling the snow melt variability and timing, which govern the run-off in high-altitude catchments (e.g. Anderton et al., 2002;DeBeer and Pomeroy, 2017). The evolution and spatial distribution of the snowpack in mountainous terrain depends on its energy budget, affected by the surface 25 radiative budget, the sensible and latent heat fluxes and the ground heat flux (e.g. Armstrong and Brun, 2008). The meteorological conditions are the main factors controlling the snow surface energy budget, with a key contribution of the radiative components (Male and Granger, 1981). For example, Cline (1997) reported a contribution of 75% of net radiative fluxes in the energy for snowmelt over the entire season at a continental midlatitude alpine site of the Colorado Front Range (3517 m), 30 while Marks and Dozier (1992) found a contribution between 66% and 90% at two alpine sites of the Sierra Nevada (2800 m and 3416 m). Therefore, incoming shortwave (SW↓) and longwave (LW↓) radiative fluxes are amongst the most significant atmospheric factors of the energy and mass budget of the snowpack, particularly during snowmelt periods. It is crucial to accurately represent them in numerical snowpack simulations, as recent works underlined the strong sensitivity of snowpack 35 Jura. The Pyrenees domain covers the latitudes from 41.6 • N to 43.6 • N and the longitudes from -2.5 • E to 3.5 • E. Hourly data, from 1 August 2010 to 31 July 2014, including in situ measurements, satellite-derived irradiance products, meteorological models and snowpack simulations were used.

Irradiance datasets 95
Several irradiance datasets were used in this study: forecasts from the NWP model AROME, reanalyses from the SAFRAN analysis system, LSA SAF irradiance products derived from remotely-sensed observations and a hybrid LW↓ irradiance product based on a combination of LSA SAF algorithms with AROME forecasts. An in situ observation dataset was built up for validation in mountains.
2.2.1 NWP system: AROME 100 AROME (Application of Research to Operations at MEsoscale) is the meso-scale NWP system of Météo-France (Seity et al., 2011), operating over France since December 2008 at 2.5 km grid spacing (1.3 km since 2015; Brousseau et al., 2016). It is a spectral and non-hydrostatic model.
The physics and data assimilation schemes are detailed in Seity et al. (2011). In particular, AROME uses the radiation parameterisations from the European Centre for Medium-Range Weather Forecasts 105 (ECMWF), with the SW scheme from Fouquart and Bonnel (1980) and the LW scheme from Mlawer et al. (1997).
In this study, we built a continuous atmospheric forcing dataset to drive snowpack simulations using hourly AROME forecasts issued from the 0 UTC analysis time, from + 6 h to + 29 h, extracted on a regular latitude/longitude grid with a 0.025 • resolution over the period and domains of study soundings. In particular, the incoming shortwave and longwave fluxes are computed with the radiation scheme from Ritter and Geleyn (1992), using as first guess vertical profiles of temperature and humidity from ARPEGE forecasts, atmospheric soundings, a guess of cloudiness based on the analysed vertical humidity profile and a cloud mask detected by satellite (Derrien et al., 1993).
In this study, we used SAFRAN reanalyses from 1 August 2010 to 31 July 2014. For compar-product also derived from MSG data (Stöckli, 2013;Castelli et al., 2014), the DSSF is not down-scaled over complex terrain, and thus not corrected for local topography effects. The target accuracy of the DSSF is 10% or 20 W m −2 for values lower than 200 W m −2 (Trigo and Viterbo, 2009). 165 • LW↓ irradiance: DSLF The algorithm to estimate the DSLF is described in details by Trigo et al. (2010). It consists in a modified version of the bulk parameterisation of Prata (1996), initially developed for clear skies only. It relies on a formulation of the effective emissivity and temperature of the atmospheric layer above the surface, using the TCWV, 2 m temperature (T 2m ) and 2 m dew 170 point (Td 2m ) forecast by the ECMWF IFS. The formulation parameters are calibrated for clear-sky and overcast conditions independently. The MSG/SEVIRI cloud mask (Derrien and Le Gléau, 2005) is thus the only observation used, to distinguish clear-sky and cloudy-sky situations. In case of partly cloudy situations, the average of both terms is taken. The DSLF can therefore be described more accurately as a longwave irradiance parameterisation using 175 satellite observations of the cloud mask rather than a satellite product. The DSLF is not downscaled over complex terrain, and thus not corrected for local topography effects. The target accuracy of the DSLF is 10% (Trigo and Viterbo, 2009).

New DSLF product using AROME forecasts
The DSLF relies on the ECMWF IFS forecasts of TCWV, T 2m and Td 2m . These atmospheric vari-180 ables have a strong dependence on altitude and a strong spatial variability in mountainous terrain.
The 16-km horizontal resolution of the ECMWF IFS hardly represents this spatial variability in the Alps and the Pyrenees, despite a constant lapse rate applied for grid elevation correction. Consequently, we developed a new DSLF product using the same algorithm (Trigo et al., 2010) depending on the cloud mask (Derrien and Le Gléau, 2005), but replacing ECMWF forecasts by AROME fore-forecasts. The algorithm was applied to the new DSLF on the LSA SAF grid over the domains of study ( Fig. 1), from 1 August 2010 to 31 July 2014. Hereafter, this product is referred to as DSLFnew.

In situ irradiance observations
To assess the distributed irradiance datasets, ground measurements of SW↓ and LW↓ were extracted 200 from Météo-France station network and additional Automatic Weather Stations (AWS). Stations with altitude higher than 1000 m were selected. As elevation influences incoming radiation (Oliphant et al., 2003), stations were not used for evaluation if the difference between the station elevation and the elevation of the four closest AROME and LSA SAF grid points was higher than 300 m.
As elevation differences up to 300 m may have an influence on the comparisons, the altitudes of 205 the grid points associated with each station are listed in Table 1 Table 1.
Irradiance measurements are scarce in mountainous terrain and their quality is often lower than plain measurements, due to the difficulty to maintain these stations and the possible occurrence of frost or snow on the sensors in winter (Lapo et al., 2015a). The pyranometers from Météo-France network (Kipp&Zonen CM5, CM6B and CM11) meet the good quality standards of the World Me-215 teorological Organization (WMO, 2014), hence an uncertainty of hourly total SW↓ irradiance of ±10% (Leroy and Leches, 2014). Due to their location in altitude, the maintenance may not be systematically weekly so that uncertainties of ±10% are probably too optimistic. The station of Carpentras in plains is equipped with the pyranometer Kipp&Zonen CM21 and the pyrgeometer Kipp&Zonen CG4. This station is a reference station for radiation measurements, as it is part of 220 the Baseline Surface Radiation Network (BSRN; Ohmura et al., 1998) (Szczypta et al., 2015), Argentière glacier and St-Sorlin glacier (data from GLACIOCLIM program, 225 https://glacioclim.osug.fr) have Kipp&Zonen CM3 pyranometers and CG3 pyrgeometers, classified as moderate quality after WMO's standards (WMO, 2014), for which the manufacturer reports a daily total accuracy of ±10%. The uncertainties have not been estimated at these stations. They are possibly higher than 10% because of the difficulty to maintain AWS in complex environment, particularly in winter. WMO (2014) indicates uncertainties up to ±20% for hourly totals for this kind 230 of instruments. The results at these stations are indicative for high altitudes but shall be considered carefully. Table 1 summarizes the measurement uncertainties at each station. SW↓ and LW↓ irradiances from LSA SAF products, AROME forecasts and SAFRAN reanalyses were evaluated using these in situ measurements. The altitude of the grid points associated to each station is reported in Table 1  The SW↓ irradiance products were only evaluated when the sun was above the horizon, or when the observed value was higher than 20 W m −2 at Andorre and Envalira stations (to discard periods when the sun is masked by the terrain). The LW↓ irradiance products were evaluated by day and night.

Snowpack datasets
The impact of the different irradiance datasets on distributed snowpack simulations is assessed using 245 the snowpack model Crocus with different atmospheric forcings. These simulations are compared to in situ measurements of snow depth.

Snowpack model: Crocus
Snowpack simulations driven by different irradiance datasets were performed with the detailed snow cover model Crocus (Brun et al., 1992;Vionnet et al., 2012) coupled with the ISBA land surface 250 model within the SURFEX simulation platform (Masson et al., 2013), to fully simulate the interactions between snowpack and soil. SURFEX/ISBA-Crocus (called Crocus hereafter) simulates the evolution of the snowpack physical properties along its stratigraphy, under given atmospheric forcing data (temperature and specific humidity at a given height above the surface, wind speed at a given height above the surface, SW↓ and LW↓ irradiance, solid and liquid precipitation).

255
The simulations were carried out over the French Alps and Pyrenees domains ( Fig. 1), on the AROME regular latitude/longitude grid at 0.025 • resolution (Sect. 2.2.1) from 1 August 2010 to 31 July 2014. The effects of aspect and slope on incoming solar irradiance were not taken into account, because the snowpack is simulated over flat terrain, and the interactions with the vegetation and the parameterisation of fractional snow cover were not activated, because the evaluation observations 260 are located in flat and open fields. This configuration has already been used in Vionnet et al. (2016) and Quéno et al. (2016).
Except incoming radiative fluxes, the atmospheric forcing of the snowpack simulations was built with AROME forecasts (Sect. 2.2.1). The radiative components of the forcings were extracted from the different irradiance datasets: a) AROME irradiance forecasts (simulations named A-Cro 265 hereafter), b) SAFRAN irradiance reanalyses (simulations named AS-Cro hereafter), c) DSSF and DSLFnew (simulations named AL-Cro hereafter), d) DSSF and AROME LW↓ irradiance (simula-tions named AL SW -Cro hereafter), e) DSLFnew and AROME SW↓ irradiance (simulations named AL LW -Cro hereafter). In order to include DSSF and DSLFnew products in AROME forcing, the interpolation on AROME grid was made to minimize the effect of elevation difference on the incoming 270 radiative fluxes. Among the four nearest LSA SAF grid points, the grid point with the minimum altitude difference with AROME grid point was chosen. Similarly to Hinkelman et al. (2015), SW↓ irradiances were not modified, whereas a vertical gradient of -29 W m −2 km −1 (Marty et al., 2002) was applied to LW↓ irradiances to mitigate the remaining differences in altitude. The different simulations are summarized in Table 2.

In situ snowpack observations
To assess the quality of Crocus simulations, an observational dataset of snow depth measurements was constituted in the French Alps and the Pyrenees, within SAFRAN massifs. Only stations with less than 150 m elevation difference to the model topography were selected, in order to use the same dataset as Quéno et al. (2016) and Vionnet et al. (2016). This dataset contains a total of 172 stations 280 (89 in the French Alps and 83 in the Pyrenees) with daily manual measurements at ski resorts (at 6 UTC) and daily automatic measurements by ultra-sonic sensors at high altitude sensors, as described in details in Vionnet et al. (2016) for the French Alps and in Quéno et al. (2016) for the French and Spanish Pyrenees.
3 Evaluation of irradiance products over the Alps and the Pyrenees 285

Comparisons with in situ measurements
The SW↓ error statistics for all stations are listed in Table 1 Overall, DSSF exhibits the best performance with a relative bias of -4% and a relative RMSE of 33%. SAFRAN has a relative bias of -7% and a relative RMSE of 40%. Finally, AROME exhibits the strongest relative bias (+ 12%) and the highest relative RMSE (43%). 300 Fig. 3 shows biases and RMSE of the different datasets of incoming LW↓ (DSLF, DSLFnew, AROME and SAFRAN) at the five LW↓ stations and the overall error statistics. In this figure, stations are ordered by altitude. In mountains, DSLF, DSLFnew and AROME have a negative bias, while SAFRAN bias tends to increase with altitude. At low elevation (Carpentras), the best performance is in favour of DSLFnew with a bias of + 4 W m −2 (+ 1%) and a RMSE of 16 W m −2 (5%), which 305 falls within the range of uncertainties of the sensor. At three mountain stations (Col de Porte, Bassiès and Argentière glacier), the lowest bias and RMSE are reached by SAFRAN, while AROME has the lowest RMSE at St-Sorlin glacier. Overall, AROME exhibits the strongest negative relative bias (-6%) and the highest relative RMSE (12%). DSLF and DSLFnew have equivalent error statistics with a relative bias of -3% and a relative RMSE of 11%. Finally, SAFRAN has a relative bias of + 1% 310 and a relative RMSE of 11%. These global error statistics are close to the sensor uncertainties in mountains, which does not enable to choose the "best product". However, some trends are identified such as an underestimation of LW↓ by DSLF, DSLFnew and AROME. The performance of LSA SAF products and models is also clearly better in terms of LW↓ than SW↓, because of lower biases and RMSE.

315
The yearly cycles of SW↓ irradiances are illustrated at Carpentras for reference (Fig. 4a) and at Péone mountain station (Fig. 4b). They show higher RMSE in Spring and Summer for each dataset, lowest RMSE for DSSF and highest RMSE for AROME during the whole year, except in December and January where the three products have equivalent RMSE. This trend was found similar at all stations. No specific trend was observed for the bias. The SW↓ daily cycles ( Fig. 4c for Carpentras 320 and Fig. 4d for Péone) show a lower RMSE for DSSF in the middle of the day. SAFRAN cycle is not marked enough (positive biases in the morning and evening, negative biases in the middle of the day). AROME overestimates SW↓ all the time. DSSF represents well the diurnal cycle, with an underestimation in the afternoon. These trends were also highlighted at the other mountain stations.
The study of the daily and yearly cycles of LW↓ irradiances did not indicate any particular trend for 325 metrics following the month or the hour (not shown).

Spatial comparisons of the distributed products
Spatial comparisons of the different irradiance products were carried out over the two domains.
DSSF and DSLF were taken as references. The spatial distributions of their annual mean computed using data from 1 August 2010 to 31 July 2014 and the differences with the other irradiance products 330 are shown in Fig. 5 for the French Alps and in Fig. 6 for the Pyrenees.
The DSLF exhibits a strong correlation with the altitude, with a decreasing LW irradiance towards the highest elevations, i.e. the East of the French Alps (Fig. 5a) and the central range of the Pyrenees (Fig. 6a). AROME presents a moderate negative bias as compared to the DSLF, both in the Alps ( Fig. 5b) and in the Pyrenees (Fig. 6b), while SAFRAN presents a strong positive bias, particularly 335 in the highest areas of the Alps (Fig. 5c) and the Pyrenees (Fig. 6c). DLSFnew presents a slight positive bias over most of the domains, except over the highest peaks where the bias is slightly negative ( Fig. 5d and Fig. 6d).
The DSSF exhibits a lower correlation with the topography (Fig. 5e and Fig. 6e). For given sky conditions, the SW irradiance increases with the elevation as the atmospheric transmissivity in-340 creases. But the annual mean of the DSSF follows more regional patterns of cloud cover than elevation patterns. For example, in the French Alps, Fig. 5e shows a North-West -South-East gradient of increasing DSSF: South-Eastern massifs are often shielded by North-Western massifs in the most frequent case of West and North-West disturbed flows. A similar gradient of precipitation was shown in Durand et al. (2009b). The heterogeneity of DSSF is even more marked in the Pyrenees (Fig. 6e)

345
where the West-East chain acts as an orographic barrier to the prevailing northwesterlies coming from the Atlantic Ocean . A clear discontinuity appears between the French Pyrenees, where the clouds are often blocked, and the Spanish Pyrenees, often affected by Foehn wind and resulting clear sky conditions. The lowest DSSF are found in the Western part of the French Pyrenees, while the Eastern part is more sunny due to the abating Atlantic influence and a 350 Mediterranean climate. AROME presents a strong positive bias ( Fig. 5f and Fig. 6f), locally higher than 30% over the highest peaks, and still higher than 15% in many plain areas. SAFRAN bias is very variable from one massif to another ( Fig. 5g and Fig. 6g). A strong negative bias for SAFRAN can be noticed in the South-Western massifs of the Spanish Pyrenees (Fig. 6g), highlighting a poor representation of the orographic blocking as already noticed in Quéno et al. (2016).

355
The dependence of the different irradiance products with the altitude was further explored with the study of altitudinal gradients. Figure 7 represents the vertical evolution of the LW↓ and SW↓ averaged over the SAFRAN massifs of the French Alps and the Pyrenees by steps of 100 m of elevation over the whole study period, together with the associated standard deviations.
The strong dependency of LW↓ irradiance with altitude is confirmed in Fig. 7a for the French 360 Alps and Fig. 7b for the Pyrenees. As a reference, the altitudinal gradient for annual LW↓ means of -29 W m −2 km −1 found by Marty et al. (2002) in the Swiss Alps is plotted in dashed line, while Table 3 lists the mean altitudinal gradient for each dataset in both domains. All datasets present a steady decrease of LW↓ with altitude, and are close to each other below 1200 m approximately. For higher elevations, SAFRAN annual mean value is significantly stronger than AROME, DSLF and 365 DSLFnew, due to a lower vertical gradient (Table 3). We showed in Sect 3.1 that AROME, DSLF and DSLFnew had a negative bias at the four mountain stations. This effect may come from a too strong vertical gradient (Table 3). DSLFnew is larger than AROME and DSLF at all altitudes below 2900 m in the French Alps (Fig. 7a) and 2200 m in the Pyrenees (Fig. 7b) approximately. It gets lower at the highest altitudes due to a stronger vertical gradient. The stronger vertical gradient of 370 DSLFnew compared to DSLF is the confirmation that the use of forecasts of higher resolution for the algorithm takes more into account the topography. The excessive vertical gradient may originate from the cold bias of AROME near-surface temperatures, enhanced with the altitude , leading to a strong underestimation of the fluxes by DSLFnew at the highest altitudes.
In terms of SW↓ irradiance, Fig. 7c and Fig. 7d highlight that AROME fluxes are significantly 375 stronger than SAFRAN and DSSF at all altitudes. SAFRAN is marked by an increase of incoming SW↓ fluxes with altitude, while AROME and DSSF present a more variable evolution, and particularly a decrease of the fluxes at the highest altitudes in the Alps (Fig. 7c). This decrease may reflect the more frequent presence of clouds blocked by the highest peaks. Furthermore, these figures underline a weaker dependency of SW↓ irradiance with altitude than LW↓ irradiance. Indeed, 380 the standard deviation of LW↓ at a given altitude is small compared to the total variation of the mean LW↓ with altitude for all products (Fig. 7a and Fig. 7b), whereas they can reach similar values for SW↓ ( Fig. 7c and Fig. 7d). This spatial variability at a given altitude is particularly marked at lowand mid-altitudes (< 1800 m) in the Pyrenees for AROME and DSSF, reflecting a good representation of the strong climate heterogeneity between French and Spanish foothills. SAFRAN, which 385 gives homogeneous analyses per massif, does not account for the spatial variability within the massif as is the case for AROME and DSSF.

Impact of the irradiance products on snowpack simulations
Snowpack simulations were performed over four winters from 2010 to 2014 to assess the impact of the different irradiance datasets as radiative forcing. Table 4   higher than the other datasets. In winter, SW↓ irradiances are low and the snow albedo is high: their contribution to the surface energy budget is much lower than in spring. Thus, LW↓ irradiances have a higher relative contribution during the accumulation period. However, during the melting period (from mid-March to mid-May here), the contribution of SW↓ irradiances is the highest, due to higher extra-terrestrial solar fluxes, longer days and lower snow albedo: because of their higher SW↓, A-Cro 415 simulations melt faster than AS-Cro, which reduces their bias.
These trends can also be observed when looking at maps of spatially distributed snowpack simulations. Figure 9 represents the SWE (snow water equivalent) simulated by A-Cro taken as a reference on 1 February 2013 during the accumulation period and on 1 May 2013 during the melting period, and the differences between AL-Cro, AS-Cro and this reference at the same dates. The differences 420 with AL-Cro are generally between -50 mm and + 50 mm on 1 February 2013. AS-Cro exhibits lower SWE values at this date, due to its higher LW↓ irradiance. However, on 1 May 2013, both simulations exhibit higher SWE values than A-Cro almost everywhere, with differences mostly higher than 200 mm, locally reaching 400 mm, due to lower SW↓ irradiances.
The impact of the radiative forcing on SWE simulations was further studied at two grid points with AROME as reference. The same evolutions at point B are represented in Fig. 11. The relative impact of DSSF and DSLFnew is represented in dashed lines (simulations AL SW -Cro and AL LW -Cro, as defined in Table 2). At point A, melting occurs during the winter. Consequently, AS-Cro and AL LW -Cro simulations lead to lower values of SWE than A-Cro, since they both exhibit higher LW↓ than AROME (+ 8 W m −2 for DSLFnew and + 9 W m −2 for SAFRAN). Thus, on 15 February 435 by the higher LW↓. A marked difference in the melt timing can be noted for AL-Cro: the lower SW↓ 445 is not counterbalanced by the slightly higher LW↓. The peak SWE is shifted by almost one month compared to A-Cro. Therefore, it leads to marked differences in terms of cumulated melting: on 1 June 2011, the cumulated melting for A-Cro reaches 1149 mm, i.e. almost the double of AL-Cro (613 mm, and 433 mm for AL SW -Cro). The simulation mixing DSSF and DSLFnew irradiances (AL-Cro) is very close to the DSSF-only simulation (AL SW -Cro). Overall, the effect of DSSF pre-450 vails at high altitude leading to a later end of the snow cover, while the effect of DSLFnew prevails at low altitude leading to an earlier end of the snow cover.

Quality of irradiance datasets in mountainous terrain
We presented an overview of the quality of several irradiance datasets through an in-depth assess-455 ment of the irradiance fields in mountainous terrain. In terms of SW↓ irradiances, DSSF exhibits best metrics in mountains, particularly below 2000 m. Above 2000 m, its RMSE is similar to SAFRAN and AROME, due to a strong negative bias. AROME presents systematic and large overestimations of SW↓ irradiances, contrarily to SAFRAN tendency to underestimate them. The spatial variations of SW↓ irradiances are better represented in DSSF and AROME than in SAFRAN. In terms of LW↓ 460 irradiances, the obtained errors are comparable and it is difficult to identify the best product. The use of forecasts at higher spatial resolution to compute DSLFnew enhances the topographic dependence, which limits the underestimation of LW↓ irradiance at low and mid-altitudes found with DSLF, but strengthens the negative bias at high altitude. The resulting altitudinal gradient is probably too strong.
It may originate from the cold bias of AROME near-surface temperatures, enhanced with the alti-465 tude , which leads to a strong underestimation of the fluxes by DSLFnew at the highest altitudes.
Several studies evaluated LSA SAF irradiance products at hourly time step (when the sun is above the horizon for SW↓) at plain stations. For DSSF, we showed in this study a bias of -14 W m −2 and a RMSE of 117 W m −2 (Fig. 2), while in plains, Geiger et al. (2008b), Ineichen et al. (2009) 470 and Cristóbal and Anderson (2013) reported biases of + 2 W m −2 , + 5 W m −2 and -5 W m −2 respectively, and RMSE of 87 W m −2 , 103 W m −2 and 65 W m −2 respectively. The higher RMSE in mountains may partly be explained by higher mean values. For DSLF, we showed in this study a bias of -8 W m −2 and a RMSE of 32 W m −2 (Fig. 3), while in plains, Trigo et al. (2010) and Ineichen et al. (2009) reported biases of + 3 W m −2 and -11 W m −2 respectively, and RMSE of 25 475 W m −2 and 29 W m −2 respectively. The error statistics in mountains are close to those in plains, and lie within the range of uncertainty of LW↓ sensors in mountains (Table 1). Thus, the performance of LSA SAF irradiance products remains satisfactory compared to previous evaluations of these products in plains, even though they generally do not reach the target accuracy (Sect. 2.2.3), derived from reference plain stations. radiative flux variability, by order of importance: slope aspect, slope angle, elevation, albedo, shad-515 ing, sky view factor, and leaf area index. These local factors are not taken into account in AROME, SAFRAN and LSA SAF irradiance products. This study aims at assessing the practical benefits of different irradiance datasets to be used as radiative forcing for distributed snowpack simulations at 2.5 km resolution in mountains. In the context of representing the mean state of the snowpack over a considered flat pixel, at a given altitude and a given location in the mountain range, the terrain 520 influence on the radiation does not need to be taken into account in the radiative forcing. However, to capture the sub-kilometric variability of the snowpack, it will be necessary to consider sub-grid effects of the surrounding terrain on the radiation, and thus a topographical correction of irradiance products (e.g. Helbig and Löwe, 2012) as done for MSG satellite-derived solar fluxes by the HelioMont method (Stöckli, 2013;Castelli et al., 2014). shadowing effect on direct solar radiation. However, this effect is not considered for snow depth comparisons. Additionally, the limited sky view and the reflection effects on diffuse solar radiation are not taken into account, as well as the limited sky view and terrain thermal radiation effects on longwave irradiance. According to the incoming radiation budget accounting for topography effects on a flat surface (e.g. Müller and Scherer, 2005), the bias for not taking into account the sky view 535 factor (SVF) in the shortwave irradiance is (1 − SV F ) * (SW dif − SW ref ), where SW dif is the diffuse shortwave irradiance and SW ref is the terrain-reflected shortwave radiation. This bias is small in clear-sky conditions, because the SVF contribution is then significantly lower than topographic shading (e.g. Olson et al., 2019), but it increases in cloudy conditions with a stronger relative contribution of SW dif . The bias for not taking into account the SVF in the longwave irradiance is

540
(1 − SV F ) * (LW atm − LW ter ), where LW atm is the atmospheric longwave irradiance and LW ter is the terrain thermal radiation. This bias is more significant when LW atm is low (Sicart et al., 2006), but is usually small compared to LW atm on horizontal surfaces (Plüss and Ohmura, 1997), which is the case of the measurement locations.

545
DSSF and DSLFnew irradiance datasets were used to replace AROME irradiance forecasts as radiative forcing of Crocus simulations. The rest of the atmospheric forcing was taken from AROME forecasts. A similar experiment was done with SAFRAN irradiances. The performance of the snowpack simulations was degraded when using DSSF and DSLFnew products, with an increased positive snow depth bias. On the contrary, the use of SAFRAN irradiances was found to decrease the 550 positive bias obtained with AROME-Crocus. Vionnet et al. (2016) and Quéno et al. (2016) already showed an overestimation of snow depth by AROME-Crocus in the French Alps and the Pyrenees respectively. In addition, Quéno et al. (2016) partly attributed this overestimation to an underestimation of strong melting. Thus, replacing AROME irradiance forecasts by lower or equivalent values (DSSF and DSLFnew) logically enhances the overestimation, despite the better quality of the new irradiance products. In this case, improving the radiative forcing leads to degraded snowpack simulations, which is consistent with the strong interaction effect in snowpack simulations highlighted by Günther et al. (2019). This effect may be attributed to error compensations within the atmospheric forcing and/or within the snowpack model: - Quéno et al. (2016) showed that the positive snow depth bias is not due to an overestimation of 560 daily snowfall by AROME-Crocus. The strong overestimation of SW↓ by AROME shown in this study would also tend to increase the melting and reduce the snow depth bias. We showed here it is not counterbalanced by the underestimation of LW↓. However, the underestimated melting may be linked to an underestimation of the turbulent fluxes. It may have several causes that need to be further explored, e.g. a possible influence of the T 2m cold bias, particularly 565 marked at the highest altitudes (-2.8 K above 2500 m ; Vionnet et al., 2016).
-Within the snowpack model Crocus, Quéno et al. (2016) showed an underestimation of snow settling, with a direct effect on snow depth bias. The snow compaction scheme used in Crocus depends on the weight of overlying snow, temperature and density of snow, liquid water content and snow grain size (Vionnet et al., 2012). The parameterisation of the albedo evolu-  (Raleigh et al., 2015;Lapo et al., 2015b;Sauter and Obleitner, 2015). First, Sauter and Obleitner (2015) studied the influence of uncertainties on atmospheric forcing variables on simulations of glacier mass-balance using Crocus in the Svalbard islands (European Arctic). They identified LW↓ uncertainty as the main source of variance (50%) of the surface energy balance throughout the year.
However, the prevailing effect of LW↓ compared to SW↓ is more marked at high latitudes, because 585 of the lack of solar insolation in winter. In our study, we showed that the new LW↓ forcing from DSLFnew (with a positive bias compared to AROME) had a significant impact on the mass budget during the whole winter at low altitudes (Fig. 10), while the impact was more limited at high altitudes ( Fig. 11). It can be explained by decreasing LW↓ irradiances with altitude together with increasing SW↓ irradiances, leading to a more significant impact of SW↓ at high altitudes. It is also due to the 590 earlier snowmelt at low altitudes, which limits the crucial role played by SW↓ in spring.
Furthermore, the differences between the different radiative forcing datasets mainly consist of biases rather than random errors: a typical example is the difference of SW↓ at high altitudes between AROME and DSSF shown in Fig. 11. Their effect is then cumulated during the whole season, rather than counterbalanced, which increases their impact. It is consistent with the outcomes of Raleigh 595 et al. (2015) who showed that snowpack models are more sensitive to biases than random errors in the forcings. It was particularly highlighted for incoming radiative fluxes by Lapo et al. (2015b).
Finally, although the SWE is not impacted by the differences in incoming radiative fluxes at high altitude during the accumulation period (Fig. 11), impacts are to be expected in terms of snow surface temperatures, with possible consequences on the snow metamorphism processes. Lapo et al. (2015b) 600 indeed showed more sensitivity of the snowpack simulations to irradiance errors at the coldest sites when evaluated in terms of snow surface temperature rather than SWE. Future works could thus focus on the impact of the different incoming radiative flux datasets on the surface energy budget and the resulting effects on the snowpack stratigraphy.

605
In this paper, we assessed the quality of satellite-derived incoming radiative flux products (DSSF for solar irradiance and DSLF for longwave irradiance) in mountainous terrain, by conducting a thorough inter-comparison study involving kilometric resolution forecasts from the NWP system AROME and fields from the SAFRAN analysis system. A new satellite-derived product for LW↓ iradiance (DSLFnew) was developed using the DSLF algorithm fed by AROME forecasts. An eval-610 uation of all available products was performed against in situ measurements using four years of data in the French Alps and the Pyrenees. The result analysis showed that DSSF products are best for solar irradiance, despite an underestimation at the highest altitudes, while AROME is associated with a strong positive bias and SAFRAN with a negative bias. In terms of longwave irradiance, contrasted results were obtained at the mountain stations, all falling within the range of uncertainty of sensors.

615
A systematic underestimation by AROME, DSLF and DSLFnew was highlighted. The negative bias of DSLF was reduced by DSLFnew up to mid-altitudes but enhanced at high altitudes due to a too strong altitudinal gradient associated with the cold bias in AROME near-surface air temperature at high altitudes. A spatial comparison of the datasets showed that AROME and DSSF better represent the spatial variability of SW↓ fluxes in mountains by comparison with SAFRAN. These results are 620 encouraging and highlight the potential benefits of using DSSF, DSLF and DSLFnew as radiative forcing for snowpack modelling in mountainous terrain. Their relatively good quality in mountains as compared to lower altitudes also supports the use of these data as climatological inputs and/or validation datasets for NWP models over complex domains such as mountains, where incoming radiative flux measurements are scarce.

625
An evaluation of distributed snowpack simulations by Crocus driven by AROME and the different irradiance datasets was then conducted in the French Alps and the Pyrenees. We showed that replacing AROME irradiances by DSSF and DSLFnew increased the positive bias of snow depth, despite an overall better performance of these datasets in terms of incoming radiative fluxes. Therefore, an improved meteorological forcing does not ensure more accurate snowpack simulations. This is 630 mostly due to error compensations within the atmospheric forcing and the snowpack model. Complementary studies are sorely needed to identify the cause of the underestimated melting, which cannot be attributed to radiative fluxes. They should tackle factors such as the turbulent fluxes simulated by AROME-Crocus and the albedo parametrisation in Crocus. Multiphysical ensemble snowpack modelling would also enable to account for simulation errors (Lafaysse et al., 2017). Until such 635 improvements are performed in the AROME-Crocus modelling context, the LSA SAF products of incoming radiative fluxes can be used to improve understanding of snowpack models as well as in other snowpack-related studies, because they provide irradiance data of reasonable quality in mountainous areas.           Fig. 9). Middle: Cumulated melting represented with the same colours. Bottom: Mean daily irradiance differences with AROME for DSSF (green), DSLFnew (orange) and SAFRAN irradiances (red). Table 1. List of ground stations, associated mountain range, altitude of the observation, altitude of the associated LSA SAF and AROME grid points, measurement uncertainties, number of hourly SW↓ observations (N), mean observation and bias/RMSE for DSSF, AROME and SAFRAN computed when the sun is not masked, from 1 August 2010 to 31 July 2014. The best metrics are given in bold. The mountain range of each station is indicated by A (Alps), J (Jura) or P (Pyrenees). Asterisks indicate official uncertainties which seem too optimistic.