The role of liquid water percolation representation to estimate snow water equivalent in a Mediterranean mountain region (Mount Lebanon)

In many Mediterranean mountain regions, the seasonal snowpack is an essential yet poorly known water resource. Here, we examine, for the first time, the spatial distribution and evolution of the snow water equivalent (SWE) during three snow seasons (2013-2016) in the coastal mountains of Lebanon. We run SnowModel (Liston and Elder, 2006a), a spatially10 distributed, process-based snow model, at 100 m resolution forced by new automatic weather station (AWS) data in three snow-dominated basins of Mount Lebanon. We evaluate a recent upgrade of the liquid water percolation scheme in SnowModel, which was introduced to improve the simulation of the snow water equivalent (SWE) and runoff in warm maritime regions. The model is evaluated against continuous snow depth and snow albedo observations at the AWS, manual SWE measurements, and MODIS snow cover area between 1200 m and 3000 m a.s.l.. The results show that the new percolation 15 scheme yields better performance especially in terms of SWE but also in snow depth and snow cover area. Over the simulation period between 2013 and 2016, the maximum snow mass was reached between December and March. Peak mean SWE (above 1200 m a.s.l.) changed significantly from year to year in the three study catchments with values ranging between 73 mm and 286 mm we (RMSE between 160 and 260 mm w.e.). We suggest that the major sources of uncertainty in simulating the SWE, in this warm Mediterranean climate, can be attributed to forcing error but also to our limited understanding of the separation 20 between rain and snow at lower-elevations, the transient snow melt events during the accumulation season, and the highvariability of snow depth patterns at the sub-pixel scale due to the wind-driven blown-snow redistribution into karstic features and sinkholes. Yet, the use of a process-based snow model with minimal requirements for parameter estimation provides a basis to simulate snow mass SWE in non-monitored catchments and characterize the contribution of snowmelt to the karstic groundwater recharge in Lebanon. While this research focused on three basins in the Mount Lebanon, it serves as a case study 25 to highlight the importance of wet snow processes to estimate SWE in Mediterranean mountain regions.


Introduction
Winter snowpack has a large influence over the hydrological processes in Mediterranean mountains and the spring meltwater from the seasonal snowpack often represent key water resource for downstream areas (Fayad et al., 2017a). In such regions, 2 large fraction of the annual precipitation falls during the winter season, which extends over a period of five months (Fayad et 30 al., 2017a). With a limited precipitation season and warm and dry summers, the Mediterranean mountainous areas are among the most vulnerable to a warming climate (e.g, Nohara et al., 2006;Loarie et al., 2009;Kumar et al., 2016).
Our understanding of the snowpack hydrological contribution in many Mediterranean mountains remains incomplete. Snow monitoring in regions such as the central Chilean Andes, Moroccan Atlas, and the Eastern Mediterranean Sea is relatively new (Fayad et al., 2017a). The increased interest, over the past decade, in understanding snow accumulation and melt in the less 35 monitored Mediterranean mountains is mainly driven by the need to meet the increased demand for water resources.
Mediterranean mountain regions, which share many characteristics of the snow in maritime regions (Fayad et al., 2017a), are characterized by their deep multilayer winter snowpack, which can exceed the 3 meters in depth, coarser grain size due to wetting, almost isothermal winter snowpack, and the high sensitivity of the precipitation phase to air temperature (Sturm et al., 1995;Sproles, et al., 2013;Pflug et al., 2019). Capturing the meltwater runoff dynamics in Mediterranean can be 40 challenging especially when transient snowmelt events are common during the winter season.
To reduce the uncertainty in simulating the snow water equivalent (SWE) in warm maritime regions, Pflug et al. (2019) implemented a discrete representation of snow liquid water percolation into the snowpack using gravity drainage to replace the standard snow density threshold routine in the distributed snowpack model SnowModel (Liston and Elder, 2006a). The new scheme by Pflug et al. (2019) showed improvement in simulating SWE compared to the default model (Liston and Elder,45 2006a) which tended to strongly overestimate the SWE at a Snow Telemetry station (SNOTEL) in the Washington Olympic Mountains, USA. The default liquid water percolation scheme did not allow the release of meltwater runoff during the accumulation season, causing a large over-accumulation of SWE. This process is typical of warm Mediterranean snowpack, such as the Mount Lebanon snowpack.
In this study, we make use of a newly available snow dataset (Fayad et al., 2017b) to evaluate this new percolation scheme 50 in SnowModel to compute the SWE distribution and its evolution, for the first time, in three basins in the warm Mediterranean mountains of Lebanon. We evaluate the model outputs using both in situ (continuous snow depth records at AWS and snow course measurements) and remote sensing data (MODIS snow cover area). We show that the improved percolation scheme led to better performance therefore, we used this model version to estimate the daily evolution of the SWE distribution over three water years (2013-2014, 2014-2015 and 2015-2016, hereafter W1314, W1415, and W1516, respectively). 55

Study area
Lebanon has relatively abundant water resources (857 m 3 per capita in 2014 according to World Bank data). Nevertheless, water shortages are frequent during the summer. The water scarcity became particularly acute in Lebanon between 1998 and 2012, when the worst drought of the past nine centuries hit the Eastern Mediterranean (Cook et al., 2016). The difficulty to meet water demand in Lebanon is also exacerbated by the poor management of the water resources (Lebanese Ministry of 60 3 Energy and Water, 2012), the increase in urbanization, and population growth over the past decades (Lebanese Ministry of Environment, 2011).
In Lebanon, the water demand for agriculture is concentrated in the irrigated Bekaa Valley between the Mount Lebanon and the Anti Lebanon mountain ranges. On the other hand, the water demand in the coastal strip of Lebanon, between the Mediterranean Sea and Mount Lebanon, is mostly associated to domestic and industrial uses. This narrow coastal plain is 65 heavily urbanized and includes the Beirut metropolitan area where approximately a third of the Lebanese population lives. In this area, the urban water supply relies on the discharge of springs and water wells from the highly karstified Jurassic limestone aquifers (Margane et al., 2013). The recharge area of this karstic system is located on Mount Lebanon, a mountain range that extends over a distance of 150 km parallel to the Mediterranean coast with a main plateau averaging 2500 m a.s.l. elevation and peak at 3088 m a.s.l. 70 Most of the precipitation occurs between January and March due to the influence of the Mediterranean climate and therefore areas above 1200 m a.s.l. receive between 50 to 67% of their total annual precipitation as snow (Shaban et al., 2004;Aouad-Rizk et al., 2005;Mhawej et al., 2014;Telesca et al., 2014). Winter precipitation events are due to the orographic lifting of the cold northerly currents from the Anatolian plateau (Fish, 1944). Snowmelt is a major contributor to the groundwater recharge (Margane et al., 2013;UNDP, 2014). In this upper cretaceous formation, usually above 1600 m a.s.l., snowmelt contributes to 75 roughly 75% of the groundwater recharge ( Fig. 1) (Margane et al., 2013;Königer and Margane, 2014).
This study focuses on the upper catchment areas of the Abou Ali, Ibrahim, and El Kelb river basins ( Fig. 1; Table 1). These basins provide key water resources to two major Governoratesthe northern districts of Mount Lebanon Governorate and the districts of the North Governorateas well as the densely populated coastal cities of Lebanon, including Beirut the capital.
The Abou Ali River provides water to the districts of Bcharre, Zgharta, Minieh, and Tripoli, the second-largest city in the 80 country. The Ibrahim River emerges from the Afqa grotto and supplies water for the districts of Kesrouan and Jbeil, a touristic area on the coast including the Byblos archaeological site. The El Kelb River basin supply water for the El Metn and Keserouan districts. The Jeita spring in the lower area of the El Kelb river basin provides 75% of the domestic water supply for the Beirut district and its northern suburbs (about 1.5 million inhabitants) (Doummar et al., 2014). The contribution of the snowfed upper cretaceous plateau to the discharge at the Jeita spring (60 m a.s.l.) was estimated at 39% (Margane et al., 2013). 85 All three basins have similar topographic, physiographic and geological characteristics. They stretch from the Mediterranean Sea to the Mount Lebanon with a mean west-facing slope. The recharge area on Mount Lebanon above 1200 m a.s.l. can be classified as mid-to high-altitude ranges with rugged topography and high plateaus (Viviroli et al., 2007). The maximum elevation exceeds 2500 m a.s.l. in the three basins. The Abou Ali river basin includes Qurnat El Sawda (3088 m a.s.l.) the highest point in Lebanon and the Levant. The land cover is mainly bare rocks and soils and short scrublands. 90 The surface drainage network is poorly developed in the upper areas of the basins as most rainfall and snowmelt-runoff infiltrates in the highly karstified limestone aquifers and aquitards beneath the surface of Mount Lebanon (Margane et al., 2013). Basin outlines used in this study ( Fig. 1; Table 1) were derived from surface topography and may differ from the actual catchments boundaries which are mostly controlled by the karst aquifer geometry and connectivity (Margane et al., 2013). The 4 topographic catchments are shown here to give an indication of the basin size and location because the hydrogeological 95 catchments are poorly known. Meteorological data at the Beirut weather station (20 m a.s.l.) operated by the national civil aviation service for the time period between 1984 and 2016 suggest that W1314 was a very dry year (fifth driest year on record between 1984 and 2016). W1415 and W1516 were less dry but had below average precipitation (Table 2). Table 1 Main characteristics of the river basins in this study (defined using the surface topography). Annual flow based on unpublished data 100 from the Litani river authority gauging stations   . Observed precipitation anomalies in Beirut (20 m a.s.l) are based on AWS meteorological data provided by the national civil aviation service.

110
We used SnowModel (Liston and Elder, 2006a), to simulate the spatial distribution of snow depth, SWE, and snowmelt over the study area. SnowModel is a physically-based snow evolution model which combines four submodels: (1) MicroMet (Liston and Elder, 2006b) spatially interpolates the meteorological forcing over the simulation domain from AWS observations and the digital elevation model (DEM) using the Barnes scheme; (2) EnBal, simulates the snowpack energy balance (Liston, 1995;Liston et al., 1999); (3) SnowPack simulates the snow depth evolution and runoff (Liston and Elder, 2006b); and (4) SnowTran-115 3D simulates blowing snow sublimation and wind-driven redistribution processes (Liston et al., 2007). Two additional submodels are available but were not used in this study: SnowAssim (Liston and Hiemstra, 2008) which allows the assimilation of SWE data and HydroFlow (Liston and Mernild, 2012), a runoff routing model. Pflug et al. (2019) implemented a gravitydriven scheme to compute the liquid water percolation through the snowpack instead of the default density threshold scheme in SnowPack submodel. In the next sections we used the default model and this new version referred to as Pflug et al. (2019). 120 6

Model setup
The model was run at half-hourly time step over a simulation domain that encompasses the three study basins (Fig. 1). The model grid resolution was set to 100 m, which is deemed sufficient to represent the main effects of the terrain on the energy and mass balance (Baba et al., 2019). The input DEM was generated by cubic resampling of the ASTER GDEM V2 in the WGS-84 UTM 36N reference system. The land cover grid was obtained from the land cover and land-use map of Lebanon 125 produced at a spatial resolution of 25 m (National Council for Scientific Research, 2015). Land classes were aggregated into 9 classes after Liston and Elder (2006a) and resampled to 100 m spatial resolution using a majority filter: Bare soils and rocks (59%), shrubland (12%), sea (8%), crops (7%), grassland (4%), clear-cut conifer (4%), urban (3%), coniferous, deciduous and mixed forests (<3%), and water (<1%). The different land classes are needed to account for the snow holding capacities and canopy radiation effects (Liston and Elder, 2006a), however bare soil is the largely dominant class in the snow dominated 130 region of the study area. First, the snow-rain threshold was set to 0°C instead of 2°C based on the analysis of the local meteorological record and other studies in similar climatic context (Harpold et al., 2017). We did not use a more complex precipitation-phase partitioning method such as those described in Harder and Pomeroy (2014) and Harpold et al. (2017) because the parameterization of such 140 methods would require introducing more parameters that we will not be able to determine for this study area. Second, we did not activate the default correction of the precipitation rate with elevation. It means that we ran the model with the minimal hypothesis that the precipitation spatial distribution is only controlled by the distances with the AWS according to the Barnes interpolation scheme (Liston et al., 2006b). This approximation is justified by the good horizontal (

Model evaluation 155
The validation dataset includes (1) half-hourly snow height (HS), snow albedo and incoming shortwave radiation collected at the AWS from W1314 to W1516, (2) bi-weekly manual SWE and snow density measurements collected near the AWS during the snow seasons of W1415 and W1516, and (3) daily snow cover area (SCA) observations from MODIS (time period from W1314 to W1516). All these data are fully described in Fayad et al. (2017b) and are available as open data in a public repository (Fayad et al., 2017c). 160 Nearly continuous snow depth is available for water years W1415 and W1516 at LAQ and for water years W1314, W1415 and W1516 at MZA and CED. Half-hourly snow depth records were averaged to the daily timestep. We computed the daily albedo of the surface below AWS from half-hourly upward and downward-looking pyranometers measurements following Stroeve et al. (2013), i.e. by computing the ratio of the daily incoming and reflected shortwave radiation totals. Summer albedo values ranged between 0.2 and 0.4 depending on the site and the year; hence, we removed albedo values lower than 0.5 to keep 165 only snow albedo measurements. The calculated daily snow albedo values and the incoming shortwave radiation were averaged to monthly values. The MODIS dataset is a daily cloud-free time series of snow cover maps providing the snow presence and absence at 500 m resolution (binary SCA). The model outputs at 100 m resolution were converted to binary SCA using a threshold of 10 mm w.e. of SWE.
We used the snow course locations to approximate the pixels that are representative for each snow course. The simulated 170 SWE were evaluated by comparing the observed SWE (revisit time ~ 11 days on average for all snow courses according to Fayad et al. (2017b)) with the simulated SWE value corresponding to the same dates when the in situ observations were taken.
Five snow courses located near the three stations are used to showcase the evolution of SWE across different elevation bands (snow courses near the three AWS have elevations at 2823 and 2834 m a.s.l. for CED, 2297 and 2301 m a.s.l. for MZA, and 1843 m a.s.l. for LAQ). 175 To further investigate the effect of the liquid water percolation scheme on the internal properties of the snowpack, we used the wet snow detection from Sentinel-1 C-band Synthetic Aperture Radar data. We use these data as a proxy of melting conditions that we compared to the simulated release of liquid water from the snowpack. The wet snow detection was done using Sentinel-1 observations over the period October 2014 to June 2016. We extracted Sentinel-1 backscatter in VV polarization mode from the Sentinel-1 Ground Range Detected (GRD) collection in Google Earth Engine (Gorelick et al., 180 2017). We spatially averaged the backscatter at each station using a buffer with a radius of 100 m. Then, we defined a reference "dry" surface backscatter for each station using the 10th percentile of the backscatter time series (CED-8.4 dB, LAQ: -10.7 dB, MZA: -6.10 dB). A negative departure of 4 dB to this reference was used to determine the occurrences of wet snow (Nagler et al., 2016).
The main reason for using daily HS, SCA, and Sentinel-1 wet snow observations in the validation of the SWE is the scarcity 185 of the spatially and temporally distributed SWE observations. However, HS is known to be a good proxy of SWE in many regions (Sturm et al., 2010). Furthermore, Fayad et al. (2017b) demonstrated that the estimation of SWE from HS could be achieved with acceptable accuracy using two years of observations in Mount Lebanon. Both the SWE and HS measurements are collected at point locations (Fayad et al., 2017c) and may not be representative at the catchment scale due to the high spatial variability of the snowpack in Mediterranean Mountains (e.g. López-Moreno et al., 2011. The use of SCA is justified 190 in the absence of distributed SWE or HS measurements to validate the model (e.g. Gascoin et al., 2015;.

SWE estimation
We used the model outputs to estimate the evolution of the SWE over the three basins (Fig. 1). The basins surface boundaries were used to estimate the SWE contribution for each basin by computing the mean daily SWE for all pixels within each basin.
To account for the hydrological contribution of SWE from the snow dominated regions, the daily distributed SWE was spatially 195 integrated for elevations above 1200 m a.s.l.. From the temporal evolution of these basin-scale SWE time series we derived the following key indicators: date and value of the peak SWE (maximum SWE during a water year), snow melt-out date (the first day of the calendar year on which SWE gets below 1 mm w.e.), distributed SWE for the peak SWE data (this map is important from a hydrological perspective as it shows the amount of snow available for melt at the end of the snowfall season for each hydrological year).  (2019) model seems to be more consistent at these stations. The model overestimates snow depth at MZA during W1516 but 215 10 the snow depth measurements are very noisy during this period suggesting that the data may be affected by a sensor issue, especially given the SWE observations presented above. Otherwise, the same positive bias as noted above (Fig. 3) can be observed at LAQ for W1415.

Figure 4
Time series of observed and modelled daily snow depth at each AWS.

220
The comparison of the modelled and observed snow cover shows that both models perform well in reproducing the snow cover evolution at the catchment scale, although the default model tend to overestimate the snow cover area during the ablation periods. This is particularly evident in spring 2016, but also during the atypical melt event in January 2014. We also evaluate the model using the snow cover duration from MODIS over the entire study area (Fig. 6). The results show that the spatial patterns of snow cover duration are generally better reproduced with the Pflug et al. (2019) scheme. The largest bias is found in the northernmost region, which corresponds to the region with the highest elevation (2700 to 3000 m a.s.l.). 230 However, the average bias for this model configuration over the entire domain remains low (2 days).
12     Table 3 shows that the correlation and RMSE are generally improved in comparison with the default model. The notable exceptions are small increases in the RMSE of snow depth and SWE at LAQ, which are hardly conclusive given the large accumulation biases at this AWS ( Fig. 3 and Fig. 4) We note that we have also evaluated the modelled shortwave radiation and snow albedo (see Figure S1 in the Supplement).
The results show that the monthly incoming shortwave is well simulated with correlation coefficients ranging between 0.94 and 0.99 and RMSE ranging between 13 W m -² (LAQ) and 32 W m -² (MZA). The snow albedo variability is less accurately captured but the RMSE remains acceptable with 0.086, 0.074, and 0.096 at CED, LAQ, and MZA respectively (Supplement 280 Fig. S1).

Estimation of SWE
Given the better performance of the Pflug et al. (2019) model in this study area, we use this model to characterize the SWE distribution and evolution over the study area. The evolution of the spatially-averaged modelled SWE over the study basins above 1200 m a.s.l. is shown in Figure 9 and the main indicators are provided in Table 4. Figure 9 shows that there is a large 285 variability in the SWE evolution during the study period hence the peak SWE is not clearly defined. However, it was consistently reached on the same date in each catchment for a given water year. During W1314 it occurred early in the snow  (Table 4). This is also illustrated in Figure 10, which presents the temporal evolution of the total mass of snow (expressed in cubic meters of water equivalent) for each catchment. The Abou Ali has the largest snow mass throughout the snow season due to its larger area and higher elevation. The temporal evolution of the three catchments was 295 similar in W1415 and W1516, however this figure also reveals that the behaviour of Abou Ali was different during W1314, with a much larger snow mass in Abou Ali in spring than in El Kelb and Ibrahim catchments. The months with the largest mean snow mass in El Kelb and Ibrahim catchments were December in W1314, and February in W1415 and W1516, while in Abou Ali it was February in W1314, March in W1415 and February in W1516. 300 Table 4 Key indicators from the spatially-averaged SWE above 1200 m a.s.l. (also see Figure 11): date and value of the peak SWE and meltout date. Here the Pflug et al. (2019)  Based on Figure 10, we selected the 01 March as a common date to represent the spatial distribution of SWE across the study domain for the three study years (Fig. 11). These maps further illustrate how the SWE in Mount Lebanon can be strongly 305 17 impacted by the inter-annual variability. On 01 March 2016, SWE exceeded 500 mm w.e. over a large portion of the Mount Lebanon range, while it was below 200 mm w.e. on 01 March 2014. Previous studies showed that SnowModel could overestimate SWE of warm maritime snowpack due to its liquid water percolation parameterization (Sproles et al., 2013;Pflug et al., 2019). We also find that the default SnowModel version 320 overestimates SWE at MZA and CED stations. The available observations show that new percolation scheme by Pflug et al. (2019) allows a clear improvement in modelled SWE but also in the snow depth and snow cover area at the catchment scale.
The improvement with the Pflug et al. (2019) model is similar to what has been reported by the same authors in the Olympic Mountains. The gravity drainage scheme allows a better representation of the liquid water percolation during the melt season, whereas the default density threshold method retains too much liquid water when applied to warm maritime snowpacks. Only 325 at LAQ the Pflug et al. (2019) version does not improve the SWE simulation. However, the snowpack is shallower and the SWE error is already large during the accumulation period at this AWS hence it is difficult to consider this AWS as good reference for evaluating a liquid water percolation scheme. This error in SWE at LAQ is probably due to errors in the precipitation input data. This issue is being investigated but so far, we do not have a clear answer. In addition, the model tends to overestimate snow depth whereas there is no clear positive or negative bias in SWE. This suggests that the modelled density 330 might be too low. High densification rates are typical of Mediterranean mountains (Fayad et al., 2017a) and this process might be insufficiently represented in SnowModel. Despite these limitations, however, the overall model evaluation with in situ and remote sensing data suggest that SnowModel modified by Pflug et al. (2019) generates a reasonable estimation of the SWE in Mount Lebanon. It is not possible to make a direct comparison of the results of this study with previous snow studies since the periods of computation do not match, but we note that the orders of magnitude of the peak SWE are compatible (Sect. 1). 335 The agreement between the measured point scale and modelled SWE in the three basins indicated good agreement, with a correlation coefficient ranging between 0.56 and 0.93 and an RMSE ranging between 16 and 26 cm w.e. (Table 3). The range of the reported RMSE values in Table 3 are not uncommon in Mediterranean like climates (Guan et al., 2013;Musselman et al., 2017). These results indicate that SnowModel modified by Pflug et al. (2019) was able to explain a large part (56-93%) of the variability in the SWE with an accuracy between 16 and 26 cm w.e.. 340 Given the fact that the three AWS are well distributed in the domain (Fig. 1), the reported RMSE values for SD, SWE, and SCA (Table 3) can be used as a possible estimate of the uncertainty of SWE in every grid cell. This uncertainty would be a conservative estimate since the random part of the error in the distributed SWE is reduced by the spatial averaging. In addition, the MODIS data showed no evidence of a large systematic SWE error over the domain except in the highest elevation region of the Abou Ali catchment, where the model tends to overestimate the snow cover duration (Fig. 6). 345 The model results tend to suggest that the spatial variability of SWE is rather low over Mount Lebanon, but this is due to the large-scale morphology of Mount Lebanon, with a strong topographic gradient from the coast to 2000 m a.s.l. followed by a high elevation plateau. This plateau is situated above the zero degree isotherm in winter therefore the snowfall regime is rather homogeneous over this area. However, Mount Lebanon is also characterized by a peculiar heterogeneous topography, visible at smaller scales (below 100 m), with many karstic surface features like sinkholes, that create highly variable snow 350 20 depth patterns due to wind-driven snow redistribution processes ( Fig. 12; Chakra et al. 2019). This variability is not resolved with this current model setup with a grid size of 100 m. This can be explained by the fact that SnowModel does not account for subgrid variability and assumes spatial uniformity over the entire grid cell (Liston and Elder, 2006a). Not accounting for snow subgrid variability is considered a fundamental limitation in the application of snow models and have implications in solving for snow ablation (DeBeer and Pomeroy, 2017). Furthermore, the presence of heterogeneous fractional snow cover 355 ( Fig. 12) have major implications on the surface energy fluxes (DeBeer and Pomeroy, 2009). If fact, the heterogeneity in snow cover fraction is known to have implications on the small-scale advection of sensible heat and other local variations in the surface energy balance terms (DeBeer and Pomeroy, 2017). This interaction between the fractional snow cover and the overlaying boundary layer is usually simplified and not captured by most models (DeBeer and Pomeroy, 2009;Mott et al., 2018). 360 Figures 5 and 9 show multiple transient melt periods during the entire snow season. These events can be clearly observed in the three basins during the months of January and February (Fig. 9), which is considered as the accumulation season.
Precipitation events during the same months did restore most of the snowpack (Fig. 9). These events, which are more common in warmer climate, such as the Mediterranean climate of Mount Lebanon, make it more difficult to distinctly separate between the accumulation and ablation seasons (DeBeer and Pomeroy, 2017). These transient snowmelt and accumulation events pose 365 major challenges in terms of simulating snow accumulation, its redistribution, and ablation due to the complex nature of these processes, which remain not fully understood (Liston, 1995;Essery et al., 2006).
Model results indicate that the largest snow reserves are found in the Abou Ali catchment (Fig. 10). This is due to its larger area but also to a larger portion of the basin is located at higher elevation bands. As a result, the SWE evolution in the Abou Ali catchment differs from the El Kelb and Ibrahim catchment. The snow season is longer and the peak SWE occurs about a 370 month later. Although the study period only spans three snow seasons, a large inter-annual variability can already be observed.
Both model simulations showed higher biases in computing SWE at lower elevation, where the default SnowModel tended to overestimate and the Pflug et al. (2019) model tended to underestimate the SWE in LAQ (Fig. 3). This can be attributed to the separation of rain from snow at lower elevations, and the fact that these regions are more sensitive to rain on snow events.
Finally, the Pflug et al. (2019) snowpack liquid water percolation algorithm may be overestimating melt at lower elevations. 375 To which extent these processes are influencing the model performance is still unknown to us.

Conclusions
The snow cover is a critical yet poorly known water resource in many Mediterranean regions (Fayad et al., 2017a). A good knowledge of the snow water equivalent (SWE) and snowmelt in Mount Lebanon is key for the evaluation and management of the water resources in Lebanon. This study presents the first simulation of the daily SWE in Mount Lebanon over three 385 contrasted snow seasons (2013)(2014)(2015)(2016) in three catchments that are critical in terms of water resources for Lebanon. This was made possible by the availability of new meteorological AWS observations at high elevation (1840-2834 m a.s.l.). The model was evaluated using a number of in situ and remote sensing observations. Recent developments with SnowModel including a gravity drainage formulation instead of the default snow density threshold to control liquid fluxes (Pflug et al., 2019) were beneficial in this study area where the snowpack is subject to multiple snowmelt episodes during the snow season. The results 390 indicated that this improved model was generally able to reproduce the observed snowpack properties with a reasonable performance despite the uncertainties in the model forcing. Although the default version of SnowModel can theoretically be applied in any climatic region since it relies mostly on physical equations to solve the energy and mass balance of the snowpack, it also includes parameterizations to replace processes that are too complex to be physically represented in the model. These parameterizations were primarily derived from studies in Arctic or subarctic regions (e.g. Bruland et al., 2004;395 Mernild et al., 2008) and should be carefully examined if the model is applied to Mediterranean mountain regions with mild, humid winters.
Over the simulation period 2013-2016, the maximum snow mass was reached between December and March and its magnitude changed significantly from year to year in the three study catchments. The results are subject to large uncertainties given that the model was only partially evaluated despite all the efforts made to collect in situ measurements. In particular 400 there is a mismatch between the spatial scale of field data and the current model extent and resolution which prevents a robust estimation of the model error. To further evaluate the model ability to resolve small-scale variability, higher resolution remote sensing products like Sentinel-2 SCA (Lebanon has been recently included in Theia Snow collection, Gascoin et al., 2019) or stereo-satellite snow depth (Marti et al., 2016) should be useful. Another important limitation of this work is the short record of the forcing data from the AWS that does not allow drawing robust conclusions on the inter-annual variability of the SWE 405 in Mount Lebanon and even less on the extremes. In particular, exceptional snowfall occurred in Lebanon during winter 2018-2019 and the snow cover was still significant in July near MZA and CED stations. Future work may focus on the application of downscaled climate reanalysis data to extend the study period (Mernild et al., 2017;. Longer simulations would also enable the use of the simulated snowmelt as input to a karst hydrogeological model in order to improve our knowledge on the snowpack dynamics influence on the spring discharges and regional water resources. 410 We found in the Mount Lebanon snowpack some typical features of the Mediterranean snow, i.e. the high temporal variability and the high densification rates (Fayad et al., 2017a). Results from this study indicate that SnowModel modified by Pflug et al. (2019) can provide realistic estimates of SWE, which are essential for understanding the hydrological contribution of the snowpack from Mount Lebanon. Yet, there are many aspects of the model that could be further improved:  The comparison with MODIS data suggests that the enhanced model performed better at mid-elevations but there 415 remains a positive bias in the snow cover duration at high elevation (Fig. 6). In this study, the model was forced with observations collected from three AWS but precipitation data were not available at the highest AWS (CED).
Since then a precipitation gauge was installed at CED and high elevation precipitation data are available. Therefore, we suggest that the simulations could be improved by using observed precipitation from all stations to have a better representation of the high elevation precipitation volumes. 420  Given that temperatures remain close to 0°C during winter precipitation events, over large areas of Mount Lebanon, a better parameterization of the precipitation phase partitioning is expected to enhance the simulations and help in better capturing the rain on snow events. In general, this is probably a key issue for Mediterranean mountain regions with mild winters.
 Blowing snow and its redistribution are important in this study domain. This is justified from the high variability 425 of snow depth at shorter distances as concluded from the field observations conducted by the authors during two winter seasons. The existence of a large number of sinkholes in the region act as a trap for blowing snow. Figure   23 11 is an example where snow depth varies from zero to more than 2.5 m over distances as short as a few hundred meters due to snow transport by wind. Previous work showed that the blowing snow and snow redistribution module in SnowModel improves the accuracy of the simulations (see Gascoin et al., 2013). Hence, we chose not 430 to focus on this aspect in this study and assume that it contributes to a realistic representation of the snowpack.
Capturing snow redistribution of snow in the future would probably require focusing on a smaller model domain with an accurate fine scale elevation DEM that is currently not available in this region. A proper representation of the snow redistribution could also help in better explaining the role of advective heat fluxes and improve the model performance in capturing SWE. 435  Snow sublimation can be a significant component of the snow mass balance as in other Mediterranean Mountains (e.g. Schulz and de Jong, 2004). In SnowModel there are two sources of sublimation: static surface sublimation and blowing snow sublimation. Figure 13 shows the distribution of the mean annual sublimation rate, which reaches 70 mm we in the upper area of the study domain, and accounts for a fraction of about 5 to 15% of the annual peak SWE. Our justifications for the relatively low sublimation rates over the study domain can be 440 associated with the (1) high relative humidity due to the proximity to the Mediterranean Sea, which reduces both sources of sublimation (2) the high snow densification rates that inhibit blowing snow sublimation (Liston and Elder, 2006). We are unable to qualify the simulated sublimation however, point (2) is compatible with our field observations that the snow surface at the mid-elevation is subject to multiple melt and refreeze events. Such snowpack surface limits the capacity for snow removal by wind and hence reduces snow sublimation. To conclude 445 the model suggests that the sublimation is probably not the main source of uncertainty in this area (precipitation remains the main uncertainty), but we suggest this aspect should be the focus of further investigation using in situ measurements by eddy covariance towers or lysimeters.

450
 The surface radiative forcing effect due to the deposition of mineral dust, such as those from the African Sahara and the Arabian Desert, is not included in SnowModel, but recent studies in regions with similar climate suggest that it can have a significant impact on snowmelt runoff (e.g. Painter et al., 2018). In fact, most of the energy used for snowmelt is shortwave radiation, with the occurrence of numerous dust storm events in this region, a proper 455 parameterization of surface albedo, through assimilation of remote sensing products, as for example, could help in better capturing the onset of snowmelt especially at mid elevations. Collecting information on the radiative and thermal fluxes and measuring the properties of the different snowpack layers could help in better representing these processes in models.

460
Finally, a specific feature of Mount Lebanon is that most of the snow accumulates on a high elevation plateau. This morphology implies that the seasonal snow cover could be highly impacted by a rise in the mean zero-degree isotherm altitude in winter, which determines the partition of precipitation into rain or snow. Given the importance of snowmelt for groundwater recharge and water supply in Lebanon it is urgent to evaluate the vulnerability of the snowfall regime to climate change in Mount Lebanon. 465