Articles | Volume 25, issue 8
Hydrol. Earth Syst. Sci., 25, 4455–4471, 2021

Special issue: Changes in the Mediterranean hydrology: observation and...

Hydrol. Earth Syst. Sci., 25, 4455–4471, 2021

Research article 19 Aug 2021

Research article | 19 Aug 2021

Snowpack dynamics in the Lebanese mountains from quasi-dynamically downscaled ERA5 reanalysis updated by assimilating remotely sensed fractional snow-covered area

Snowpack dynamics in the Lebanese mountains from quasi-dynamically downscaled ERA5 reanalysis updated by assimilating remotely sensed fractional snow-covered area
Esteban Alonso-González1, Ethan Gutmann2, Kristoffer Aalstad3, Abbas Fayad4, Marine Bouchet5, and Simon Gascoin5 Esteban Alonso-González et al.
  • 1Instituto Pirenaico de Ecología, Spanish Research Council (IPE-CSIC), Zaragoza, Spain
  • 2Research Application Laboratory, National Center for Atmospheric Research (RAL-NCAR), Boulder, CO, USA
  • 3Department of Geosciences, University of Oslo, Oslo, Norway
  • 4Centre for Hydrology, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
  • 5Centre d'Etudes Spatiales de la Biosphère (CESBIO), UPS/CNRS/IRD/INRA/CNES, Toulouse, France

Correspondence: Esteban Alonso-González (


The snowpack over the Mediterranean mountains constitutes a key water resource for the downstream populations. However, its dynamics have not been studied in detail yet in many areas, mostly because of the scarcity of snowpack observations. In this work, we present a characterization of the snowpack over the two mountain ranges of Lebanon. To obtain the necessary snowpack information, we have developed a 1 km regional-scale snow reanalysis (ICAR_assim) covering the period 2010–2017. ICAR_assim was developed by means of an ensemble-based data assimilation of Moderate Resolution Imaging Spectroradiometer (MODIS) fractional snow-covered area (fSCA) through an energy and mass snow balance model, the Flexible Snow Model (FSM2), using the particle batch smoother (PBS). The meteorological forcing data were obtained by a regional atmospheric simulation from the Intermediate Complexity Atmospheric Research model (ICAR) nested inside a coarser regional simulation from the Weather Research and Forecasting model (WRF). The boundary and initial conditions of WRF were provided by the ERA5 atmospheric reanalysis. ICAR_assim showed very good agreement with MODIS gap-filled snow products, with a spatial correlation of R=0.98 in the snow probability (P(snow)) and a temporal correlation of R=0.88 on the day of peak snow water equivalent (SWE). Similarly, ICAR_assim has shown a correlation with the seasonal mean SWE of R=0.75 compared with in situ observations from automatic weather stations (AWSs). The results highlight the high temporal variability in the snowpack in the Lebanese mountain ranges, with the differences between Mount Lebanon and the Anti-Lebanon Mountains that cannot only be explained by hypsography as the Anti-Lebanon Mountains are in the rain shadow of Mount Lebanon. The maximum fresh water stored in the snowpack is in the middle elevations, approximately between 2200 and 2500 m a.s.l. (above sea level). Thus, the resilience to further warming is low for the snow water resources of Lebanon due to the proximity of the snowpack to the zero isotherm.

1 Introduction

The hydrological processes related to mountain areas are essential for the water supplies for a large part of humanity (Viviroli et al., 2007). Despite the relatively mild temperature of the Mediterranean climates, mountains there often exhibit deep and long-lasting snowpacks accumulating more than 3 m and an average snow season of 5 months at the summit areas (Alonso-González et al., 2020b; Fayad et al., 2017b). Thus, as most of the annual precipitations falls during the winter season (García-Ruiz et al., 2011) the mountain snowpack strongly reshapes the hydrographs to sustain high flows until the end of the spring, permitting better synchronization of water demand and availability during the dry season (García-Ruiz et al., 2011). Mediterranean snowpacks are characterized by a high interannual variability, which affects the amount and seasonality of river flows (López-Moreno and García-Ruiz, 2004). Despite this variability, the thickness and high density exhibited by the snowpack in the Mediterranean climate (Fayad et al., 2017b) makes them an effective water storage system. In addition, high sublimation rates are associated with Mediterranean snowpacks (Fayad and Gascoin, 2020; Herrero et al., 2016; Schulz and de Jong, 2004). The fact that snowpack conditions are close to isothermal during most of the snow season makes them highly sensitive to the current climate warming (Alonso-González et al., 2020a; López-Moreno et al., 2017; Yilmaz et al., 2019).

The Lebanese mountains are a clear example of Mediterranean mountains where snow exerts a key control on the hydrology, and water resources are critically dependent on the interannual fluctuations of the snowpack (El-Fadel et al., 2000). Despite their importance, snow observations in the region are scarce (Fayad et al., 2017a), making the study of distributed snow dynamics challenging. Recently, Fayad and Gascoin (2020) developed distributed snowpack simulations over key areas of Mount Lebanon, forcing the model by interpolating observations of the few existing automatic weather stations (AWSs) using the SnowModel by Liston and Elder (2006). They showed the importance of the liquid water percolation scheme given the isothermal condition of the snowpack and estimated the snow water equivalent over three key catchments in the windward western divide of Mount Lebanon. However, due to the lack of meteorological data outside this area, these simulations did not cover the whole mountain area of the country and were limited to three snow seasons.

Remote sensing and numerical modeling have become reliable tools to generate useful meteorological information for mountain regions (Lundquist et al., 2019) and also to generate robust snow data worldwide. Atmospheric reanalyses are a valuable source of long-term (multidecadal) climatological information, especially at planetary scales (e.g., Wegmann et al., 2017; Wu et al., 2018). However, spatially downscaling such products is mandatory to derive relevant snow information over complex terrain (Baba et al., 2018b, and Mernild et al., 2017, among others). Dynamical downscaling has been shown to outperform statistically gridded products for meteorological variables in complex terrain (Gutmann et al., 2012b). More specifically, high-resolution fully dynamical meteorological models can reproduce the snowfall patterns over complex terrain (Ikeda et al., 2010; Rasmussen et al., 2011). However, the computational cost of full dynamical downscaling solutions becomes prohibitive for large domains at high spatial resolutions. To reduce the computational cost, different solutions of varying complexity have been developed using statistical interpolations corrected with the topography or using simplifications of the atmospheric dynamics (Fiddes and Gruber, 2014; Gutmann et al., 2016a; Liston and Elder, 2006). In this way, energy and mass balance snowpack models have been coupled with atmospheric models to develop multidecadal snow simulations (Alonso-González et al., 2018, and van Pelt et al., 2016, among others). In addition, remote sensing products have been widely used to study the duration of and variability in the snow cover (Gascoin et al., 2015; Saavedra et al., 2017; Yilmaz et al., 2019). Assimilation of remotely sensed snow cover observations has shown considerable potential for improving numerical snowpack models outputs in both distributed (e.g., Baba et al., 2018a; Margulis et al., 2016) and semi-distributed simulations (Cluzet et al., 2021; Fiddes et al., 2019). These approaches are particularly promising in data-scarce regions to reduce the biases in atmospheric forcing.

In this work, we have simulated the snowpack of the Lebanese mountains, as an alternative to sparse snowpack observations. We have generated a 1 km resolution snowpack reanalysis, using an ensemble-based assimilation of a fractional snow-covered area (fSCA) obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite sensor. More specifically, the ERA5 reanalysis (Hersbach, 2020) was dynamically downscaled using regional atmospheric models in two steps. First, a 10 km resolution atmospheric simulation, using the Weather Research and Forecast model (WRF; Skamarock et al., 2008a), was performed covering the period between 2010 and 2017. Then, a finer 1 km simulation, using the Intermediate Complexity Atmospheric Research model (ICAR; Gutmann et al., 2016a), was nested inside the previous WRF simulation covering the same time period. To improve the ICAR snowpack outputs, the simulated meteorological data were used to force an energy and mass balance snowpack model, the Flexible Snow Model (FSM2; Essery, 2015a), while perturbing the meteorological fields to generate an ensemble of snowpack simulations. Then, the particle batch smoother (PBS; Margulis et al., 2015), a Bayesian data assimilation scheme, was applied to assimilate daily remotely sensed fSCA information. We tested the generated snow products in the mountains of Lebanon with independent observations. Finally, the dynamics of the snowpack in the mountains of Lebanon are studied from the generated multi-year snow time series. The objectives of this paper are (i) to explore the potential of a methodology to develop a snowpack reanalysis over data-scarce regions and (ii) to describe the main snowpack dynamics over the Lebanese mountains. This is the first use of ICAR for generating a snow reanalysis.

2 Study area

Lebanon is a country located on the eastern Mediterranean Sea between 33 and 35 N. Its climatology is typically Mediterranean (Peel et al., 2007) and is influenced mainly by its proximity to the Mediterranean Sea and its complex topography (Fig. 1). There are two main mountain ranges that run in parallel to the Mediterranean coast from north to south. These mountain ranges are Mount Lebanon and the Anti-Lebanon Mountains, reaching 3088 m a.s.l. (above sea level; Qurnat as Sawdā peak) and 2814 m a.s.l. (Mount Hermon peak), respectively. The Lebanese mountains are highly karstified, encouraging the infiltration of rainfall and snowmelt. The land cover is mostly composed of bare rocks and soils with irregularly distributed patches of shrubland, as well as oaks and pine forest.

Despite Lebanon having more available water resources than its neighboring countries, it is considered a water-scarce region (El-Fadel et al., 2000) where droughts are frequent and are expected to increase due to climate change (Farajalla et al., 2010). The particular spatial distribution of its mountain ranges constitutes an effective topographical barrier to humidity advected from the Mediterranean Sea, enhancing the winter precipitation as a consequence of orographic effects (Jomaa et al., 2019). In these mountain ranges, the combined effects of the orography and Mediterranean climate result in yearly seasonal snowpack over a large part of the country (Mhawej et al., 2014).

It was estimated from satellite retrievals of snow cover that 31 % of the spring discharge of Lebanon is associated with snowmelt (Telesca et al., 2014). In addition, the groundwater dynamics of Lebanon are mainly controlled by the snowmelt as consequence of its karstic nature (Bakalowicz et al., 2008; El-Fadel et al., 2000). Thus, the water resource provided by the snowpack is crucial for the Lebanese society, with this need becoming more acute during the recent drought in the eastern Mediterranean (Cook et al., 2016). In addition, the water stress increased notably in recent years partially due to the increase in domestic water demand, agricultural water use and the Syrian refugee crisis (Jaafar et al., 2020) but also due to the poor management of the water resources and water pollution.

Figure 1Atmospheric models domain configuration (a) and a localization map of Lebanon (b). The red dots represent the automatic weather station (AWS) positions.

3 Data and methods

3.1 Regional atmospheric simulations configuration

To generate the meteorological forcing, we used the ICAR atmospheric model nested inside a WRF simulation forced by the ERA5 reanalysis. The WRF model was used to generate a regional atmospheric simulation on a 10 km × 10 km grid, covering the eastern part of the Mediterranean Sea, with 179 × 179 grid cells centered over Lebanon's mountains (Fig. 1). In the vertical dimension, the domain was composed of 35 levels, with the top set to 50 hPa, similar to other studies over Mediterranean regions (Arasa et al., 2016). The simulation covers the period from 1 January 2010 to 30 June 2017, using the first 9 months as spin-up period allowing for physical equilibrium between the external forcings and the land model (Montavez et al., 2017). We used the ERA5 reanalysis data set at an hourly frequency as the boundary and initial conditions of the WRF (version 3.8) model. The ERA5 data set is an atmospheric reanalysis, which replaces the widely used ERA-Interim reanalysis (Dee et al., 2011). It has a spatial resolution of 30 km, with 138 vertical levels and the top at 80 km. It has been shown to outperform ERA-Interim in many climatological applications and as a forcing data set for different modeling applications (Albergel et al., 2018; Tarek et al., 2020, and Wang et al., 2019, among others). The parameterization schemes used in the WRF simulation include the Thompson cloud microphysics scheme (Thompson et al., 2008), the NCAR Community Atmosphere Model (CAM) scheme for both shortwave and longwave radiations (Neale et al., 2004), the Noah-MP scheme for the land surface physics (Niu et al., 2011), the Mellor–Yamada–Janjic scheme for the planetary boundary layer (Janjic, 2002) and the Betts–Miller–Janjic scheme (Betts and Miller, 1986; Janjic, 1994) for deep and shallow convection. This WRF configuration has shown its consistency in previous studies simulating the seasonal snowpack over complex terrain (Ikeda et al., 2010; Rasmussen et al., 2011). In addition to the described parameterization, we applied the spectral nudging technique to satisfy the large-scale atmospheric conditions at the higher altitudes, while allowing the model to have its own dynamics inside the planetary boundary layer (Von Storch et al., 2000; Waldron et al., 1996). The spectral nudging technique was applied for the wind vectors, temperature and geopotential, with a wave number of one in each direction, based on the parameters recommended by Gómez and Miguez-Macho (2017), and nudging the waves above ∼1000 km wavelength.

Next, the ICAR model was used to obtain a finer 1 km × 1 km spatial grid atmospheric simulation nested in the aforementioned WRF simulation domain. This enabled us to significantly reduce the high computational cost compared to a long-term high-resolution WRF simulation. ICAR is a 4D meso-atmospheric model designed for downscaling purposes and is based on the linear mountain wave theory. The linear theory allows ICAR to compute the main dynamical effect of topography on the atmosphere using an analytical solution, thus avoiding the need to solve the Navier–Stokes equations and reducing the computational cost by a factor of 100. The center of the ICAR simulation was established in the center of the WRF simulation, using 179 × 179 grid cells in both latitude and longitude directions and preventing the boundaries from intersecting complex terrain. The model top was situated at 4150 m above the topography with 12 vertical levels, using the default model level heights (Horak et al., 2019). The model configuration used the Thompson cloud microphysics scheme (Thompson et al., 2008), the Noah land surface model (Chen and Dudhia, 2001) and the multidimensional positive definite advection transport algorithm (MPDATA) for the advection (Smolarkiewicz and Margolin, 1998). Convection schemes were not implemented for this simulation, and the radiative fluxes at the surface were prescribed by WRF. The lack of convection could have some impact on the total amount of precipitation and, therefore, on the seasonal snowpack. However, such deviations in the total amount of precipitation are partly compensated by the PBS (as described in Sect. 3.3.2).

3.2 Ensemble-based fractional snow cover assimilation

3.2.1 MODIS fractional snow cover area data estimation

For this study, we used satellite observations of fSCA, assimilated in an ensemble of snow simulations, to improve the snow water equivalent products (SWE) of ICAR. The daily fSCA information was obtained by means of the MODIS sensor, which is orbiting the Earth on board two satellites, Terra and Aqua. We have chosen MODIS because of its daily revisit time combined with a spatial resolution of 500 m, which is higher than our ICAR simulation. More specifically, we have used the normalized difference snow index (NDSI) retrievals of collection 6 of the NASA snow-cover products MOD10A1 (Terra) (Hall et al., 2006) and MYD10A1 (Aqua) (Hall and Riggs, 2016) distributed by the National Snow and Ice Data Center. To estimate the fSCA from the MODIS NDSI we have used a linear function following Salomonson and Appel (2004). The coefficients of the function were optimized using a series of 20 m resolution snow products from Theia Snow collection (Gascoin et al., 2019). The Theia Snow collection provides snow cover area maps derived from Sentinel-2 observations. The revisit period of Sentinel-2 is, at most, 5 d since the launch of Sentinel-2B (i.e., after March 2017). It can be shorter in areas where successive swaths overlap laterally. We downloaded 645 Theia Sentinel-2 snow products acquired between 3 September 2017 and 24 December 2018 over Lebanon. For every Sentinel-2 image we can match a MODIS image since there is a MODIS image over Lebanon for every day during the same period. Theia binary snow maps were resampled to 500 m fSCA in the same grid as the MODIS products by averaging the contributing pixels. By comparing these fSCA Theia maps with the MOD10A1 products, we could find 5.84×104 cloud-free pixels which corresponded to MOD10A1 snow-covered pixels on the same date. A subset of 40 % of the NDSI fSCA were used to fit the linear function using the least squares method. The optimized function was tested against the remaining data and yielded an fSCA RMSE of 11 % and a mean absolute error of 5.7 %. The same analysis was done with MYD10A1 (Aqua) products, but we opted not to use them in the remainder of the analysis because they exhibited a lower agreement with the Theia Sentinel-2 snow cover products (RMSE of 21 %). The lower agreement of MYD10A1 is likely due to degraded sensors (Wang et al., 2012) but may also be related to the difference between the overpass time of Sentinel-2 (10:30 local time – LT) and Aqua (13:30 LT), while Terra shares the same overpass time as Sentinel-2.

We reprojected the generated MODIS fSCA products to the spheroid datum (6370 km Earth radius) Lambert conformal projection used in the ICAR simulation. To avoid artifacts as a consequence of the data gaps in MODIS imagery caused by the cloud cover, we performed the aggregation when the majority of the MODIS cells used to calculate each new resampled cell was cloud free (less than 25 % cloud cover); otherwise, the cell was considered empty or missing for the scene in question. In previous studies, the MODIS fSCA products have shown to have a good performance in retrieving fSCA information compared with field observations – even considering its moderate resolution (Aalstad et al., 2020). Thus, they are a robust resource to use when developing regional-scale snow reanalysis.

3.2.2 Particle batch smoother implementation

The assimilation procedure was implemented using the PBS scheme (Margulis et al., 2015). The PBS assigns a weight to each ensemble member according to its agreement with the observations through Bayes' theorem. The most obvious advantage of this technique is its computational efficiency, as it avoids the resampling step common in other assimilation algorithms. A complete description of the PBS can be found in Margulis et al. (2015). It is also summarized in Aalstad et al. (2018) and Fiddes et al. (2019). The PBS has been shown to perform well relative to other assimilation algorithms when used to assimilate fSCA information (Aalstad et al., 2018; Margulis et al., 2015), even though it can suffer from particle degeneracy as consequence of a highly inhomogeneous distribution of weights (Van Leeuwen, 2009). In this context, the PBS has been successfully used to develop a series of snowpack reanalyses (Cortés et al., 2016; Fiddes et al., 2019; Margulis et al., 2016).

For the prior of the PBS implementation, we generated an ensemble of snowpack simulations, forcing the FSM2 (Essery, 2015a) with the ICAR predicted surface meteorology. The configuration of the FSM2 model includes an albedo correction, as snow ages differently in time for melting and cold snow and increases with snowfall, with a maximum of 0.9. The compaction rate was calculated based on overburden and thermal metamorphism (Verseghy, 1991). The turbulent exchange coefficient was stability corrected based on the bulk Richardson number. The thermal conductivity was calculated based on snow density. Finally, the FSM2 configuration accounted for the retention and refreezing of water inside the snowpack. Such a configuration has been shown to properly simulate the inter- and intra-annual variability in the snowpack dynamics over mountains with a similar Mediterranean climate (Alonso-González et al., 2018).

To generate the ensemble of forcing data sets, we perturbed the precipitation and the 2 m air temperature surface fields of the ICAR output using log-normal and normal (Gaussian) probability density functions, respectively. We choose the mean of the probability functions from the averaged biases of the ICAR simulation, estimated from independent observations provided by three mountain AWSs at the locations shown in Fig. 1 (Fayad et al., 2017a). The variance of the probability distribution functions was calculated by doubling the variance of the errors to increase the spread of the ensemble to cover the apparent uncertainty in the ICAR outputs. The precipitation phase had to be recalculated for the new synthetic temperatures for each ensemble member. Due to the strong dependency of the snowpack over Lebanon on precipitation phase, simple temperature-threshold-based precipitation-phase partitions are not recommended (Fayad and Gascoin, 2020). Instead, we have used the psychrometric energy balance method approach proposed by Harder and Pomeroy (2013), where the precipitation phase is estimated by means of the estimation of the temperature of the falling hydrometeor calculated from the air temperature and relative humidity. A total of 400 ensemble members per ICAR cell were independently generated by randomly drawing multiplicative time-constant parameters from the log-normal probability function for precipitation and additive parameters from the normal probability function for the 2 m air temperature.

To estimate the fSCA of each ensemble member, we used the probabilistic snow depletion curve proposed by Liston (2004). This model simulates the subgrid peak SWE distribution using a log-normal probability density function. Then, the fSCA is diagnosed using the accumulated melt depth estimated from the energy balance outputs of the FSM2, the peak mean SWE and the peak subgrid coefficient of variation (CV) in the log-normal probability density function, assuming a constant melt over the grid cell. The CV used in this model is strongly controlled by the characteristics of the terrain. We have included the CV parameter as part of the assimilation, perturbing its value inside the recommended values in Liston (2004), using a mean of 0.4 and a variance of 0.01 (Aalstad et al., 2018). The PBS was implemented over the fSCA ensemble over each grid cell and season independently, using the values of the melting season corresponding to the months of March through June. Finally, the generated SWE products (hereafter ICAR_assim) were estimated from the weighted mean of the SWE of the ensemble members, where the weights were obtained using the PBS. A schematic description of the whole process is presented in Fig. 2.

Figure 2Schematic flow chart of the ICAR_assim snow product development.


3.3 Validation procedure and analysis of the SWE products

The ICAR atmospheric simulation and the ICAR_assim products were compared against independent observations. First, the ICAR atmospheric simulation was compared with three AWSs located in the main mountain range of the domain (Fayad et al., 2017a; Fig. 1). Temperature and precipitation measurements were aggregated to the hourly model output frequency from the original 30 min time resolution. Then, the temperature and precipitation biases were estimated. The precipitation data were available only in two of the AWSs. The error values and its variance were used to define the shape of the probability density functions of the perturbation parameters described above to generate each ensemble.

Table 1AWS geographical coordinates and elevations. Elevation of the ICAR cell that contains each AWS.

Download Print Version | Download XLSX

After the PBS implementation, we compared the ICAR and ICAR_assim snow products with the snow depth observed information derived from a Campbell Scientific SR50A acoustic gauge at the three AWSs. The observed snow depth was transformed into SWE by assuming a constant snow density value of 467 kg m−3 estimated from observations in the area (Fayad et al., 2017a). This was necessary to make the AWS data comparable with the ICAR snow outputs as they are provided only as SWE. Even if it is commonly implemented in operational atmospheric forecast models, the assumption of a constant density could introduce obvious bias in the SWE estimation (Dawson et al., 2017). In the Mediterranean snowpacks, such biases are partially reduced as consequence of the high densification rates of the snowpack (Bormann et al., 2013; Fayad et al., 2017b). However, we introduced a sensitivity analysis in the comparison, varying the density value in the range of ±15 % to illustrate such uncertainty. To compensate the big shift between the ICAR and ICAR_assim resolutions (1 km × 1 km) and the point scale nature of the AWS observations, we interpolated a new SWE series from the four nearest cells of the simulations using the inverse distance method.

The spatial accuracy of the SWE products was compared to satellite observations. First, we developed a daily gap-filled snow cover time series covering the time period of the ICAR simulation from the MODIS snow cover products using the methodology proposed by Gascoin et al. (2015). Then, the products were aggregated to estimate the averaged snow presence over each cell in percentage (P(snow)). The MODIS P(snow) product was aggregated to the ICAR grid to make it comparable. Then, we calculated the P(snow) for the ICAR and ICAR_assim simulations. We chose a SWE MODIS detection threshold of 20 mm to calculate the P(snow) from the simulated SWE series inside the range recommended by Gascoin et al. (2015). All of the spatial analysis and the data assimilation were computed over the areas that had exhibited a P(snow)>5 %, which amounts to a total surface of 4412 km2.

Figure 3ERA5 (green), WRF (blue), ICAR (red) and AWS (black) daily temperature data. The box plots represent the distribution of the errors, and the gray shadows represent the data gaps in the observations.


4 Results and discussion

4.1 Atmospheric simulation results

The use of ICAR is justified as it is computationally inexpensive compared to similar WRF simulations, while retaining a physical basis to enable simulations in regions lacking observations. The speeding up factors can range from 140, in its more complex configurations (as chosen for this study), to 800, in its simpler configurations (Gutmann et al., 2016a). However, the linear theory simplification presents some limitations when predicting the motion of the atmosphere, such as interactions between waves and turbulence (Nappo, 2012) or the lack of explicit convection. Despite these limitations, ICAR has been shown to be a valuable tool for downscaling, showing a good consistency with observations (Horak et al., 2019) and with fully dynamical WRF simulations (Gutmann et al., 2016a). Figure 3 shows how the ICAR model was able to improve the 2 m air temperature data, compared with the ERA5 reanalysis (ICAR mean error =2.8C compared with 8.5 C in ERA5), showing performances comparable to the WRF coarser simulation (WRF mean error =2.3C). It was not expected that the parent WRF simulation would need to be improved with ICAR, but the increase in resolution was necessary as the snowpack simulations require higher resolutions. The coarser ERA5 resolution that smooths the terrain leads to warm biases. This is particularly evident in the Lebanese ranges where the elevation gradient ranges from 0 to 3000 m a.s.l. in approximately 25 km (Fig. 1). Despite the clear improvement in the temperature performance, the simulation is biased towards slightly higher temperatures than the AWS data. However, the main temporal patterns and the magnitude of the temperature are well represented.

Similarly, precipitation outputs of ICAR were compared with the gauges deployed in two of the AWS sites. ICAR reduces the spread of the daily precipitation errors of ERA5, as shown in Fig. 4 (standard deviation of 11.5 mm in ERA5 compared with the 8.4 mm of ICAR), even though the ERA5 errors are already surprisingly low considering the spatial resolution and the fact that precipitation is challenging to simulate by numerical models, especially over complex terrain (Legates, 2014). This validation provides a range of uncertainty estimates to help generate the probability density functions for the perturbations of the ensemble. The selected parameters to define the shape of the normal probability density function which defines the additive perturbation to the temperature were set to a mean of 3.0 C and a variance of 1.8 C. Similarly, the parameters of the log-normal probability density function used to obtain the multiplicative perturbation factors for the precipitation were a mean of 2.0 and a variance of 0.75. Even though the parameters were designed to model the uncertainty of ICAR, they are similar to comparable implementations of the PBS (Cortés et al., 2016). Through the forced increase in the variance of the probability density functions, we ensure that the ensemble of snow simulations covers the expected uncertainty space of ICAR, while the PBS has proved to be robust to progressive variations of the perturbation parameters (Cortés et al., 2016).

Figure 4ERA5 (green), WRF (blue), ICAR (red) and AWS (black) daily precipitation data. The box plots represent the distribution of the errors, and the gray shadows represent the data gaps in the observations.


4.2 Fractional snow cover assimilation

The new proposed linear relationship function to derive fSCA from NDSI has improved the MODIS fSCA products when compared with the relationship function by Salomonson and Appel (2004; Fig. S1 in the Supplement). Using the relationship function by Salomonson and Appel (2004) resulted in a larger mean absolute error (MAE; 6.2 % compared to 5.7 %) and a root mean squared error (RMSE; 12 % compared to 11 %). The equation of the linear fit is as follows:


The performance of ICAR_assim was compared against snow depth measurements at the AWS locations (Fig. 5) and MODIS gap-filled products (Figs. 6 and 7). In general, ICAR has a tendency to underestimate the SWE compared with ICAR_assim. This is likely related to the warm biases detected in the simulation, combined with the limitations in the snow model implemented in the Noah land surface model used by ICAR (Barlage et al., 2010). Thus, future versions of ICAR with better representations of the snow processes through the implementation of more complex land surface parameterizations like Noah-MP (Niu et al., 2011), as used in the parent WRF simulation, could potentially improve the accuracy of ICAR's SWE outputs (Suzuki and Zupanski, 2018). This effect could be particularly enhanced in the mild climatic conditions of Lebanon, as larger disagreements in the SWE outputs between Noah and Noah-MP occur under warm conditions (Kuribayashi et al., 2013). However, the improvement of the snow representations of ICAR is clear when compared with ERA5 reanalysis, which was not able to reproduce the snowpack at all due to its coarse resolution.

Figure 5Comparison between observed, ICAR, FSM and ICAR_assim SWE products. The green background indicates the time steps when ICAR_assim improves the performance of ICAR.


The use of FSM to generate the ensemble of simulations introduced some uncertainty in the snow simulations. Some water years showed earlier snowmelt. As the uncertainty of the snow models associated to the forcing is higher than the uncertainty associated by the use of different model parameterizations and model structures (Günther et al., 2019), we hypothesize that such differences were caused by the differences in the precipitation phase partitioning that is challenging to simulate in the areas that remain close to 0 C during the snow season (Fayad and Gascoin, 2020). The lack of spring snowfalls in some years may have deep implications in the snowpack simulation that are not limited to its effect in the mass balance and the releasing of latent heat by refreezing the liquid precipitation. It leads to lower albedos which, combined with the high shortwave radiation of Lebanon due to its latitude, causes earlier snowmelt. However, such discrepancies are greatly minimized in ICAR_assim by the assimilation of the fSCA retrievals.

The results of the validation of ICAR_assim show a good agreement with the observations. For the estimated SWE, the RMSE and the MAE relative to the AWSs were 189.2 and 104.52 mm, respectively, after removing the summer from the analyses, with a correlation coefficient (R) of 0.75 for the annual mean SWE accumulation. Even though ICAR_assim generally shows a good agreement with the observations (especially considering the scale mismatch between the stations and ICAR_assim), some clear differences were found. Figure 5 exhibits a surprisingly high difference in the magnitude of the observed SWE and the ICAR_assim output for the 2011–2012 winter season in the third AWS. However, independent observations in the area have described an exceptional snowpack in this season, with snow depths of more than 6 m and even reaching up to 10 m locally (Koeniger et al., 2017). Such disagreements between the AWS information and the independent observations can be explained by the high spatial heterogeneity of the snow depth at point scales (López-Moreno et al., 2011). This effect was studied in the Atlas Mountains, where the agreement of the snow simulations rapidly drops when using resolutions over 250 m (Baba et al., 2019). Such spatial heterogeneity has been shown to be particularly high over Mount Lebanon due to the important role of the wind redistribution as a consequence of the geomorphology (Fayad and Gascoin, 2020). For example, Fayad and Gascoin (2020) reported large differences in the AWS data from in situ measurements on 15 January 2016 when they measured snow depths up to 258 cm in the surroundings of the third AWS location (Fig. 5; bottom panel), while the AWS sensor itself detected 7.5 cm. However, the comparison between the temporal patterns of the snow cover over Lebanon from MODIS gap-filled daily products and ICAR_assim have shown good levels of agreement with a RMSE equal to 270.2 km2, a MAE equal to 124.1 km2 over a total surface of 4412 km2 (Fig. 6) and a Pearson correlation value of R=0.88 in the annual maximum of the snow cover extent (Fig. 6). The larger spatial support of the MODIS products permits a more representative and extensive validation of ICAR_assim. Thus, the good agreement between both snow cover products and the general SWE magnitudes with the AWS observations shows the temporal consistency of the ICAR_assim reanalysis.

Figure 6Daily snow cover extent comparison between MODIS gap-filled products and ICAR_assim.


The spatial patterns of ICAR_assim were also compared with the MODIS gap-filled products (Fig. 7). The spatial comparison of the P(snow) showed a very good level of agreement, demonstrating the potential of fSCA assimilation through the PBS in improving the ICAR SWE products. The comparison showed a correlation value of R=0.98, a RMSE of =3.0 % and a MAE of =2.3 %, improving the ICAR simulation that exhibited values of R=0.79, RMSE of =14.3 % and MAE of =12.3 %. There was a general tendency by ICAR_assim to slightly underestimate the P(snow) values, especially at the lower elevations. We hypothesize that this effect could be caused by the selection of a constant SWE depth to calculate the snow cover from the ICAR_assim product. Thus, the shallow snowpacks for which the SWE values are under the selected threshold are not recorded as snow presence in the ICAR_assim, even though they could potentially be detected as snow by the MODIS sensor. In addition, the MODIS snow cover products should be considered less accurate over areas of rapid melt (Gascoin et al., 2015). Such a mismatch between ICAR_assim and MODIS, combined with the fact that the 2011–2012 snow season showed persistent cloud cover related to its exceptional snowpack, could explain the biases in Fig. 6. During the 2011–2012 snow season, the gap-filling algorithm had less information to fill the MODIS snow cover time series, while the PBS had propagated the fSCA information through the whole season from the few available observations. In summary, our results have shown how ICAR_assim can accurately reproduce the interannual and intra-annual spatiotemporal patterns of the snow cover, with a SWE magnitude comparable with independent observations that agree well in terms of temporal patterns.

Figure 7Snow probability spatial comparison between observed MODIS products and ICAR_assim.

5 Snowpack dynamics over the Lebanese mountains

ICAR_assim exhibits some limitations that should be considered. First, despite the high resolution of the reanalysis, the regional nature of the simulations prevent the representation of some processes like wind or avalanche snow redistribution. In addition, there are some other sources of uncertainty involved in the development of the reanalysis, like the depletion curve, the fSCA derived from MODIS or the structural uncertainty associated with each model. However, ICAR_assim has been shown to be consistent with the limited observations, providing a valuable resource in the data-scarce context of the Lebanese mountains.

Figure 8 shows the spatial distribution of the mean peak SWE values and its temporal coefficient of variation for the 2010–2017 time period. Such values can be influenced by the fact that the study period is relatively humid compared with the previous years (Cook et al., 2016), showing slightly higher values than a long-term climatology. However, the length of the reanalyses constitutes a reasonably representative sample of the main snowpack dynamics over the region. The snowpack over Lebanon has exhibited the high temporal variability that is characteristic of the Mediterranean snowpacks (Fayad et al., 2017b), with similar values of the coefficient of variation to those observed in other Mediterranean mountain ranges (Alonso-González et al., 2020b). The maximum accumulations reach 2000 mm of SWE and are located at the higher elevations of Mount Lebanon where there is a plateau over the elevation of the winter zero isotherm (Fayad and Gascoin, 2020). The temporal coefficient of variation of the annual peak SWE follows unequal spatial patterns. It tends to exhibit higher values over the areas sheltered from direct interaction with the warm and moist Mediterranean air. In addition, it exhibits a decreasing trend with elevation (Fig. 9), as found in other Mediterranean ranges (Alonso-González et al., 2020), reaching a minimum of 15 %.

Figure 8Averaged annual peak SWE (a) and annual coefficient of variation (b).

There are clear differences between Mount Lebanon and the Anti-Lebanon Mountains that can be just partially explained by their different orography. Despite the closeness of both Mount Lebanon and the Anti-Lebanon Mountains, they exhibit different relationships between the values of the mean peak SWE (Fig. 9; top panel) and snow duration (Fig. 9; bottom panel) and with the elevation, showing that the differences are not just related to the particular orography of each range but also with its climatological characteristics. Thus, at comparable elevations, Mount Lebanon tends to show higher values of P(snow) and mean peak SWE, with lower values of the coefficient of variation, suggesting a thicker, longer-lasting and seasonally stable snowpack. The orographic precipitation caused by the uplift of the Mediterranean moisture is a major source of precipitation in the area (Jomaa et al., 2019). That is probably why the Anti-Lebanon Mountains show lower peak accumulations than Mount Lebanon, with the Anti-Lebanon Mountains in the rain shadow, leading to lower precipitation and snow accumulation. However, despite the differences in the coefficient of variation values, they tend to become similar at the higher elevations. The same coefficient of variation occurs in the elevations where the precipitation leads the snow accumulation, while they differ at the lower elevations, where the accumulation is conditioned by the temperature. This effect suggests warmer conditions on the Anti-Lebanon Mountains as a consequence of lee-side wind effects (Foëhn-type effect) and confirms the sensitivity of the snow simulation to the chosen partition phase method over Mediterranean mountains (Fayad and Gascoin, 2020).

Figure 9Relationship between annual peak SWE and elevation (a), coefficient of variation and elevation (b) and snow duration and elevation (c).


Figure 10Mean annual evolution of SWE at different elevation bands. The dark blue line represents the Anti-Lebanon Mountains, the black line the Mount Lebanon range and the red line the relative areal coverage of each elevation above 1300 m a.s.l.


Figure 11Averaged annual water stored in the snowpack at different elevation bands. The dark blue line represents the Anti-Lebanon Mountains, the black line the Mount Lebanon range and the red line the relative areal coverage of each elevation above 1300 m a.s.l.


Figure 10 shows the averaged seasonal SWE accumulation at different elevations over both Mount Lebanon and the Anti-Lebanon Mountains. Each elevation represents the aggregated pixels of the elevation, with a range of ±50 m a.s.l. For reference, they show on average a peak SWE of 306 mm at the elevation band of 2000 m a.s.l., which is comparable to those found in the Iberian Peninsula mountain ranges (Alonso-González et al., 2020b). More specifically, the peak SWE and duration values show intermediate values between the central Iberian and Pyrenees ranges at 2000 m a.s.l but with a peak SWE coefficient of variation of 53 % that is greater than the highest values of Iberia located at Sierra Nevada with 34 %. The relative area lying at each elevation compared to the total elevation over 1300 m a.s.l. is represented to highlight the importance of the hypsography from the hydrological point of view. Thus, Lebanon exhibits a deep and long-lasting snowpack, with up to 1000 mm of peak SWE on average particularly over 2500 m a.s.l., but the relative areal coverage of such elevations is very low. This suggests that the mean peak SWE series at lower elevations could hide a large variation in mass due to the wider areas at lower elevations where many different peak SWE values can coexist, as Alonso-González et al. (2020) found in the Iberian mountain ranges.

The thick snowpacks found at the higher elevations are not necessarily the biggest freshwater resources available due to the hypsometry of the mountain area. Figure 11 shows the average amount of freshwater stored in the snowpack per elevations band. It is obvious that the maximum amount of freshwater is stored between 2100 to 2500 m a.s.l., despite the fact that thicker snowpacks are at higher elevations. The cumulative water storage in the snowpack is more than double in the medium elevation zone (average maximum up to 552 Hm3 from 1300 to 2300 m a.s.l.) when compared to the higher areas (average maximum up to 189 Hm3 at 2400 m a.s.l. and onward). This is an important part of the yearly water budget, as mean annual precipitation was estimated in to be 7200 Hm3 for the period (2010–2016; Jaafar et al., 2020). Note that this in contrast to the fact that the orography of Lebanon encourages the storage of snow in the upper areas because of the existence of a high-elevation plateau (Fayad et al., 2017a; Fayad and Gascoin, 2020). These results suggest new challenges for the water management of Lebanon in the future as a consequence of warming climate. The snowpack at low-elevation areas is more sensitive to warming (Jefferson, 2011; Marty et al., 2017; Sproles et al., 2013), particularly over areas with mild winter conditions, as has been shown in other Mediterranean regions (Alonso-González et al., 2020a).

6 Conclusions

The assimilation of MODIS fSCA through the use of the PBS has proven to be a cost-effective way to use remote sensing data in snow simulations and is particularly appropriate for simulating snow in data-scarce regions. Thus, the generated SWE products show good agreement with MODIS snow cover gap-filled data, with R=0.98, RMSE equal to 3.0 % and MAE equal to 2.3 % for the spatial map of the probability of snow. The time series of snow cover showed a R=0.88, RMSE equal to 270.2 km2, and MAE equal to 124.1 km2 over a total surface of 4412 km2. The performances, in terms of SWE magnitude with the few available point-scale observations, was R=0.75, RMSE equal to 189.2 mm and MAE equal to 104.5 mm after removing the summer from the analyses.

The snowpack over Lebanon is characterized by a high temporal variability. Some differences exist between its two main mountain ranges. Mount Lebanon exhibits thicker, longer and more regular snowpacks compared to the Anti-Lebanon Mountains. Such differences cannot only be explained by the elevation difference but also reflect the drier conditions on the lee side of the Mount Lebanon range due the rain shadow effect. The hypsometry of Lebanon results in the most important snow freshwater reservoir being in the middle elevations (2200–2500 m a.s.l.). Snowpacks at these elevations close to the 0 C isotherm are highly vulnerable to climate warming. As such, our findings suggest big challenges for the future management of water resources over the Lebanon region.

Code and data availability

WRF code can be downloaded from (Skamarock et al., 2008b). ICAR code can be found at (Gutmann et al., 2016b). FSM2 is archived at (Essery, 2015b). The meteorological data can be found at (Fayad et al., 2017).


The supplement related to this article is available online at:

Author contributions

EAG led the conceptualization, developed the methodology, wrote the original draft, developed the software, curated and validated the data and visualized the project. EG and SG supervised, and EG, KA, AF and SG assisted with the methodology, supervised and reviewed and edited the paper. EG and MB developed the software, and MB curated the data and also supervised.

Competing interests

The authors declare that they have no conflict of interest.


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

Special issue statement

This article is part of the special issue “Changes in the Mediterranean hydrology: observation and modeling”. It is not associated with a conference.


Kristoffer Aalstad has been funded by the SatPerm project (grant no. 239918; Research Council of Norway), a personal overseas research grant (Research Council of Norway), and the European Space Agency Permafrost CCI ECV project (, last access: 14 July 2021) and acknowledges support from the Land–ATmosphere Interactions in Cold Environments (LATICE) strategic research area at the University of Oslo.

Financial support

This research has been supported by the Spanish Ministry of Economy and Competitiveness (grant no. BES-2015-071466), the Spanish Ministry of Economy and Competitiveness (grant no. CGL2014-52599-P10), and the Spanish Ministry of Economy and Competitiveness (grant no. CGL2017-82216-R).

The publication fee for this work has been covered by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Review statement

This paper was edited by María José Polo and reviewed by three anonymous referees.


Aalstad, K., Westermann, S., Schuler, T. V., Boike, J., and Bertino, L.: Ensemble-based assimilation of fractional snow-covered area satellite retrievals to estimate the snow distribution at Arctic sites, The Cryosphere, 12, 247–270,, 2018. 

Aalstad, K., Westermann, S., and Bertino, L.: Evaluating satellite retrieved fractional snow-covered area at a high-Arctic site using terrestrial photography, Remote Sens. Environ, 239, 111618,, 2020. 

Albergel, C., Dutra, E., Munier, S., Calvet, J.-C., Munoz-Sabater, J., de Rosnay, P., and Balsamo, G.: ERA-5 and ERA-Interim driven ISBA land surface model simulations: which one performs better?, Hydrol. Earth Syst. Sci., 22, 3515–3532,, 2018. 

Alonso-González, E., López-Moreno, J. I., Gascoin, S., García-Valdecasas Ojeda, M., Sanmiguel-Vallelado, A., Navarro-Serrano, F., Revuelto, J., Ceballos, A., Esteban-Parra, M. J., and Essery, R.: Daily gridded datasets of snow depth and snow water equivalent for the Iberian Peninsula from 1980 to 2014, Earth Syst. Sci. Data, 10, 303–315,, 2018. 

Alonso-González, E., López-Moreno, J. I., Navarro-Serrano, F., Sanmiguel-Vallelado, A., Aznárez-Balta, M., Revuelto, J., and Ceballos, A.: Snowpack sensitivity to temperature, precipitation, and solar radiation variability over an elevational gradient in the Iberian mountains, Atmos. Res., 243, 104973,, 2020a. 

Alonso-González, E., López-Moreno, J. I., Navarro-Serrano, F., Sanmiguel-Vallelado, A., Revuelto, J., Domínguez-Castro, F., and Ceballos, A.: Snow climatology for the mountains in the Iberian Peninsula using satellite imagery and simulations with dynamically downscaled reanalysis data, Int. J. Climatol., 40, 477–491,, 2020b. 

Arasa, R., Porras, I., Domingo-Dalmau, A., Picanyol, M., Codina, B., González, M. Á., and Piñón, J.: Defining a Standard Methodology to Obtain Optimum WRF Configuration for Operational Forecast: Application over the Port of Huelva (Southern Spain), Atmos. Clim. Sci., 6, 329–350,, 2016. 

Baba, M. W., Gascoin, S., and Hanich, L.: Assimilation of Sentinel-2 data into a snowpack model in the High Atlas of Morocco, Remote Sens., 10, 1982,, 2018a. 

Baba, M. W., Gascoin, S., Jarlan, L., Simonneaux, V., and Hanich, L.: Variations of the snow water equivalent in the ourika catchment (Morocco) over 2000–2018 using downscaled MERRA-2 data, Water, 10, 1120,, 2018b. 

Baba, M. W., Gascoin, S., Kinnard, C., Marchane, A., and Hanich, L.: Effect of Digital Elevation Model Resolution on the Simulation of the Snow Cover Evolution in the High Atlas, Water Resour. Res., 55, 5360–5378,, 2019. 

Bakalowicz, M., El Hakim, M., and El-Hajj, A.: Karst groundwater resources in the countries of eastern Mediterranean: The example of Lebanon, Environ. Geol., 54, 597–604,, 2008. 

Barlage, M., Chen, F., Tewari, M., Ikeda, K., Gochis, D., Dudhia, J., Rasmussen, R., Livneh, B., Ek, M., and Mitchell, K.: Noah land surface model modifications to improve snowpack prediction in the Colorado Rocky Mountains, J. Geophys. Res.-Atmos., 115,, 2010. 

Betts, A. K. and Miller, M. J.: A new convective adjustment scheme. Part II: Single column tests using GATE wave, BOMEX, ATEX and arctic air-mass data sets, Q. J. Roy. Meteor. Soc., 112, 693–709,, 1986. 

Bormann, K. J., Westra, S., Evans, J. P., and McCabe, M. F.: Spatial and temporal variability in seasonal snow density, J. Hydrol., 484, 63–73,, 2013. 

Chen, F. and Dudhia, J.: Coupling an advanced land surface-hydrology model with the Penn-State-NCAR MM5 modeling system. Part II: Preliminary model validation, Mon. Weather Rev., 129, 587–604,<0587:CAALSH>2.0.CO;2, 2001. 

Cluzet, B., Lafaysse, M., Cosme, E., Albergel, C., Meunier, L.-F., and Dumont, M.: CrocO_v1.0: a particle filter to assimilate snowpack observations in a spatialised framework, Geosci. Model Dev., 14, 1595–1614,, 2021. 

Cook, B. I., Anchukaitis, K. J., Touchan, R., Meko, D. M., and Cook, E. R.: Spatiotemporal drought variability in the mediterranean over the last 900 years, J. Geophys. Res., 121, 2060–2074,, 2016. 

Cortés, G., Girotto, M., and Margulis, S.: Snow process estimation over the extratropical Andes using a data assimilation framework integrating MERRA data and Landsat imagery, Water Resour. Res., 52, 2582–2600,, 2016. 

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597,, 2011. 

Dawson, N., Broxton, P., and Zeng, X.: A new snow density parameterization for land data initialization, J. Hydrometeorol., 18, 197–207,, 2017. 

El-Fadel, M., Zeinati, M., and Jamali, D.: Water resources in Lebanon: Characterization, water balance and constraints, Int. J. Water Resour. Dev., 16, 615–638,, 2000. 

Essery, R.: A factorial snowpack model (FSM 1.0), Geosci. Model Dev., 8, 3867–3876,, 2015a. 

Essery, R.: Flexible Snow Model – a multi-physics energy balance model of accumulation and melt of snow on the ground and in forest canopies, available at: (last access: 14 July 2021), 2015b. 

Farajalla, N. and Ziade, R.: Drought frequency under a changing climate in the Eastern Mediterranean: the Beka'a Valley, Lebanon, in: EGU General Assembly Conference Abstracts, p. 11653, available at: (last access: 14 July 2021), 2010. 

Fayad, A. and Gascoin, S.: The role of liquid water percolation representation in estimating snow water equivalent in a Mediterranean mountain region (Mount Lebanon), Hydrol. Earth Syst. Sci., 24, 1527–1542,, 2020. 

Fayad, A., Gascoin, S., Faour, G., Fanise, P., Drapeau, L., Somma, J., Fadel, A., Al Bitar, A., and Escadafal, R.: Snow observations in Mount Lebanon (2011–2016), Earth Syst. Sci. Data, 9, 573–587,, 2017a. 

Fayad, A., Gascoin, S., Faour, G., López-Moreno, J. I., Drapeau, L., Page, M. L., and Escadafal, R.: Snow hydrology in Mediterranean mountain regions: A review, J. Hydrol. 551, 374–396,, 2017b. 

Fayad, A., Gascoin, S., Faour, G., Fanise, P., and Drapeau, L.: Snow dataset for Mount-Lebanon (2011–2016),, 2017. 

Fiddes, J. and Gruber, S.: TopoSCALE v.1.0: downscaling gridded climate data in complex terrain, Geosci. Model Dev., 7, 387–405,, 2014. 

Fiddes, J., Aalstad, K., and Westermann, S.: Hyper-resolution ensemble-based snow reanalysis in mountain regions using clustering, Hydrol. Earth Syst. Sci., 23, 4717–4736,, 2019. 

García-Ruiz, J. M., López-Moreno, I. I., Vicente-Serrano, S. M., Lasanta-Martínez, T., and Beguería, S.: Mediterranean water resources in a global change scenario, Earth-Sci. Rev., 105, 121–139,, 2011. 

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

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

Gómez, B. and Miguez-Macho, G.: The impact of wave number selection and spin-up time in spectral nudging, Q. J. Roy. Meteor. Soc., 143, 1772–1786,, 2017. 

Günther, D., Marke, T., Essery, R., and Strasser, U.: Uncertainties in Snowpack Simulations – Assessing the Impact of Model Structure, Parameter Choice, and Forcing Data Error on Point-Scale Energy Balance Snow Model Performance, Water Resour. Res., 55, 2779–2800,, 2019. 

Gutmann, E. D., Rasmussen, R. M., Liu, C., Ikeda, K., Gochis, D. J., Clark, M. P., Dudhia, J., and Thompson, G.: A comparison of statistical and dynamical downscaling of winter precipitation over complex terrain, J. Climate, 25, 262–281,, 2012. 

Gutmann, E. D., Barstad, I., Clark, M. P., Arnold, J. R., and Rasmussen, R. M.: The Intermediate Complexity Atmospheric Research Model, J. Hydrometeorol., 17, 957–973,, 2016a. 

Gutmann, E. D., Barstad, I., Clark, M. P., Arnold, J. R., and Rasmussen, R. M.: The Intermediate Complexity Atmospheric Research model (ICAR), available at: (last access: 14 July 2021), 2016b. 

Hall, D. K. and Riggs, G. A.: MODIS/Aqua Snow Cover Daily L3 Global 500 m SIN Grid, Version 6, NASA Natl. Snow Ice Data Cent. Distrib. Act. Arch. Center, Boulder, Coloroda, USA,, 2016. 

Hall, D. K., Riggs, G. A., and Salomonson, V.: MODIS/Terra Snow Cover 8-day L3 Global 500 m Grid V005, Natl. Snow Ice Data Cent., Colorodo, USA, 2006. 

Harder, P. and Pomeroy, J.: Estimating precipitation phase using a psychrometric energy balance method, Hydrol. Process., 27, 1901–1914,, 2013. 

Herrero, J. and Polo, M. J.: Evaposublimation from the snow in the Mediterranean mountains of Sierra Nevada (Spain), The Cryosphere, 10, 2981–2998,, 2016. 

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

Horak, J., Hofer, M., Maussion, F., Gutmann, E., Gohm, A., and Rotach, M. W.: Assessing the added value of the Intermediate Complexity Atmospheric Research (ICAR) model for precipitation in complex topography, Hydrol. Earth Syst. Sci., 23, 2715–2734,, 2019. 

Ikeda, K., Rasmussen, R., Liu, C., Gochis, D., Yates, D., Chen, F., Tewari, M., Barlage, M., Dudhia, J., Miller, K., Arsenault, K., Grubišić, V., Thompson, G., and Guttman, E.: Simulation of seasonal snowfall over Colorado, Atmos. Res., 97, 462–477,, 2010. 

Jaafar, H., Ahmad, F., Holtmeier, L., and King-Okumu, C.: Refugees, water balance, and water stress: Lessons learned from Lebanon, Ambio, 49, 1179–1193,, 2020. 

Janjic, Z.: Nonsingular Implementation of the Mellor-Yamada Level 2.5 Scheme in the NCEP Meso model, NCEP Off. Note, 437, 61 pp., 2002. 

Janjic, Z. I.: The step-mountain eta coordinate model: further developments of the convection, viscous sublayer, and turbulence closure schemes, Mon. Weather Rev., 122, 927–945,<0927:TSMECM>2.0.CO;2, 1994. 

Jefferson, A. J.: Seasonal versus transient snow and the elevation dependence of climate sensitivity in maritime mountainous regions, Geophys. Res. Lett., 38, L16402,, 2011. 

Jomaa, I., Saab, M. T. A., Skaf, S., El Haj, N., and Massaad, R.: Variability in Spatial Distribution of Precipitation Overall Rugged Topography of Lebanon, Using TRMM Images, Atmos. Clim. Sci., 9, 369–380,, 2019. 

Koeniger, P., Margane, A., Abi-Rizk, J., and Himmelsbach, T.: Stable isotope-based mean catchment altitudes of springs in the Lebanon Mountains, Hydrol. Process., 31, 3708–3718,, 2017. 

Kuribayashi, M., Noh, N. J., Saitoh, T. M., Tamagawa, I., Wakazuki, Y., and Muraoka, H.: Comparison of snow water equivalent estimated in central Japan by high-resolution simulations using different land-surface models, Sci. Online Lett. Atmos., 9, 148–152,, 2013. 

Legates, D.: Climate models and their simulation of precipitation, Energy Environ., 25, 1163–1175,, 2014. 

Liston, G. E.: Representing subgrid snow cover heterogeneities in regional and global models, J. Climate, 17, 1381–1397,<1381:RSSCHI>2.0.CO;2, 2004. 

Liston, G. E. and Elder, K.: A distributed snow-evolution modeling system (snowmodel), J. Hydrometeorol., 7, 1259–1276,, 2006. 

López-Moreno, J. I. and García-Ruiz, J. M.: Influence of snow accumulation and snowmelt on streamflow in the central Spanish Pyrenees (Influence de l'accumulation et de la fonte de la neige sur les écoulements dans les Pyrénées centrales espagnoles), Hydrol. Sci. J., 49, 802,, 2004. 

López-Moreno, J. I., Fassnacht, S. R., Beguería, S., and Latron, J. B. P.: Variability of snow depth at the plot scale: implications for mean depth estimation and sampling strategies, The Cryosphere, 5, 617–629,, 2011. 

López-Moreno, J. I., Gascoin, S., Herrero, J., Sproles, E. A., Pons, M., Alonso-González, E., Hanich, L., Boudhar, A., Musselman, K. N., Molotch, N. P., Sickman, J., and Pomeroy, J.: Different sensitivities of snowpacks to warming in Mediterranean climate mountain areas, Environ. Res. Lett., 12, 074006,, 2017. 

Lundquist, J., Hughes, M., Gutmann, E., and Kapnick, S.: Our skill in modeling mountain rain and snow is bypassing the skill of our observational networks, B. Am. Meteorol. Soc., 100, 2473–2490,, 2019. 

Margulis, S. A., Girotto, M., Cortés, G., and Durand, M.: A particle batch smoother approach to snow water equivalent estimation, J. Hydrometeorol., 16, 1752–1772,, 2015. 

Margulis, S. A., Cortés, G., Girotto, M., and Durand, M.: A landsat-era Sierra Nevada snow reanalysis (1985–2015), J. Hydrometeorol., 17, 1203–1221,, 2016. 

Marty, C., Schlögl, S., Bavay, M., and Lehning, M.: How much can we save? Impact of different emission scenarios on future snow cover in the Alps, The Cryosphere, 11, 517–529,, 2017. 

Mernild, S. H., Liston, G. E., Hiemstra, C. A., Malmros, J. K., Yde, J. C., and McPhee, J.: The Andes Cordillera. Part I: snow distribution, properties, and trends (1979–2014), Int. J. Climatol., 37, 1680–1698,, 2017. 

Mhawej, M., Faour, G., Fayad, A., and Shaban, A.: Towards an enhanced method to map snow cover areas and derive snow-water equivalent in Lebanon, J. Hydrol., 513, 274–282,, 2014. 

Montavez, J. P., Lopez-Romero, J. M., Jerez, S., Gomez-Navarro, J. J., and Jimenez-Guerrero, P.: How much spin-up period is really necessary in regional climate simulations?, in: Geophysical Research Abstracts, EGU General Assembly, Vienna, Austria, 19, 2017–15806, 2017. 

Nappo, C. J.: The Linear Theory, 2nd edn., International Geophysics, Academic Press,, 2012. 

Neale, R. B., Chen, C., Lauritzen, P. H., Williamson, D. L., Conley, A. J., Smith, A. K., Mills, M., and Morrison, H.: Description of the NCAR Community Atmosphere Model (CAM 5.0), Ncar/Tn-464+Str 214, available at: (last access: 14 July 2021), 2004. 

Niu, G. Y., Yang, Z. L., Mitchell, K. E., Chen, F., Ek, M. B., Barlage, M., Kumar, A., Manning, K., Niyogi, D., Rosero, E., Tewari, M., and Xia, Y.: The community Noah land surface model with multiparameterization options (Noah-MP): 1. Model description and evaluation with local-scale measurements, J. Geophys. Res.-Atmos., 116, D12109,, 2011. 

Peel, M. C., Finlayson, B. L., and McMahon, T. A.: Updated world map of the Köppen-Geiger climate classification, Hydrol. Earth Syst. Sci., 11, 1633–1644,, 2007. 

Rasmussen, R., Liu, C., Ikeda, K., Gochis, D., Yates, D., Chen, F., Tewari, M., Barlage, M., Dudhia, J., Yu, W., Miller, K., Arsenault, K., Grubišić, V., Thompson, G., and Gutmann, E.: High-resolution coupled climate runoff simulations of seasonal snowfall over Colorado: A process study of current and warmer climate, J. Climate, 24, 3015–3048,, 2011. 

Saavedra, F. A., Kampf, S. K., Fassnacht, S. R., and Sibold, J. S.: A snow climatology of the Andes Mountains from MODIS snow cover data, Int. J. Climatol., 37, 1526–1539, 2017. 

Salomonson, V. V. and Appel, I.: Estimating fractional snow cover from MODIS using the normalized difference snow index, Remote Sens. Environ., 89, 351–360,, 2004. 

Schulz, O. and de Jong, C.: Snowmelt and sublimation: field experiments and modelling in the High Atlas Mountains of Morocco, Hydrol. Earth Syst. Sci., 8, 1076–1089,, 2004. 

Skamarock, W. C., Klemp, J. B., Dudhia, J. B., Gill, D. O., Barker, D. M., Duda, M. G., Huang, X.-Y., Wang, W., and Powers, J. G.: A description of the Advanced Research WRF Version 3, NCAR Technical Note TN-475+STR, Tech. Rep., 113,, 2008a. 

Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Barker, D. M., Duda, M. G., Huang, X.-Y., Wang, W., and Powers, J. G.: The official repository for the Weather Research and Forecasting (WRF) model, available at: (last access: 14 July 2021), 2008b. 

Smolarkiewicz, P. K. and Margolin, L. G.: MPDATA: A Finite-Difference Solver for Geophysical Flows, J. Comput. Phys., 140, 459–480,, 1998. 

Sproles, E. A., Nolin, A. W., Rittger, K., and Painter, T. H.: Climate change impacts on maritime mountain snowpack in the Oregon Cascades, Hydrol. Earth Syst. Sci., 17, 2581–2597,, 2013. 

Suzuki, K. and Zupanski, M.: Uncertainty in solid precipitation and snow depth prediction for Siberia using the Noah and Noah-MP land surface models, Front. Earth Sci., 12, 672–682,, 2018. 

Tarek, M., Brissette, F. P., and Arsenault, R.: Evaluation of the ERA5 reanalysis as a potential reference dataset for hydrological modelling over North America, Hydrol. Earth Syst. Sci., 24, 2527–2544,, 2020. 

Telesca, L., Shaban, A., Gascoin, S., Darwich, T., Drapeau, L., Hage, M. El, and Faour, G.: Characterization of the time dynamics of monthly satellite snow cover data on Mountain Chains in Lebanon, J. Hydrol., 519, 3214–3222,, 2014. 

Thompson, G., Field, P. R., Rasmussen, R. M., and Hall, W. D.: Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part II: Implementation of a new snow parameterization, Mon. Weather Rev., 136, 5095–5115,, 2008. 

Van Leeuwen, P. J.: Particle filtering in geophysical systems, Mon. Weather Rev., 137, 4089–4114,, 2009. 

van Pelt, W. J. J., Kohler, J., Liston, G. E., Hagen, J. O., Luks, B., Reijmer, C. H., and Pohjola, V. A.: Multidecadal climate and seasonal snow conditions in Svalbard, J. Geophys. Res.-Earth, 121, 2100–2117,, 2016. 

Verseghy, D. L.: Class – A Canadian land surface scheme for GCMS. I. Soil model, Int. J. Climatol., 11, 111–133,, 1991. 

Viviroli, D., Dürr, H. H., Messerli, B., Meybeck, M., and Weingartner, R.: Mountains of the world, water towers for humanity: Typology, mapping, and global significance, Water Resour. Res., 43,, 2007. 

Von Storch, H., Langenberg, H., and Feser, F.: A spectral nudging technique for dynamical downscaling purposes, Mon. Weather Rev., 128, 3664–3673,<3664:ASNTFD>2.0.CO;2, 2000. 

Waldron, K. M., Paegle, J., and Horel, J. D.: Sensitivity of a spectrally filtered and nudged limited-area model to outer model options, Mon. Weather Rev., 124, 529–547,<0529:SOASFA>2.0.CO;2, 1996. 

Wang, C., Graham, R. M., Wang, K., Gerland, S., and Granskog, M. A.: Comparison of ERA5 and ERA-Interim near-surface air temperature, snowfall and precipitation over Arctic sea ice: effects on sea ice thermodynamics and evolution, The Cryosphere, 13, 1661–1679,, 2019.  

Wang, D., Morton, D., Masek, J., Wu, A., Nagol, J., Xiong, X., Levy, R., Vermote, E., and Wolfe, R.: Impact of sensor degradation on the MODIS NDVI time series, Remote Sens. Environ., 119, 55–61,, 2012. 

Wegmann, M., Orsolini, Y., Dutra, E., Bulygina, O., Sterin, A., and Brönnimann, S.: Eurasian snow depth in long-term climate reanalyses, The Cryosphere, 11, 923–935,, 2017. 

Wu, X., Che, T., Li, X., Wang, N., and Yang, X.: Slower Snowmelt in Spring Along With Climate Warming Across the Northern Hemisphere, Geophys. Res. Lett., 45, 12331–12339,, 2018. 

Yilmaz, Y., Aalstad, K., and Sen, O.: Multiple Remotely Sensed Lines of Evidence for a Depleting Seasonal Snowpack in the Near East, Remote Sens., 11, 483,, 2019. 

Short summary
Snow water resources represent a key hydrological resource for the Mediterranean regions, where most of the precipitation falls during the winter months. This is the case for Lebanon, where snowpack represents 31 % of the spring flow. We have used models to generate snow information corrected by means of remote sensing snow cover retrievals. Our results highlight the high temporal variability in the snowpack in Lebanon and its sensitivity to further warming caused by its hypsography.