Impact of distributed meteorological forcing on snow dynamic and induced water fluxes over a mid-elevation alpine micro-scale catchment

. From the micro to mesoscale, water and energy budgets of mountainous catchments are largely driven by topographic features such as terrain orientation, slope, steepness, elevation together with associated meteorological forcings such as precipitation, solar radiation and wind. This impact the snow deposition, melting and transport, which further impact the overall water cycle. However, this microscale variability is not well represented in Earth System Models due to coarse resolutions, and impacts of such resolution assumptions on simulated water and energy budget lack quantification. This study aims 5 at exploring these effects on a 15.28 ha small mid-elevation (2000-2200 m) alpine catchment at Col du Lautaret (France). This grass-dominated catchment remains covered with snow for 5 to 6 months per year. The surface-subsurface coupled hyper-resolution (10 m) distributed hydrological model ParFLOW-CLM is used to simulate the impacts of meteorological variability at spatio-temporal micro-scale on the water cycle. These include 3D simulations with spatially distributed forcing of precipitation, solar radiation and wind compared to 3D simulations with non-distributed forcing simulation. Our precipitation 10 distribution method encapsulates the spatial snow distribution along with snow transport. The model simulates the snow cover dynamics and spatial variability through the CLM energy balance module and under the different combinations of distributed forcing. The resulting subsurface and surface water transfers are solved by the ParFLOW module. Distributed forcing induce a snowpack with a more spatially heterogeneous thickness, which becomes patchy during the melt season and shows a good agreement with the remote sensing images. This asynchronous melting results in a longer melting period and smoother hy- 15 drological response than the non-distributed forcing, which does not generate any patchiness. Amongst the tested distributed meteorological forcing that impacts the hydrology, precipitation distribution, including snow transportation, is the most important. Solar insolation distribution has an important impact in reducing evapotranspiration depending on the slope orientation. For the studied catchment mainly facing east, it adds small differential melting effect. Wind distribution in the energy budget calculation has a more complicated impact on our catchment as it participate to accelerate the melting when meteorological 20 conditions are favourable but does not generate patchiness at the end in our test case. yearly subsurface stock difference spin-up storage to evapotranspiration regime. We showed that the precipitation distribution had the largest impacts on hydrological behaviour because it favours the appearance of no-snow patches in the melting season. Shortwave distribution had a smaller effect at creating no snow patches in our simulation but enhanced the differential melting when combined with precipitation distribution. Our wind distribution simulation also induced melting spatial variability in the core of the melting period but reduces at the end of the melting period, as only taken into account for evapotranspiration processes. However, wind has a major role on snow re-distribution. This was well accounted for in our precipitation distribution algorithm using a laser scan 415 of the snow mantle during the core accumulation period which registered snow variability aligned with the prevailing wind direction. Furthermore, accounting for distributed solar incidence reduced incoming radiation in our catchment subsequently reduced the evapotranspiration. This led to higher runoff coefficients at catchment scale. In conclusion, this study shows the relative importance of small-scale processes in earth system models when run at hyper-resolution. The terrain substantially controls the hydrological behaviour in the mid-elevation alpine catchment for runoff generation but also for evapotranspiration,

(b) (a) summer winter (c) Figure 1. (a) overview of the study area at Col du Lataret, France, the small sub-alpine catchment is delineated in red with the outlet at the blue point. Green dot is the Flux-alp micro-meteorological station. (b) landscape views of the Lautare pass area in winter (January) and summer (July). (c) Catchment domain (84×42 grid cells at 10m resolution) with river branches (violet), elevation contours (green) and tree patches. Coloured pixels represent the distributed snow coefficients. The dotted area is the approximated footprint for the for daily spring wind directions.

Climate
The study area is a typical mid-latitude alpine environment. Figure 2 shows the meteorological observation for the simulated hydrological year starting the 11 November 2017 at the first snowy day to 10 November 2018. The catchment has a long winter season with 5 to 6 months of snowfall (Fig. 2a) and a complete snow cover. Flux'Alp meteorological station records a total of 100 1530 mm year-1 precipitation, out of which 970 mm was snow (dry) in the studied period. According to 2017-2018 weather data (simulation year for this study), the site-average temperature is 4°C. The site has below-zero winter conditions, with a -7.4°C mean February temperature recorded and a 14°C mean July temperature (Fig. 2b). Strong winds are common throughout the year, usually from the South-West direction along the mountain pass ( Fig. 2c and 4a). The temperature and specific humidity annual cycle are well phased (Fig. 2d). March is the most humid period of the year, while July is the driest. The studied 105 catchment has significant differences between summer and winter solar radiation (Fig. 2f). Additionally, mountains around the catchment reduce the sunshine period with projected shadow, especially during winter. These observations also served as forcing to the model.

Monitoring
Most of the monitoring on the site has started in 2012. It includes the temperature and humidity (CS215, Campbell Sc.), at-110 mospheric pressure (Setra CS100, Campbell Sc.), Wind speed and wind direction (Vector anemometer A100LK and W200P, Campbell Sc.), 4 components of net radiation (CNR4, Kipp and Zonen), snow height (SR50A, Campbell Sc.), and NDVI (Normalised Difference Vegetation Index) measured through Skye Instruments SKR1800. Since 2015, the site received eddy covariance sensors composed of the LI-COR LI-7200 close-path gas analyser and the HS50 Gill 3D sonic anemometer. Since 2017 the site has also add an OTT Pluvio rain gauge. Site setup, monitoring and data processing follows the ICOS (https://www.icos-115 ri.eu/) standards, and the site has entered this program since 2018. All measured variables are recorded at 15 min time step and then reduced to a 30 min time series for this study. The EddyPro Software is used to process the turbulent fluxes at the same 30 min time step following the ICOS recommendations (Hellström et al., 2016). The precipitation data from the rain gauge have been first processed to account for the lack of gauge shield (Klok et al., 2001).
The adopted algorithm follows as: where, P corr is the corrected precipitation (mm), P is measured precipitation (mm) at observation station, u is the wind speed at observation station in m s-1, a and b are correction factors, and are different for rain (a = 1.04, b = 0.04) and snow (a 155 = 1.18, b = 0.20).
To account for snow coverage spatial variability on the catchment domain, S c (x, y), the snow precipitation at (x, y) location was calculated using a snow coefficient map C s (x, y) from the ratio between the measured snow height at the gauge to the snow height measured through the laser scan on the same day (21/02/2018) at the end of the accumulation period (Fig. 1). The snow height calculated from the laser is the difference of apparent snow height (from laser scan) at the end of the accumulation 160 period and the digital elevation model (DEM) for the surface without snow. The snow DEM and surface DEM were prepared at the resolution of 2 m and were upscaled to 10 m resolution for modeling purpose. The S c (x, y) calculation hypothesizes that distributed snow cover on that date aggregates all spatial heterogeneity of the snow deposition including the snow transportation (redistribution). It also neglects the snow compaction between date of deposition and the laser scan. Then the corrected Snow precipitation is calculated according to: where, S m (x, y) is the measured snow precipitation at the observation station. Unfortunately, it can be noted that the laser scan didn't cover the upper part of the catchment. Zones not covered by the scan were each given a fixed value, according to our field observation. Moreover, due to micro-scale catchment and a low altitudinal range, elevation differences between upper and lower locations ( 211 m) have not been considered and the rain has not been distributed according to an altitudinal gradient.

Shortwave radiation
The shortwave radiation, SW c (x, y) has been distributed from the observed shortwave measurement, SW m (x 0 , y 0 ) at the meteorological station considering the sun position and the terrain effect (Liston and Elder, 2006). shortwave was taken from the CLM technical setup (Oleson et al., 2004).

Wind speed
Wind speed were spatialized to better account for the estimation of turbulent fluxes (Liston and Elder, 2006). The wind speed was distributed as, map. Snowmelt dynamics were compared to Sentinel 2A and Sentinel 2B products from Sentinel satellite obtained through 210 USGS Earth Explorer (https://earthexplorer.usgs.gov/).
The lateral and bottom boundary conditions were set to no flow and the surface boundary condition was an atmospheric pressure condition that allow surface runoff (Kollet and Maxwell, 2006). Hence, the inflow and outflow are restricted to exchange only through the surface. Subsurface has been made heterogeneous with three layers of soil, regolith and flysch divided into 11 numerical cells (Fig. 3). The bottom of the domain was set deep enough to accommodate various subsurface movement 215 (118 m deep from the surface). The soil physical parameters used in this study includes porosity, permeability, soil horizons and Van Genuchten parameters. Soil horizons distribution were determined from an electromagnetic survey measuring apparent electrical conductivity (related to water and clay content) and ground penetrating radar (GPR) measuring soil thickness along key-profiles, while the soil properties were determined by field permeability experiments (Reynolds and Elrick, 1991;Vandervaere et al., 2000a, b) and laboratory mercury porosity experiments. Elaboration about the detailed hydro-geological 220 characterization is beyond the scope of this study and will be detailed in a companion paper. This study is more focused on surface dynamics due to meteorological distribution. the effect of spatially distributed forcing all together or individually. Meteorological distribution algorithms described above serves the purpose of catching the slope, curvature and aspect effect in spatial distribution. Even at a micro-scale, one can observe the spatial meteorological variability along the grid after applying equation 2-8 ( Fig. 4). In wind distribution, there is not so much variation. The mean wind speed before and after the spatial distribution is 3.6 m/s and 3.5 m/s, respectively. melting is simulated because of warmer positive temperature ( Fig. 2b) but also rain on snow events. This produces the highest river discharge peaks on the hydrograph following the daily cycle (between min and max discharge values at daily timestep) and the rain timing (highest discharge peaks). This period also increases the subsurface stock (Fig. 5b) which produces a base flow that combines with the snowmelt contribution to streamflow. One of the most important and noticeable point using the non-distributed forcing is the sudden disappearance of the snow 260 at the end of the snowmelt season, which is usually not observed on actual field observation. This means that all the pixels behaved the same way and there is no noticeable impact of the subsurface hydrology on the catchment spatial variability when considering a uniform forcing. From June to the beginning of the next snow period, summer rains produce almost instantaneous river responses and subsurface stock sustains stream discharge for several months. During this period, one can note a radical change of net radiation because of the change of the albedo from snow cover ( 0.8) to herbaceous ( 0.2) (see Fig. 9(a)). Earlier 265 in spring net radiation contributes to snowmelt (the two variables are correlated) because of higher sun elevation more clear sky conditions and higher daily temperature.

Results and Discussions
During winter and spring monthly cumulated ET is very small (Fig 5b) because of low available energy and complete snow coverage. After the complete snowmelt the model simulates much higher monthly cumulated ET according to the prescribed LAI cycle. ET at this period is higher than the monthly cumulated rain (June, July, September), which means that ET partic-

270
ipates to the extraction of shallow water storage during the summer. This can be seen by the difference of subsurface storage decline between the summer (higher water storage diminution) and the winter (lower water storage diminution). In October one can notice a small subsurface stock increment when ET reduces because of vegetation decay. At the end of the hydrological year, the subsurface water storage has a low negative value of 0.15 mm which is the hydrological unbalance much smaller than the annual cycle.
275 Figure 5c-d present the same time series and monthly water budgets for the 1D-AM simulation. The major difference is that precipitation, solar radiation and wind velocity have been prescribed using the spatial average of the distributed field. The major difference comes from solar radiation reduction from 190.8 W m-2 to 152.1 W m-2 on average, which reduces melting and ET. It leads to a 9 days longer snow period, an increased runoff and an increased infiltration of 34.61 mm (Table 1). This unbalance vanishes after several years leading to higher water table and then higher runoff with a 0.8 runoff coefficient. This 280 means that for similar geomorphology, any reduction in input solar radiation because of orientation or else will lead to higher water table and then higher runoff coefficient.

Distributed forcing simulations
In the following section, we will discuss the differences in surface fluxes for all combined and individually distributed simulations, along with their role in streamflow generation.

285
Looking at the surface fluxes on Fig. 6a, one can see that the simulation with all distributed forcing (2D-AD) have longer snowmelt period with a long steady decline during the streamflow recession till mid-July. Compared to the non-distributed forcing simulation 1D-PM, the distributed forcing is causing a smoother snowmelt dynamic which last till July and correspondingly impacts the net radiation recharge and streamflow discharge dynamic with significant differences for streamflow values during the core of the melting period up to 50 % more for 1D-AM compared to 2D-AD. However, budget terms are 290 pretty much the same between 1D-AM and 2D-AD. From this point of view distribution at that scale seems not to add much information compared to a 1D forcing. Looking at the surface fluxes on Fig. 6a, one can see that the simulation with all distributed forcing (2D-AD) have longer snowmelt period with a long steady decline during the streamflow recession till mid-July. Compared to the non-distributed forcing simulation 1D-PM, the distributed forcing is causing a smoother snowmelt dynamic which last till July and corre-295 spondingly impacts the net radiation recharge and streamflow discharge dynamic with significant differences for streamflow values during the core of the melting period up to 50 Looking at individually distributed simulation results, it seems that this smoothing effect is caused by both the precipitation and short wave radiation distribution (Fig. 6b, c). At the contrary only wind distributed (2D-WD) results seems very similar to the non-distributed results. One can even refer to (Fig 6d) for easier inter-comparison with (Fig. 5a). The melting period tail 300 length seems to be controlled by the snow pack depth variability (Fig. 9a, b) and higher accumulated snow on some pixels. This is combined with solar radiation effects which also produces melting spatial variability on the catchment even if the snow precipitation is uniformly distributed (Fig 6c). Only wind distribution (2D-WD) simulation, showed the highest melting regimes from mid-March to mid-May when temperatures and incoming radiations are favourable for melting with daily melting peaks upper than 4mm h-1 (Fig. 6d). In detail, wind distribution increase melting rates which led to higher runoff (Fig. 8d) 305 when compared to 1D-PM.
Streamflow differences between simulations basically follow the melting differences. The 20 day April rain-on-snow period is particularly marked on streamflow on Fig. 6a, b. It has to be reminded that there are differences in term of incoming solar radiation between simulations. Unsurprisingly 2D-WD and 2D-PD simulation show larger streamflow values compared to 2D-the summer period, it is difficult to see any differences on streamflow.  The pixel run (saffron curve and saffron dots) clearly overestimated observed evapotranspiration as one can see on both the time series and the scatter plots. This is also the case for the Non-distributed simulation (green curve and green dots) and for 2D-PD and 2D-WD which have comparable evapotranspiration amplitude (Fig 7c, e). However, 2D-AD and 2D-SD have 320 reduced simulated ET which better matches the observations (Fig 7d-b). The cause is that the average shortwave after the distribution is less than the shortwave without distribution (section 4) and the catchment is facing east which reduces direct incoming solar radiation from noon to sunset.
The evapotranspiration of both Pix-PM and 1D-PM over-estimate ET compared to observations. First the pixel run (Pix-PM) is supposed to simulate a catchment border (Flux'Alp location) with dryer soil/ground condition (top of a catena) and the ET observations are supposed to average both the wet zone close to the river and dryer zones. To the first order, it is not the case in our pixel and 1D-PM simulations. Shortwave distribution (Fig. 7d) seems to have the most important impact in our measurement area. The corresponding reduced ET in 2D-SD (and 2D-AD), averaged on the footprint area, also corresponds much better to the Eddy-Covariance observations.

Water budget 330
Annual water budgets in (Table 1) shows that the main impact of forcing distribution is caused by the shortwave distribution and subsequent ET calculation. This increases the runoff coefficient from 0.73 to 0.77. One can note the water storage change over the year. As we started from the same initial conditions for all simulations, it reaches more than 30 mm when SW is reduced (1D-AM, 2D-AD and 2D-SD) and 15.5 mm for 2D-WD simulation. It remains smaller than the ET difference and it minimizes runoff coefficient hence not change our conclusions.
335 Figure 8 shows monthly water budgets for 2D-AD and individually distributed simulations. Snow precipitation from November to March doesn't infiltrates much to fill up the subsurface stock (dotted line). January rain on snow events slightly reduces the subsurface drainage. Very similar runoff values are observed whatever the forcing. In contrast, from mid-March to June the subsurface stock is replenished by melting (Fig. 5a, Fig. 6) and runoff increased. The 2D-WD forcing produced the largest values of recharge (430 mm) and the 2D-SD the largest values of streamflow. From May to October, streamflow at the outlet 340 and ET drawdown subsurface stock. Higher SW radiations (2D-PD and 2D-WD) led to longer ET periods. One can finally note that reduced ET because of vegetation senescence in November leads to increased subsurface stock. for compaction. The snowmelt dynamic is particularly well simulated, especially along the dry period at the end of April. In early May one can note some discrepancies again probably because of our limited ability to separate rain and snow in the precipitation forcing, close to the phase change temperature range. This can be seen on the pixel simulated Albedo which 350 return to its snow value at the end of the melting season (0.8), which is not the case on the observations. Concerning simulated albedo, it globally follows the observations, however it is clear that the snow age parameterisation in the model is not adequate enough to simulate the albedo where observation does not show albedo decrease during dry periods.

Snow dynamics
Blue line on figure 9a shows the average snow depth for the all-distributed simulation (2D-AD) and the sky blue shading in the same figure shows spatial variability between min and max snow depth over all pixels. For this simulation patchiness starts 355 early in May and some pixels snow coverage last up more than one month later than the 1D-PM simulation. Depth variability is mainly increasing during accumulation events and is slightly reduced during melting. This effect can also be seen on the 2D-PD simulation but not on the other ones (Fig. 9b). As 2D-AD and 2D-PD simulations were prescribed the same input radiation and temperature, this means this effect (the deeper the snow, the faster the melting) is intrinsic to the snow scheme.
At the contrary, 2D-SD simulation shows a slight increase of depth variability during the melting period.

360
It can be observed in (Fig. 9b) that none of the individually distributed simulations have longer snow coverage compared to the all distributed run (Fig. 9a). It indicates that snow variability during accumulation events is not enough to capture the actual behaviour of snow dynamics. The combination of both forcing results in longer snow periods with more patchiness accounting for precipitation spatial variability during accumulation events and shortwave spatial distribution dominance for differential melting. The 2D-WD simulation does not show much impact on the snow depth variability which is very similar to the 1D-PM 365 simulation at the end of the snow period (Fig. 9b). However along spring (mid-March to end of April) it produced the same snow depth spatial variability as 2D-SD and higher snowmelt regimes ( Fig. 9b and 6d). Wind distribution is then impacting snow patchiness through wind transport (accounted for in the snow distribution algorithm) and is increasing melting through the energy budget.
Spatial distribution of snow coverage during the melting period is shown on figure 10 for all simulations. For comparison 370 and validation, we used four cloud-free 'Sentinel-2' images (Table 2). In the model, the snow pixels were selected for snow depth threshold over 1 cm. The Sentinel snow pixels were selected with normalised snow difference index (NDSI) > 0.4 (Riggs et al., 1994).
The 21st of November, first snow events were followed by a partial melting over the catchment (1st sentinel image). Our 2D-AD and 2D-PD simulations show a partial success in representing this feature, but the simulated melting was overall too   small. It is interesting to see that, apart from the upper part of the catchment where snow distribution is not well controlled, the first pixels that get out of snow are the eastern edge of the catchment, a central area aligned with the river, and the outlet area.
The 2D-AD simulation has less green pixels because reduced incoming radiation caused by reduced solar angle, reduced the melting. On 6th of December, the catchment is completely covered by snow for all simulation. This date still also corresponds to early season snow events. The 2D-AD and 2D-PD simulations are able to represent the snow dynamic even for very low 380 snow depth at the beginning of the season. This means in particular that 1) our spinup process has initiated well the ground temperature profile and distribution and 2) our distribution algorithm is well adapted, especially for snow deposition. The 25th of May is located when snow partially cover the catchment during the melting period. Again 2D-PD simulation represent very well the snow pixels to non-snow pixels ratio and the snow distribution. One can see on both Sentinel image and 2D-PD simulation some SW-NE alignment slightly present on the snow distribution coefficient map (green/blue pixels on Fig. 1c), 385 the timing of disappearance is remarkably well simulated for this simulation. The 2D-AD is not as good because of more snow coverage than the sentinel image. Because of a reduced melting when taking into account the solar incidence, this run is