Mass balance and hydrological modeling of the Hardangerjøkulen ice cap in south-central Norway

15 A detailed, physically based, one dimensional column snowpack model (Crocus) has been incorporated into the hydrological model WRF-Hydro to allow for direct surface mass balance simulation of glaciers and subsequent modeling of meltwater discharge from glaciers. The new system (WRF-Hydro/Glacier) system is only activated over a priori designated glacier areas. This glacier area is initialized with observed glacier thickness and assumed to be pure ice (with corresponding ice density). This allows for melt of the glacier to continue after all accumulated snow has melted. Furthermore, the simulation of surface 20 albedo over the glacier is more realistic as surface albedo is represented by snow where there is accumulated snow, and glacier ice when all accumulated snow is melted. To evaluate the WRF-Hydro/Glacier system over a glacier in Southern Norway, WRF atmospheric model simulations were downscaled to 1 km grid spacing. This provided meteorological forcing data to the WRF-Hydro/Glacier system at 100 m grid spacing for surface and streamflow simulation. Evaluation of the WRF downscaling showed a well comparison with in situ meteorological observations for most of the simulation period. The WRF-Hydro/Glacier 25 system reproduced the glacier surface winter/summer and net mass balance, snow depth, surface albedo and glacier runoff well compared to observations. The improved estimation of albedo has an appreciable impact on the discharge from the glacier during frequent precipitation periods. We have shown that the integrated snow pack system allows for improved glacier surface mass balance studies as well as hydrological studies.

of user-set layers, the thickness of the snow layers and the snow grain characteristics), these snow layers will merge into single 125 snow layers.
The Crocus module is added to the Noah MP land surface model inn WRF-Hydro to act as a glacier mass balance model. Over designated glacier grid points, the Crocus snow model represents both snow and ice, while outside of the designated glacier grid points, the regular three-layer snow model in Noah-MP is used. Since the current Crocus implementation in WRF-Hydro only acts over designated glacier grid points, we follow Gerbaux et al. (2005), and assume that the temperature at the bottom 130 of the glacier and the ground below are both at 0 ºC. Note that we have no yet incorporated fluxes between the glacier and the ground below, thus there is a constant-temperature boundary conditions. Both Crocus and Noah-MP (for the non-glacier grid points) output runoff from snowmelt (and precipitation). This runoff is provided to the terrain routing models in the hydrological model system WRF-Hydro. WRF-Hydro is a model-coupling framework designed to link multiscale process models of the atmosphere and terrestrial hydrology (Gochis et al., 2015;Yucel 135 et al., 2015). In coupled mode it includes the full functionality of the atmospheric Weather Research and Forecasting (WRF) modelling system. WRF-Hydro enables simulation of land surface hydrology and energy states and fluxes at high spatial resolutions (typically 1 km or less) using a variety of physics-based and conceptual approaches (Yucel et al., 2015;Senatore et al., 2015). It contains horizontal routing processes and water management modules and is linked to the Noah-MP land surface module (among others). The added capability of running Crocus as a glacier mass balance module in WRF-Hydro is 140 called the "WRF-Hydro/Glacier" system from here onward.
Note that our implementation of Crocus as glacier mass balance model does not address glacier movement (i.e. plastic flow) nor lateral wind (re)distribution of snow. Being a relatively flat dome glacier, the Hardangerjøkulen glacier that we focus on here does not move much in a year, and therefore the lack of dynamical movement of the glacier is not expected to have a major impact on the results in this paper, as we only consider four simulation years. On the other hand, the lateral snowdrift 145 and wind-driven redistribution of snow on Hardangerjøkulen can be significant, and our results are likely impacted by the lack of this physical process in the model. It is worth mentioning that including lateral movement of snow due to snowdrift in the model system is not a trivial task, and is therefore not currently included. However, there are two options to include impacts on the snow due to wind. One of the options impacts the snow density during blowing snow events (Brun et al., 1997). This option is important in polar environment (Brun et al., 1997), and we found it necessary in our simulations as well. The other 150 option is the sublimation due to snow drift, which was implemented by Vionnet et al. (2012) and which is in the Crocus version that is used in this study. This option was also crucial to include in order to accurately simulate the glacier mass balance over Hardangerjøkulen. It is especially important for reproduction of the observed heterogeneous snow distribution.

Initialization of glacier module in WRF-Hydro/Glacier
To run the WRF-Hydro/Glacier system, the glacier to be evaluated must be initialized with its thickness and extent. Here we 155 focus on the Hardangerjøkulen; its extent and height were obtained from the NVE (Melvold et al., 2011). Figure 1 shows the glacier thickness and extent of Hardangerjøkulen at 100 m grid spacing (for which the entire WRF-Hydro/Glacier model is run). At initialization, it is assumed that the glacier consists of only ice, and the density is that of pure ice (900 kg m -3 ). In the simulations presented here, the user-defined maximum layers are set to 40 layers, and the glacier is initialized with all the layers having the same assumed density and snow grain properties. As new snow accumulates during the simulations, the 160 layers representing the glacier will start to merge since all 40 layers are occupied with the initialized ice. Here, the thickest parts of the glacier merged to an average of 8 layers after 5 months of simulation and remained fairly constant for the remainder of the simulation. Revuelto et al. (2018) also used Crocus for surface mass balance studies. They initialized their glacier with the same glacier thickness (40 m) over the entire glacier and in the six lowest layers with the thickness progressively increasing with depth. In contrast to this study, they reinitialized the glacier to 40 m every new season (August 1) so that the glacier would 165 not decrease in extent, while here the glacier is only initialized once at the beginning of the simulation with the observed glacier thickness.
As implemented, if the glacier completely melts over a user defined glacier grid point, the original Noah-MP module is used from this point on. Therefore, as currently implemented, the glacier cannot grow horizontally in extent, it can only decrease in extent, as no dynamic response of the ice mass is included in the model. Over short model time periods, the lack of increase 170 in glacier extent might impact a few grid points at the edges of the glacier. However, given the expected increase in temperature in the future we do not expect that limiting glacier horizontal growth will have a major impact over most studied glaciers as most are likely to decrease in mass and extent.

Experimental Description
For meteorological input data to the WRF-Hydro/Glacier system, dynamically downscaled data from the WRF model version 175 3.9 (Skamarock et al., 2008) was created over Southern Norway. WRF was run with an outer domain (Domain1) with a 3 km grid spacing and an inner domain (Domain 2) with a 1 km grid spacing (see Figure 2, top) with 51 stretched vertical levels (lowest prognostic level is 25 m). The ECMWF Re-Analysis Interim (ERA-I) dataset was used for input and boundary conditions and the model was run from August 1, 2014 to January 1, 2019. The microphysics scheme used was the Thompson-Eidhammer aerosol aware scheme (Thompson and Eidhammer, 2014), the Yonsei University (YSU, Hong et al., 2006) scheme 7 for the boundary layer (Hong et al., 2006), the rapid radiative transfer model (RRTMG) for longwave and shortwave radiation calculations (Iacono et al., 2008), and the Noah-MP land surface model (Niu et al., 2011). See Table 1 for the configuration.
We did not use any lake models, thus, in WRF, the skin temperature of lakes is typically set to the same temperature as the nearest grid point that is defined as sea surface. With this setting, the lakes will not reach freezing temperatures since the oceans surrounding southern Norway typically do not freeze in the winter. To rectify this problem, we assign a 10-day moving 185 average skin temperature from the associated ERA-Interim grid point onto the lake grid points to allow a representation of freezing lakes. The fjords still use the sea surface temperature from ERA-Interim. We acknowledge that a large step from ~75 km (ERA-I) to 3km in WRF is of concern. However, we follow findings by Liu et al., (2016) where they state: "Tests showed that one-way nesting WRF, at 4-km grid spacing, with the ~75 km reanalysis was an adequate configuration without the need for a coarse grid that intermediates the ERA-Interim data and the WRF domain." What is important is that the area of interest 190 must be sufficient large enough for mesoscale spin up. Our domain of interest (Domain 2) is slightly closer to the boundary than what is in Liu et al., (2016). However, as shown below, the model results are quite reasonable, thus we believe that the jump from ~75 km to 3 km is adequate.
The 1 km WRF (inner-nest) simulation results are used as input to run WRF-Hydro/Glacier. The WRF-Hydro/Glacier domain has 100 m grid spacing (Figure 2, bottom) and covers a smaller area compared to that of the 1 km domain. The precipitation,195 2 m temperature, 10 m wind speed, 2 m water mixing ratio, surface pressure, and long and shortwave radiation outputs from the WRF 1 km simulations were bi-linearly re-gridded to the 100 m grid spacing in the high-resolution domain. We note that we did not account for variability in terrain in the re-gridding process. Thus the atmospheric forcing is still "smooth" as regards to a 100 m grid. However, the region of interest (Hardangerjøkulen and surrounding terrain) is an open, mostly flat area and we therefore believe that for this specific case, disregarding the variation in terrain does not have much impact on mass balance 200 calculations. For partition of rain/snow from the input precipitation, WRF-Hydro/Glacier use the rain/snow partition from Jordan (1991). The streamflow routing is run at the same resolution as Crocus (i.e. at 100 m). Finally, no calibrations were applied to the routing model.

Glacier observations
Hardangerjøkulen is a well observed glacier, with several decades of mass balance observations, and several field campaigns. 205 In the following, data from field observations, the ongoing mass balance observations and remote sensing are used to evaluate the WRF-Hydro/Glacier system.

Glacier mass balance
Glacier mass balance is the amount of mass a glacier gains or loses over a year (sum of accumulation and ablation). NVE has gathered winter, summer and annual (net) mass balance (winter + summer mass balance (where summer mass balance is 8 negative)) observations over Rembesdalskåka since 1963(e.g. Andreassen et al., 2020 where Rembesdalskåka is the west/south-west glacier outlet of Hardangerjøkulen (see Figure 1). The observations are gathered at several locations on the glacier, from the lower to the upper parts (1066-1854 m a.s.l.) by using stakes (aluminum poles inserted into the glacier) and in spring by probing (using thin metal rods to measure snow depth) to the previous year's summer surface. Winter mass balance is found from the observed snow depth (by stakes and rod soundings at approximately 60 locations (see

Radar-derived snow thickness
Variations in snow accumulation were measured over Hardangerjøkulen in April 2017 and 2018, using ground-penetrating 225 radar (GPR). In 2017, surveys were conducted with a MALA Geosciences GPR system; in 2018, this system was not available hence a Sensors & Software pulseEKKO PRO model was used instead. However, data from these two systems are directly compatible since both were acquired with antennas of 200 MHz center-frequency. The GPR systems were towed behind a snowmobile at ~15-20 km/h. The interval between successive GPR recordings is ~0.2 s, giving a distance sampling interval of ~ 1 m (regularized to exactly 1 m in processing). A GPS system was also mounted on the snowmobile, recording positions 230 every 1-2 s, to locate the GPR recordings. The positional accuracy of the GPS is ~ 5 m. A total of 116 km and 27.4 km of measurements were acquired in 2017 and 2018, respectively, and both acquisitions featured numerous crossing-points such that the internal consistency of accumulation estimates could be ensured.
GPR systems record the travel-time of a radar pulse through the ground, therefore estimates of winter accumulation requires some measure of the GPR propagation velocity to covert time to depth. This was obtained using so-called common midpoint 235 (CMP) data (e.g., Booth et al., 2011Booth et al., , 2013, in which GPR responses suggested an average velocity of 0.218 ± 0.001 m/ns for the upper ~ 2.8 m of the snowpack. With no other velocity control available, this value is applied to convert all GPR traveltime estimates to a snow depth.
The base of winter snow accumulation was taken to be the first prominent reflective horizon within the GPR record. This is straightforward in the Hardangerjøkulen ablation zone, typically at elevations ~< 1600 m, where winter snow directly overlies 240 the glacier surface, typically at a depth of 2-3 m. Here, the only significant reflection is from the glacier surface itself.
Consequently, depth errors are expected to be less than ± 0.1 m. However, areas of firn accumulation have a more complex pattern of reflectivity and it is not always possible to guarantee accurate snow thicknesses, and errors may here be up to ± 1 m. However, given the crossing points in the GPR records, depth estimates are at least internally-consistent and errors are expected to vary systematically across the entire record.

Streamflow
Discharge measurements were obtained at two rivers. One river (here named "Middalselvi") is fed by meltwater from the Midtdalsbreen (a glacier arm of Hardangerjøkulen) where the catchment is about 12 km 2 and 60% glacierized (see Figure 2).
The other river is Finseelvi, where the catchment is about 16 km 2 and 14.7% glacierized, and not much impacted by glacier melt. Two Hobo Water Level loggers were installed in each catchment in Fall of 2016 and we have data until November 2018. 255

WRF standalone verification
The and Geilo) are located at high-altitude exposed locations. At these stations there is a large under-catch of observed snow in the wintertime when the snow can blow past the precipitation gauges instead of falling into the gauges (Rasmussen et al., 2012).
The stronger the wind, the larger the under-catch. The data obtained from the Norwegian Meteorology Institute for these stations have not been corrected for any under-catch. Figure 5 show the effect of under-catch of snow at Finse, which is the 275 station that is located closest to Hardangerjøkulen, about 4 km north-northeast of the edge of Midtdalsbreen and about 11 km from Rembesdalskåka. In Figure 5, the accumulated precipitation and temperature for the 2016 summer and winter season is shown at both Finse, Fet (a station about 25 km south west of Finse, but about 14 km away from Rembesdalskåka) and Evanger (a non-exposed inland station with little snow). Similar results are seen for 2015, 2017 and 2018 (not shown). As can be seen, WRF compares well with observations at Evanger and Fet almost the entire period (Figure 5a-5d). WRF precipitation also 280 compares well with observations at Finse in the summer season (Figure 5f), but has much more precipitation than the observations in the winter season ( Figure 5e) where wind speeds are often more than 5 m s -1 (Figure 5i) and temperatures are below freezing. Furthermore, WRF does simulate the storm sequences well as seen in Figure 5g and 5h. Thus we attribute the low observed precipitation compared to WRF during the winter season at Finse to the under-catch problem as precipitation modeled at Fet compares great with observations. Note that during September to 17 to September 26, the Finse station did not 285 provide any data (Figures 5f, h and j). However, during this time period, WRF did not predict any precipitation, and Fet did not observe any precipitation. Thus, the cumulative precipitation shown in Figure 5f is still valid.
The World Meteorology Organization (WMO) Solid Precipitation Intercomparison Experiment (SPICE) was set up to evaluate the under-catch of snow and develop transfer functions to correct for the under-catch of solid precipitation. (Smith et al., 2020). One location for these studies is Haukeliseter, which is about 20 km from Røldal (see Figure 2) and is within the 1 km 290 domain. In these studies, several different precipitation gauges and wind shield combinations were used. The Double Fence Automated Reference (DFAR) was deployed as the reference and is used as the "truth" precipitation. We compared the WRF model results with the DFAR data (Smith et al. 2019), and WRF is predicting more precipitation compared to these observations, with a bias typically at ~30% (not show). About 10% of this bias could potentially be attributed to underestimation with the DFAR (Rasmussen et al. 2012). The bias in WRF is opposite and higher compared to what is found 295 at locations with little impact of snow. In regards to transfer functions (correcting for under catch in observations), Smith et al (2020) stated this is their study: "Although the application of transfer functions is necessary to mitigate wind bias in solid precipitation measurements, especially at windy sites and for unshielded gauges, the inconsistency in the performance metrics among sites suggests that the functions be applied with caution." We are therefore not adjusting the observed observations on Finse for our evaluation, and rather stress the well comparison between model and observations at Fet and summer season 300 precipitation on at Finse. We conclude, as shown above, most of the non-exposed inland stations are relatively well simulated by the WRF-model and the seasonal cycle of precipitation is captured. We note, however, one time period where WRF is underpredicting the precipitation relative to locations near Finse (Hardangerjøkulen). Several stations near Finse do not catch a precipitation period in the middle of January 2017, which has an impact on Finse as well (this underprediction in WRF is difficult to directly 305 evaluate at Finse, but the two stations near Finse (Fet and Eidfjord) clearly have an underprediction at this time period (not shown). We note that the observed storm sequences were captured in the simulation, wind direction was well simulated as well as the wind speed during this precipitation event (not shown), just not the precipitation amount. The effect of this precipitation time period will be discussed in relation to the mass balance and streamflow results is Section 5. Figure 6 shows a scatterplot of observed and simulated 2 m temperature at the Finse AWS. As can be seen, the modeled 310 temperature compares well with the observed, but with a small negative bias in the winter. At the very low observed temperatures ( < -15ºC), WRF often has a positive bias. This is likely due to WRF not capturing the strong inversions that often occurs in the winter months (Mölders and Kramm, 2010;Hines et al., 2011). Figure 5e also shows this positive bias at the very low temperatures. However, in general, WRF compares well with observations with a correlation coefficient near 0.9. Figure 7 shows the simulated and observed 10 m wind speed and direction for the entire simulation period at Finse. Although 315 the wind direction is not an input variable in the land surface model in WRF-Hydro/Glacier, the wind direction can dictate precipitation amount and type, thus wind direction is important for obtaining correct simulation of precipitation and subsequent glacier mass balance as shown in Bhatt et al (2020) 1 . As can be seen, the simulated wind direction compares well with the observed. Mean bias is -0.13 m s -1 and coefficient of correlation is 0.8. Overall, the major simulated wind directions and speed are simulated quite well (see also Figure 5 for wind speed). 320 winter (and summer) balance is plotted as averages of all grid points over Rembesdalskåka within intervals of 30 m. As can be seen, the modeled mass balance (winter, summer and net, Figure 8a, d, g, and j) is generally comparable with the observations. The observed winter mass balance shows a small decrease at the top of the glacier, with a slight increase about 200 meters below the top (see left panels). This is most likely due to redistribution of snow that is common at Finse and its surroundings, as strong winds occur often there. Since lateral redistribution of snow is not included in the Crocus and Noah-330 MP models, the modeled mass balance increases more linearly towards the top compared to the observations. The winter mass balance is about 29 % (20 %) underestimated in 2017 by Crocus (Noah-MP), which is likely due to the underprediction of winter precipitation at Finse and nearby stations in this year (discussed in section 3.1). For the three other years (2015,2016 and 2018), the modeled winter balance is within 9.3, 13.8 and 5.1 % (8.7, 14.7 and 2 %) of the observed winter mass balance for Crocus (Noah-MP) respectively (see Figure 9). The reason for the slightly larger bias for Crocus is not known at this point. 335

Glacier mass balance and snow height
However, note that Crocus matches better with the observations at lower levels below about 1500-1600 m in 2017 and 2018 (see Figure 8 left and middle panels). This does not show up in the total mass balance comparisons since most of the mass is above 1600 m due to higher snow thickness and larger area (see Rembesdalskåka, Figure 1). Another point to add here is the improved simulation of winter mass balance of Rembesdalskåka compared to Engelhardt et al (2012). They used gridded interpolated precipitaion data from seNorge (http://senorge.no) and obtained a mean negative bias of 28% from 19 modeled 340 winter mass balance years for Rembedalskåka, while we use high-resolution regional scale modeling to obtain horizontally distributed precipitation.
For the summer and net mass balance we only consider Crocus since Noah-MP cannot melt more than the accumulated snow amount from the model start (see Figure 8, left panels). As the winter balance, the modeled summer mass balance is also in general agreement with the observations with a bias of 6.67%, -1. 5%, -7.4 % and 5.9 % in 2015, 2016, 2017 and 2018 345 respectively ( Figure 9). However, the modeled 2018 summer balance curve in Figure 8j  The modeled distribution of snow thickness and snow water equivalent (SWE) over Rembesdalskåka is comparable to 355 measurements by NVE ( Figure 10). The NVE data are the same as used in Figure 8, while not averaged in elevation transects as in Figure 8. There is more heterogeneity in the observations, as would be expected since the models does not account for lateral redistribution of snow. Despite this, Crocus does account for sublimation of snow drift, which results in more heterogeneity in the spatial snow distribution over the glacier compared to Noah-MP. The relatively large underestimation of modeled SWE (i.e. the winter balance) in 2017 is also seen in the snow thickness. This is also evident in Figure 8 Noah-MP are also lower, but higher than the Crocus snow densities. In the first modeling year of 2015 Figure 8c), the snow 365 density with Noah-MP is uniform over all elevations. However, as the simulations continued over several seasons (Figure (8f,   8i and 8l), the snow density increases with height for the following years. The reason for this increase is that Noah-MP only has three model layers. When new snow in the winter accumulates, some of this snow is merged with higher density multiyear snow from previous years in the accumulation zone of the glacier. It is therefore difficult to estimate the actual modeled snow density of the new seasonal snow layer in the accumulation zone with Noah-MP in a multiyear simulation. The underestimation 370 of density with Crocus is in agreement with Queno et al., (2016) who found that the bulk density in their study using Crocus also was underestimated and that Crocus tends to underestimate the snow compaction. On the other hand, since the density used to approximate the observed mass balance is taken from only one location, there are likely some uncertainties in the representability of this density over the entire glacier.
Ground penetrating radar observations were gathered in April 2017 and 2018. Note that this time period is not the same as 375 when winter mass balance is determined with stake observations which occurred at the end or after the accumulation season ( Figure 10). Figure 11 shows both modeled and observed spatially distributed accumulated snow for the respective winter season, and the snow depth as a function of height. Figure 12 shows a 2-D histogram of the same data. Although these data show some scatter, and isolated observations significantly diverging from the model output (blue), the occurrence of such points is between 10-100 times lower than those which plot close to the one-to-one comparison (red). As such, the majority 380 of depths observed within the GPR dataset match the model outputs to within 1 m. Crocus has more variation in snow depth over the glacier compared to Noah-MP (see Figure 12), but the modeled glacier in both snow models has less heterogeneity compared to the observations (Figure 11). For 2017, we can see that there are areas over Rembesdalskåka where the model clearly underestimates the snow thickness, as can also be seen with the point observations from NVE in Figure 10. However, at other locations on Rembesdalskåka the comparison is favorable. Also, in the northern and western part of Hardangerjøkulen 385 (the entire glacier complex) the model matches the observed snow thickness quite well. In the middle part of the glacier, Crocus estimates deeper snow depth compared to Noah-MP and is closer in agreement with the observations both in 2017 and 2018.
The reason for the lower snow depth in Noah-MP is likely due to the high snow density at higher elevations (see Figure 8) In section 2.1 we mentioned the importance of adding sublimation of blowing snow in our simulations. Figure 13 shows the snow thickness and the respective scatterplot for when the sublimation of blowing snow is not included (which is the default 390 in the SURFEX V8.0 setup (as downloaded)). As can be seen, the simulated snow thickness is slightly higher than the observations with the GPR, and this is especially true at the eastern part of Hardangerjøkulen. During the 2017 winter season, this region had on average the strongest winds, causing more sublimation from snowdrift than at other locations on the glacier.
Without including sublimation of blowing snow, the simulation overestimates snow thickness. However, for Rembesdalskåka, pixel-to-pixel variations in the Crocus output. This is not due to variations of atmospheric forcing (which has a 1 km grid spacing compared the 100 m grid spacing of the WRF-Hydro/Glacier simulations) or blowing snow. We suspect the pixel-to-pixel variations arise from small vertical resizing errors of the very thick glacier layers, for where we relaxed some of the test requirements when resizing. This does not change the conclusions in this paper, and work is in progress to address this issue.

Moved (insertion) [1]
much better with the observations than the Noah-MP albedo in the summertime when the snow melts and the ice becomes apparent due to the assumption of glacier land surface category in Noah-MP. 435 In regards to the observed peak flow at the beginning of melt season in May, the model does not predict such a peak flow.

Discharge
These peaks mark the onset of the hydrologic season in the rivers. They occurred almost at the same time in the two catchments 445 in two consecutive years under similar preconditions. Photographs taken on site prior the event indicate that water was flowing over the snow pack and was carving down to the river bed during the event. Probably some part of the peak is due to a pulsed meltwater flux and/or the associated pressure build up to that time which is subsequently released. Therefore, when evaluating the accumulated discharge, we start the accumulated period directly after the large spike in the observations. The accumulated discharge for 2017 shows that the observations are slightly higher than the discharge simulated with WRF-Hydro/Glacier (with 450 Crocus). They still follow closely, but WRF-Hydro/Glacier flattens out in the fall earlier compared to observations. This is likely due to lack of using the baseflow/groundwater module in these specific WRF-Hydro/Glacier simulations, which could add some water to the surface streamflow. The WRF-Hydro simulations with only Noah-MP are consistently much lower compared to both observations and WRF-Hydro/Glacier at Middalselvi. For 2018, simulations with both Noah-MP over the glacier and Crocus over the glacier follow the observations reasonably well, however they both overpredict the discharge in 455 the end of May and early June. Around July 30, the Crocus simulations increase slightly more than the observations while the Noah-MP simulations reduce discharge considerably at the end of August. The Crocus simulations still have discharge comparable with the observations until September. As shown in Figure 8j, Crocus has a larger negative summer mass balance compared to the observations which, for which the large discharge in early June in both Noah-MP and Crocus is likely connected to. Thus, even though Noah-MP compares well with observations at the end of the melt season in 2018, this is most 460 likely due to too much melt in the beginning of the melt season and not the skill of the Noah-MP simulations.
One of the reasons for the lower discharge in Middalselvi with the Noah-MP compared to Crocus is likely the lack of glacier runoff once the seasonal snowpack has melted. The Middalselvi is about 60% glacierized, while Finseelvi is only 14.7% glacierized. As can be seen in the discharge data for Finseelvi (Figure 17), the Noah-MP compares better with the observations included lateral movements of the glacier or any treatment of blowing snow in WRF-Hydro/Glacier. The impacts of including these additional physical processes may afford additional improvements in annual snowpack and snowmelt representation, 525 particularly with respect to their spatial distribution. Finally, the forcing at 1 km does not account for any topographic variations in the 100 m domain, thus snowpack evolution at 100 m scale is not included. More development and testing of WRF-Hydro/Glacier are certainly required and are currently underway. In particular the system should be applied and assessed over other glaciated regions of the world. It is a tool that can contribute to better understanding of future hydroclimate; especially in light of the continued widespread glacier retreat and decline in mass balance, which has profound implications for future 530 water resources in many parts of the world.
Tables 710 Table 1: WRF and WRF/Hydro model configuration

Model Features Configuration
Horizontal