On the performance of the snow model Crocus driven by in situ and reanalysis data at Villum Research Station in northeast Greenland

Reliable and detailed snow data are limited in the Arctic. We aim at overcoming this issue by addressing two questions: (1) Can the reanalysis ERA5 replace limited in situ measurements in high latitudes to drive snow models? (2) Can the Alpine model Crocus simulate reliably Arctic snow depth and stratigraphy? We compare atmospheric in situ measurements and ERA5 reanalysis and evaluate simulated and measured snow depth, density and specific surface area (SSA) in northeast Greenland (October 2014–October 2018). To account for differences between Alpine and Arctic region, we introduce a new 15 parametrisation for the density of new snow. Our results show a good agreement between in situ and ERA5 atmospheric variables except for precipitation, wind speed and direction. ERA5’s resolution is too coarse to resolve the topography in the study area adequately, leading presumably to the detected biases. Nevertheless, measured snow depth agrees better with ERA5 forced simulations than forced with in situ measurements. 20 Crocus can simulate satisfactory the evolution of snow depth, but simulations of SSA and density profiles for both forcings are biased compared to field measurements. Adjusting the new snow density parametrisation leads to improvements in the simulated snow stratigraphy. In conclusion, ERA5 can be used instead of in situ measurements to force snow models but the use of Crocus in the Arctic is affected by limitations likely due to the missing vertical water vapour transport and snow redistribution during strong winds. These limitations strongly affect the accuracy of the vertical profiles of physical snow 25 properties.

The overall objective of this paper is to evaluate the performance of the snow model Crocus driven by in situ and reanalysis data for Greenland. In particular, we aim to answer two main questions: (1) Can the reanalysis ERA5 replace in situ measurements where in situ data is limited? (2) Can we apply the Alpine snow model Crocus in the Arctic to obtain reliable snow depth evolution and other physical snow properties as, e.g. the vertical profiles of snow density and SSA. 70 For this study we use in situ observations from Greenland at Villum Research Station, where both meteorological data and snow depth were measured. In addition, during a campaign in spring 2018 numerous snow parameters, as snow depth, snow density and SSA were measured.

Methods
We investigated the potential of using the atmospheric reanalysis ERA5 to force the snow model Crocus by comparing the 75 needed near surface data to meteorological in situ measurements at Villum Research Station. We further evaluated the performance of Crocus to simulate physical snow properties at this high Arctic location.
The study area is located in the surroundings of the atmospheric measurement mast of Villum Research Station (81°34' N, 16°38' W, 37 m above sea level) (hereafter Villum Research Station) in northeast Greenland (Fig. 1). The Wandel Sea in the north, a local ice cap called Flade Isblink in the south and east and fjords in the west surround the low land, which consists of 80 continuous permafrost with an active layer of 20 cm in maximum. The observed annual mean temperature is -21°C and the annual precipitation is 188 mm (Rasch et al., 2016).

The snow model Crocus
We used the multilayer snow model Crocus (Brun et al., 1992;Vionnet et al., 2012) in single-column mode and evaluated its performance in the Arctic on bare soil. The snow model is embedded in the surface scheme SURFEX and is coupled with the 90 soil model ISBA-DF (Boone et al., 2000). We use SURFEX version 8.1. The model can deal with up to 50 vertical dynamical snow layers. Layer thickness, heat content, density and age characterise each snow layer. The number of the snow layers is dynamical, i.e. layers can run empty (zero thickness).
Crocus describes the evolution of the snowpack by taking the energy-and mass balance of the snowpack into account.
Implemented processes are surface melting, internal melting and refreezing, compaction, snow metamorphism, providing a 95 tracer for snow age, near surface densification and enhanced sublimation as well as fragmentation and compaction due to snowdrift and solar absorption in the snowpack. Snow compaction and microstructural changes due to wind drift in Crocus https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. occur when wind speed exceeds a transport threshold depending on snow properties. Commonly, this threshold is at 6 m s -1 (Vionnet et al., 2013). These wind drift effects are passed on to the underlying layers with an exponential decay until a layer's transport threshold is lower than the wind speed. Snow redistribution is not taken into account (Vionnet et al., 2013). 100 Crocus is driven by surface fluxes calculated from the atmospheric forcing variables air temperature, specific humidity, wind speed, incoming direct and diffuse shortwave and longwave radiation, rain-and snowfall rate and surface air pressure. In addition to the atmospheric forcing, the model uses the terrain information aspect, slope and altitude as well as optionally wet and dry deposition coefficients of black carbon or other light absorbing impurities (Tuzet et al., 2017). A detailed description of the model can be found in Vionnet et al. (2012). 105

ERA5 reanalysis data
The global atmospheric reanalysis ERA5 is available in the Climate Data Store of Copernicus on a regular latitude-longitude grid at 0.25° x 0.25° resolution, with atmospheric parameters on 37 pressure levels and has a temporal resolution of one hour (Hersbach et al., 2020). For Villum Research Station the grid cell size is 5 km along the longitude and 31 km along the latitude. 110 The output of the simulations driven by ERA5 are representable on spatial scales of at least one grid cell size, i.e. on 5 km x 31 km. Delhasse et al. (2020a) found a good performance for near-surface variables over the Greenland Ice Sheet as well as other studies in snow-covered locations worldwide (Albergel et al., 2018;Wang et al., 2019;Urraca et al., 2018).
Instead of interpolating spatially to the location of Villum Research Station, we took the nearest available grid cell (mid-point at 81.5° N 16.75° W, 185 m above sea level), which is about 10 km south from Villum Research Station. This is done to avoid 115 spatial interpolation of the needed atmospheric variables, which might lead to destruct their physical consistency due to nonlinearities in the underlying atmospheric process formulation in the model. We ran the model with atmospheric input from four additional nearest grid cells to Villum Research Station located on land ( Fig. 1) to account for biases resulting from taking the nearest neighbour. The results of the additional runs were taking as a metric for the co-location error. i.e. were deemed to give information of the representativeness of the ERA5 data (see Loew et al., 2017 for an introduction of the terminology 120 used).
Several data conversions were needed to drive Crocus: (i) ERA5 provides the variables snowfall and precipitation. However, model tests showed a better agreement between simulated and observed snow depths when splitting ERA5's total precipitation rate into liquid and solid precipitation rates after Jennings et al. (2018). Therefore, precipitation at temperatures of 1°C or warmer is handled as liquid precipitation and the residual as snowfall. (ii) Incident shortwave radiation was split into diffuse 125 and direct shortwave radiation using the zenith angle depending on the solar declination, the local time and latitude. (iii) Wind speed and wind direction were calculated from the wind vector and (iv) specific humidity from surface air pressure, temperature and dew point temperature according to Willett et al. (2008).

Meteorological observations
Villum Research Station provided meteorological data obtained by an automatic weather station installed in 2014 (Villum 130 Reserach Station, 2021). Table 1 gives an overview about the installed meteorological sensors and their accuracy. To force Crocus, we resampled the data to obtain hourly mean and accumulated data, respectively, for our study period from 27 November 2015 to 08 August 2018. We used ERA5 data to fill in missing data for all variables, except for precipitation, as especially the amount and timing of strong precipitation events were too diverse between ERA5 and in situ measurements. In addition, an overestimation of precipitation was introduced when filling missing precipitation data with ERA5 precipitation, most likely because of the different spatial scales of precipitation in ERA5 and the in situ data (the latter containing much smaller scales). Therefore, we 140 set missing precipitation data to zero deliberately introducing a potential underestimation of the accumulated precipitation.
Most of the surface input variables required by Crocus were measured directly. However, some conversions had to be carried out to get all required forcing variables in the appropriate form: (i) Precipitation was split in liquid and solid phase after Jennings et al. (2018), as it was done for the ERA5 precipitation (Sect. 2.3.1). (ii) The zenith angle depending on the solar declination, the local time and latitude is used to separate incident shortwave radiation into diffuse and direct shortwave 145 radiation. (iii) The Magnus-equation (Magnus, 1844) is applied to calculate specific humidity from the measured surface air pressure, air temperature, and relative humidity.

Validation data
We evaluated the model performance by comparing simulated and measured snow depth, snow density and SSA ( Vertical profiles of snow density and SSA were retrieved from a SnowMicroPen (SMP), a snow penetrometer measuring the 155 bonding force between snow grains with a vertical resolution of 1.25 mm (Schneebeli et al., 1999). While the density retrieved from SMP measurements was calculated following Calonne et al. (2020), Proksch et al. (2015) and King et al. (2020), SSA was obtained based on Calonne et al. (2020) and Proksch et al. (2015), as King et al. (2020) did not provide any SSA parametrisation. After comparing the results for all parametrisations, we used results after Calonne et al. (2020) (see Sect. 4.2.2 for discussion). Due to thick ice layers within the snowpack, the SMP could not penetrate the full vertical profile and snow 160 depth measurements were not possible. To keep, a comparable snow depth between SMP measurements and simulation results, we used the latter as reference height (maximum 2.4 m during PAMARCMiP campaign) also for the SMP measurements.
In addition, on 03 April 2018, snow depth was measured with a ruler, and every 10 cm of the profile, snow density with a density cutter (60x30x56 mm) and SSA with an IceCube3. The optical system IceCube measures the hemispherical infrared reflectance of snow and converts the reflectance in SSA (Zuanon, 2013). 165

Adaptation of Crocus: New snow density parametrisation
The original parametrisation of the density of new (freshly fallen) snow, ρnew within Crocus (hereafter V12) formulated for 170 the European Alps is: where Ta is the air temperature, Tfus is the melting point temperature for water, U is the wind speed and aρ is 109 kg m -3 , bρ is 175 6 kg m -3 K -1 , and cρ is 26 m -7/2 s -1/2 (Vionnet et al., 2012). Figure 2a shows the dependency of the new snow density using V12 on wind speed and air temperature. With this linear dependence on temperature the new snow density becomes negative for temperatures below -30°C. To avoid negative and unrealistic small new snow densities, the minimum new snow density is originally set to 50 kg m -3 . However, this parametrisation was deduced from measurements at Col de Port, French Alps, where https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License.
temperatures below -20°C occur seldom (Vionnet et al., 2012). For the Arctic, where temperatures below -20°C are frequent, 180 this parametrisation leads to numerous new snow densities of the minimum density (50 kg m -3 ) deemed unrealistic for high Arctic conditions. The lowest measured snow density during PAMARCMiP, for instance, was 117 kg m -3 . Also other studies found too low new snow densities for the Arctic simulated with Crocus (Essery et al., 2016;Sauter and Obleitner, 2015). Therefore, we aimed to find a new snow density parametrisation suitable for the cold temperatures in the Arctic preferable without an "unphysical" ad hoc cut-off density at low temperatures. We compared the performance of Crocus for alternative new snow density parametrisations from Anderson (1976), Liston et al. (2007), Sauter and Obleitner (2015) and Van Based on that, we use for this study a combination of Liston et al. (2007) and Van Kampenhout et al. (2017) as it represents best our observations and needs. 195 The new snow density parameterisation by Liston et al. (2007) is defined for ≥ 258.16 K and is given as: for wind speed below 5 m s -1 , where is the wet-bulb air temperature in Kelvin. For wind speed ≥ 5 m s -1 the new snow density is defined as: with 205 where 1 is the density offset of 25.0 kg m -3 for a 5.0 m s -1 wind, 2 is the maximum density increase of 250.0 kg m -3 caused by wind and 3 sets the progression from low to high wind speed and is 0.2 m s -1 . Wt is the terrain-modified wind speed at 2 m 210 above the surface in m s -1 (Liston et al., 2007).
To simplify Eq. (4) we used directly measured wind speed for in situ and given ERA5 wind speed without any terrain modifications. Wind speed was measured at 9 m and ERA5 wind speed was calculated for 10 m height. We tested the log law and the power law approach to downscale wind speed to 2 m with almost identical results in snow density and snow SSA for both approaches. Therefore, we present only results for 9 m (in situ) and 10 m (ERA5) wind speeds in this study. 215 For temperatures below -15°C, we extend the parameterisation by the new snow density parametrisation from Van Kampenhout et al. (2017) but keep the wind depended part from Liston et al. (2007): where is the near surface air temperature and is the freezing temperature of water 273.15 K.
The parametrisation for temperatures below -15°C from Van Kampenhout et al. (2017) does not consider any dependence on 225 humidity. However, the influence of humidity below -15°C should be negligible, as humidity at temperatures below -15°C is low. Density inversions at temperatures below -15°C as observed by Van Kampenhout et al. (2017) are taken into account with this parametrisation. Figure 2 visualise the different new snow densities obtained with the original parametrisation and our adapted parametrisation (hereafter K21).

Model simulations setup 230
An overview of all model simulations is listed in Table 3. We forced several model simulations with ERA5 data over the period from 2010 to August 2020. Due to a constant negative offset throughout the years in the reanalysis surface air pressure, https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License.
probably caused by orographic effects (see Sect. 3.1.1), we adjusted ERA5 surface air pressure by adding the mean difference between in situ air pressure and ERA5 air pressure. To reduce spin-up effects, introduced through the initial conditions of the soil, we used the initial conditions on 31 December 2019 after running one ERA5 forced simulation from 2010 to December 235 2019. This second pass of the forcing period simulation is called ERA5 control run (ERA5-CTRL). ERA5-NSD used the same forcing as ERA5-CTRL but with the introduced new snow density parametrisation K21. Simulations to analyse the impact of simulated snow depth to forcing variables (ERA5-sens) (see Sect. 2.6) used ERA5-CTRL forcing but with a disturbance on each individual forcing variable.
The in situ control simulation (Insitu-CTRL) is forced with the in situ measurements from Villum Research Station from 240 27 November 2015 to 08 August 2018, except from erroneous measured longwave radiation, which was replaced by the corresponding ERA5 variable for the whole study period. To reduce spin-up effects for the in situ simulation, we used the archived initial conditions from 26 November 2015 of ERA5-CTRL, as initial conditions for the in situ simulation. The in situ simulation using the introduced new snow density parametrisation K21 is called Insitu-NSD.
For all simulations we set the absorption enhancement parameter B0 to 1.6 in line with other studies (Tuzet et al., 2020;Libois 245 et al., 2013). As we had no measurements for wet and dry black carbon deposition rates, we used results from Tuzet et al.
(2020) for the French Alps as indication, but assumed lower rates for northeast Greenland. Therefore, we set wet and dry constant deposition rates to 6.67 x 10 -11 g m -2 s -1 and 3.33 x 10 -11 g m -2 s -1 , respectively.

Sensitivity simulations
To understand the difference between the simulated snow depth obtained from Insitu-CTRL and ERA5-CTRL, respectively, we used a linear approach of partial disturbance. With this approach we were able to figure out the contribution of every forcing variable to the difference in the simulated snow depth.
We used ERA5-CTRL as baseline simulation because of its consistency and completeness of all variables for our sensitivity 255 survey. From the ERA5 data, we calculated the daily standard deviation from 2010 to August 2020 for every forcing variable.
We added one-tens of the standard deviation to the individual variable as typical disturbance di while all other variables were https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. unchanged in each simulation of ERA5-sens. Then, we determined the difference between snow depth values retrieved from ERA5-sens and ERA5-CTRL for every simulation with one disturbed forcing variable. The sensitivity e(xi) for every forcing variable is defined as: 260 with di being one-tens of the standard deviation from the ERA5 ith forcing variable.
We estimated the influence of each forcing variable on the snow depth difference between Insitu-CTRL and ERA5-CTRL as: 265 with m(xi') -m(xi) being the mean bias of in situ and ERA5 simulated snow depth and with xi' -xi being the mean difference between the ith in situ and the corresponding ERA5 forcing variable. For a perfect linear model, the sum over all forcing 270 variables xi of m(xi')-m(xi) would be equal to the snow depth difference between Insitu-CTRL and ERA5-CTRL. However, because of the non-linear model the sum is an approximation of the "real" difference of both simulations but allows us to identify the main reasons of the difference between Insitu-CTRL and ERA5-CTRL.

Comparison of reanalysis data and in situ measurements from Villum Research Station
Villum Research Station provided meteorological measurements that we compared with ERA5 data to evaluate the performance of the reanalysis for Greenland. We removed the annual cycle from the time series by subtracting the observed monthly mean "climatology" from 27 November 2015 to 08 August 2018 from the observed in situ data and the ERA5 reanalysis, respectively. By analysing these anomalies, we identified the mean biases and correlations between the in situ data 280 and the reanalysis. Figure 3 shows the comparison of daily anomalies of the atmospheric forcing variables from ERA5 with the in situ data.
Observed direct shortwave radiation, diffuse shortwave radiation, specific humidity and air temperature agreed well over the period from 27 November 2015 to 08 August 2018. Mean biases were small and the Pearson correlation coefficients were in the range of 0.84 to 0.91. ERA5 underestimated the air temperature by 1.0°C and the specific humidity by 0.0002 kg kg -1 . 285 Longwave radiation showed a poor correlation (r: 0.28) and a high root-mean-square difference (RMSD: 60.5 W m -2 ) explained by a tilt of one of the sensors used to obtain longwave radiation (per. Comm. with Andreas Massling).
Atmospheric surface pressure correlated well with in situ surface pressure (r: 0.99) but showed a systematic offset of 19 hPa.
We attribute that bias to the altitude difference of the ERA5 nearest grid cell to Villum Research Station of about 150 m due to the relatively coarse resolution of ERA5. ERA5 surface air pressure of the four grid cells closest to Villum Research Station 290 showed a large spread (mean surface air pressure difference -9 hPa, standard deviation: 13 hPa). Grid cells further north had a lower elevation than the grid cells at 81.5° N (Fig. 1).
To account for the offset, we corrected ERA5 surface pressure by the mean difference between in situ and ERA5 over the period. As we used surface pressure to calculate the specific humidity from ERA5's dew point temperature, changes in pressure influenced the specific humidity. However, effects were small. 295 Rain-and snowfall showed mean biases of 0.1 mm day -1 and 0.6 mm day -1 w.e. (water equivalent) and poor correlations of 0.15 and 0.13, respectively. Additionally, precipitation measurements were incompletefor about 35 % of the time no data were available. Observed and ERA5 precipitation differed strongly not only in the amount but also in the timing of extreme events. In situ precipitation was generally higher than ERA5 precipitation. ERA5 rainfall tended to be delayed and weaker than in situ rainfall.  ERA5 and measured wind direction show large deviations. Measured preferred wind direction was from the south-west while ERA5's preferred wind direction was from the south (Fig.4). Also, wind direction for wind speeds higher than 5 m s -1 , which https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. potentially lead to snow redistribution due to wind drift, differed considerable. The mean wind speed from ERA5 over the 310 whole study period was slightly higher than measured (0.2 m s -1 ). However, ERA5 underestimated the occurrence of wind speeds higher than 5 m s -1 .

Atmospheric conditions during and shortly before the PAMARCMiP campaign
In the second half of February 2018 there was a remarkable period of exceptional weather conditions lasting about two weeks ( Fig. 5), where even air temperatures above freezing temperature were present in hourly data (Moore et al., 2018;Ludwig et al., 2019). The mean daily air temperatures between 17 Feb 2018 and 27 Feb 2018 were about 20°C warmer than in the weeks 320 before and after. During PAMARCMiP itself, air temperature was about -20°C.
The warm event was accompanied by remarkable strong winds, which were above the wind drift mark of 5 m s -1 . Hourly in situ wind speeds exceeded 20 m s -1 . During the time of the snowpit survey observed wind speeds were slightly below 5 m s -1 with some hours above 5 m s -1 at the beginning of the campaign. During the warm and strong wind event, wind direction was mainly south to southeast, shifting towards southwest afterwards. 325 ERA5 underestimated snowfall during the first warm event by more than five times. Due to the short periods of air temperature above the freezing temperature there was also some rainfall measured. ERA5 rainfall during that time was close to zero at Apart from the first warm event the amount of precipitation in both datasets were similar, but the timing of precipitation events 330 differed clearly. During PAMARCMiP, in situ measurements showed a snowfall event over several days, which is not present in ERA5 instead ERA5 showed about a week before the campaign snowfall. Field observations during PAMARCMiP confirmed strong snowfall events together with strong winds and snow drift.  (Fig. 6). In these years, snow depth was within the co-location error of the ERA5 data until early spring. The snow depth measurements for PAMARCMiP showed distinctive differences in snow depth even over small distances. We measured snow depths from 60 to 117 cm with a standard deviation of 21.6 cm from the snow depth measurements during the detailed snowpits on 03 April 2018. 345 Melting started too early in the years 2015/16 and 2016/17 in ERA5-CTRL compared to the measurements (Fig. 6). Measured and simulated snow depth at the beginning and end of the snow season differed slightly, with the overall timing of the buildup of snow cover matching better than the precise timing of no snow cover. Date and amount of the maximum snow depth showed some moderate differences. Over the whole study period, the mean bias between ERA5-CTRL snow depth and the observation was only -0.17 m with a 360 daily RMSD of moderate 0.29 m (Fig. 7a). In 2017/18, Insitu-CTRL overestimated the maximum snow depth by almost a meter and showed a later increase in snow depth than measured and ERA5-CTRL. From late February to early March, there was a distinct increase in snow depth of about 375 1.5 m within a few days. Even if the mean bias between Insitu-CTRL and observed snow depth was low (0.06 m) and RMSD was 0.39 m over the whole study period, explained variance, which is the proportion of variability in the data explained by the model, was only 0.23 (Fig. 7b). For the same period ERA5-CTRL snow depth showed a much higher explained variance of 0.73.
As our sensitivity survey showed (Table 4), especially snowfall and specific humidity but also air temperature caused 380 differences between simulated snow depth in Insitu-CTRL and ERA5-CTRL. The remaining forcing variables had clearly a lower impact on snow depth differences. The large response to the differences in specific humidity turned out to be an artefact of the linearisation. The experiment was rerun with a disturbance of 0.0002 kg kg -1 and resulted in a difference in snow depth of 16.7 cm showing the limits of our approachthe sensitivity calculated for a disturbance of one tens of the standard deviation of the variability of the specific humidity was overestimated. Adding up the contributions of snowfall, specific humidity 385 (0.0002 kg kg -1 ) and air temperature amounted to 0.29 m and compared well with the mean difference between Insitu-CTRL and ERA5-CTRL of 0.23 m.

Snow density profiles
Comparisons between different approaches to receive density and SSA from SMP penetration force are shown in Fig. 8. While 395 density calculated after Calonne et al. (2020) and King et al. (2020) showed some similarities, densities after Proksch et al. (2015) were considerably higher in the upper part of the profiles and generally more variable expressed in larger standard deviation ranges (Fig. 8 shaded areas). SSA after Calonne et al. (2020) differed strongly from Proksch et al. (2015) and is clearly higher and more variable.  (Fig. 9a, b). Additionally, we compared the simulations to a density profile obtained with a density cutter on 03 April 2018 (Figure 10a, b). Overall, simulated and measured density profiles showed a poor agreement for the CTRL simulations. 405 Simulated snow densities, both in ERA-CTRL and Insitu-CTRL were much lower than measured with the SMP or with the density cutter. This was especially true for the upper part of the profile.
Density variations with depth within the simulated snowpack were mostly small. Simulated snow densities showed a clear densification with depth and did not capture the snow density decrease from 388 kg m -3 to 270 kg m -3 measured with the density cutter in the lower part of the snowpack (e.g. below 0.3 normalised snow depth, Fig. 10a). SMP measurements showed 410 a pronounced increase in the density from 180 kg m -3 to 300 kg m -3 over a few centimetres below the snow surface, which https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. density cutter measurements could not resolve. Both simulations were not able to capture this density increase near the surface (Fig. 9a, b).

SSA profiles
Overall, SSA measured by the SMP and simulated were in good agreement when the sample variability of the measurements is taken into account (the simulated SSA were within the range spanned by ±1 standard deviation of the sampling variability). SSA in ERA5-CTRL was much more variable on 23 March 2018 than in Insitu-CTRL (Fig. 9c, d). SSA in ERA5-CTRL 430 showed three peaks while SSA in Insitu-CTRL is similar to ERA5-CTRL but showed less variability. Measured SSA with the SMP showed a huge standard deviation, which led to some match with simulated SSA even if both simulated profiles differed by about 10 m 2 kg -1 in average.
The vertical pattern of the SSA profiles of both simulations were similar to observations taken with the IceCube on 03 April 2018 (Fig. 10). All profiles showed a decrease in SSA from the surface towards the middle part of the snow continued by an 435 almost constant course to greater depth. However, amplitudes of the decrease and the height of SSA differed considerably.
ERA5-CTRL on 03 April 2018 started at the surface with a reliable SSA of 47 m 2 kg -1 but decreased earlier and stronger than the measurements with depth (Fig. 10d). Therefore, ERA5-CTRL underestimated SSA over the whole snow profile and showed mostly values below Insitu-CTRL. SSA in ERA5-CTRL decreased from the top of the snowpack until 0.6 m, followed https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. by a more or less constant middle part with SSA lower than 10 m 2 kg -1 . SSA increased in the lower part of the profile until the 440 ground, which contradicts the measurements. SSA in In-situ-CTRL (Fig. 10e) showed an overestimation of surface SSA of about 15 m 2 kg -1 . Insitu-CTRL SSA decreased stronger than measured SSA before it became stable from 0.8 to 0.1 normalised snow depth while measured SSA showed a slight decrease over the whole profile with higher SSA than Insitu-CTRL.
Overall, ERA5-CTRL and Insitu-CTRL could capture the gross profile of measured SSA. ERA5-CTRL started with more reliable surface SSA than Insitu-CTRL. However, both simulations failed to reproduce the details in the measured SSA 445 profiles. Moreover, the SSA decreasing rates from top to bottom of the profiles were different for simulations and measurements.

Impacts of the new snow density parametrisation
With the introduced new snow density parametrisation K21, we could reduce mean bias and RMSD of snow depth in ERA5-NSD by 0.04 m and 0.07 m, respectively, for the simulated period 2015 to 2018 (Fig. 7c). Furthermore, the explained variance 450 increased by 0.11 to 0.84. The Insitu-NSD's snow depth mean bias and RMSD could be improved for snow depth by 0.04 m and 0.05 m, respectively. The explained variance increased by 0.17 to 0.40 (Fig. 7d). Figure 11 visualises the measured and simulated snow depth time series for CTRL and NSD simulations. ERA5-NSD's snow depths were slightly higher towards the end of the snow season than ERA5-CTRL's (Fig.11a). Insitu-NSD showed clearly lower snow depths than Insitu-CTRL for 2015/16 and 2017/18 (Fig. 11b). 455 Applying K21, we obtained more reliable density simulations than using the original parametrisation V12. The simulated snow densities showed a much higher variability with several distinct density peaks, which were also visible in the SMP profiles and the density cutter profiles (Fig. 9a, b and Fig. 10a, b). The densities with K21 were also considerably higher than with V12. Furthermore, the snow density not only increased with depth, but also showed more pronounced variability.  With K21, we could occasionally reach the sample variability ranges of the measurements. However, it was still not possible to simulate the sharp density increase near the snow surface observed with the SMP. Also, the model could not capture the 465 decrease in snow density towards the ground as measured with the density cutter. The underestimation of the density in the uppermost layers and the overestimation of the basal layers remained.
K21 affected the simulated SSA to a smaller extent than the snow density. The overall SSA profile of ERA5-NSD was similar to ERA5-CTRL on 23 March 2018 (Fig. 9c, d). SSA in Insitu-NSD changed mainly in the upper and in the lower part of the profile while the middle part remained: After the decrease of SSA to 20 m² kg -1 , SSA remained constant as in Insitu-CTRL, 470 but increased slightly before the decrease at the lower part of the profile. Changes due to K21 on 03 April 2018 compared to SSA obtained with the IceCube were similar to SSA obtained by the SMP on 23 March 2018 (Fig. 10c, d).

Forcing data
One of the advantages of ERA5 is its physical consistency and completeness of time series, which is not the case with the in 475 situ measurements. However, the reanalysis is only able to resolve processes, which are within the scales of motion of the model resolution or even larger (Minola et al., 2020). Precipitation on a kilometre scale can already vary significantly (Betts https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. et al., 2019). Therefore, ERA5 cannot resolve e.g. local precipitation events also visible in the atmospheric data during PAMARCMiP where ERA5 underestimates precipitation during the warm event in the second half of February 2018.
Atmospheric variables measured at Villum Research Station are point measurements representing local conditions. ERA5 on 480 the other hand represents mean values over a whole grid cell. Due to the relatively coarse resolution of ERA5, the relief is much smoother than in reality. This causes altitudinal differences and leads to different values e.g. in surface air pressure than observed. After correcting for this altitudinal differences ERA5 air pressure agrees well with measurements, which is in line with the findings of Delhasse et al. (2020a).
We find a good correspondence between measured and ERA5 air temperatures, shortwave and longwave radiation with only 485 minor differences in agreement with other studies (Betts et al., 2019;Wang et al., 2019;Delhasse et al., 2020a). Differences in observed and ERA5 surface air temperature can also result from the elevation difference between Villum Research Station and the ERA5 grid cell. In our study, the use of longwave radiation from ERA5 instead of observations for the in situ forcing enabled us to circumvent measurement errors caused by a tilted sensor.
Differences in local topography due to the relative coarse resolution of ERA5 also affect wind speed and wind direction 490 (Delhasse et al., 2020a). Higher obstacles to the southeast from Villum Research Station, not resolved correctly by ERA5, can alternate and influence wind direction and wind speed. In addition, the vertical interpolation from higher levels with negligible influence of its underlying terrain causes differences in wind speed and wind direction between the reanalysis and in situ measurements. This can be problematic in terms of snow modelling, as wind is an important variable to correctly simulate snow density due to wind compaction, enhanced sublimation and wind redistribution at higher wind speeds. 495 In our study, ERA5 overestimates mean wind speed over the whole study period while Betts et al. (2019) found an underestimation for the Canadian Prairies. Delhasse et al. (2020a) identified an underestimation over stations located mainly in the ablation area of the Greenland Ice Sheet but mostly overestimation for stations located mainly in the accumulation area (Delhasse et al., 2020b). We identify a higher wind speed underestimation for higher ERA5 wind speeds in agreement with Betts et al. (2019) and Delhasse et al. (2020a). Betts et al. (2019) identified an quasi-linearly increasing underestimation of 500 wind speed for ERA5 while Delhasse et al. (2020a) described an enhanced underestimation for wind speeds higher than 3 m s -1 .
Our analyses of the wind direction show strong differences between reanalysis and in situ data, although our simulations are not impacted because of no aspect ratio prescribed. However, the huge difference between measured wind direction and reanalysis might be of interest for others using ERA5. Especially when it comes e.g. to modelling of snow redistribution by snowdrift, wind direction plays an important role and has to be handled with care. 505 Crocus requires a separation in liquid and solid precipitation. We used the approach after Jennings et al. (2018) and defined all precipitation as solid for surface air temperature lower than 1°C. Such a simple approach does not represent the reality as the phase transition spread out over about -2°C and +4°C over land area (Dai, 2008). As snowfall is one of the major factors determining simulated snow depth and other snow properties, the separation in liquid and solid precipitation adds to the uncertainty. Another source of uncertainty is the conversion from relative humidity and dew point into specific humidity using 510 https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License.
surface air pressure and air temperature. Due to different available variables (ERA5: dew point, in situ: relative humidity), we had to choose different approaches to obtain specific humidity.
Furthermore, large measurement gaps in total precipitation (35.2 % missing values) are problematic. We introduce an underestimation of precipitation by setting these gaps to zero precipitation, which leads to underestimation of in situ precipitation. However, in situ simulations clearly overestimate the snow depth although missing precipitation values are set 515 to zero, i.e. our simulation suggest that the in situ precipitation is not consistent to the measured snow depth. False precipitation measurements can be caused by snowdrift, which is falsely detected as precipitation.
ERA5 is able to catch the heat waves and strong winds in the second half of February 2018. Most variables match well during this period apart from precipitation and some minor deviations in wind direction and hourly wind speed during peak times, where hourly ERA5 wind speed underestimates. In agreement with Delhasse et al. (2020a), Betts et al. (2019) and Wang et al. 520 (2019) and despite differences due to topographic and vertical resolution of ERA5, atmospheric variables used in Crocus agree well with in situ measurements except during local precipitation events.

Snow depth
Snow depth varies considerably even on small spatial scales. With the SA50A-sensor at Villum Research Station, we can only 525 measure local snow depth. On 03 April 2018, when the detailed snow pit measurements were performed, snow depth measurements varied in the vicinity of Villum Research Station between 60 and 117 cm with a standard deviation of 21.6 cm. This might be caused by several small depressions and bumps that influenced local snow depth. Also other studies found differences of 0.45 m to more than 1 m in snow accumulation over the course of a year within a radius of several kilometres (Domine et al., 2016a;Pedersen et al., 2016). Reasons for this include snow redistribution by wind, microtopographic relief 530 and vegetation (Barrere et al., 2017;Liston and Sturm, 2002). Surprisingly ERA5 simulated snow depth clearly outperforms in situ forced snow depth. Taking into account the considerable spatial snow depth variability, simulated ERA5 snow depth, which represents a mean of the grid cell, agrees unexpectedly well with the point measurements at Villum Research Station. The simulated in situ snow depth shows larger deviations from the measured snow depth than the ERA5 simulation. The simulated in situ snow depth considerably overestimates. We did not 535 expect that result as atmospheric forcing is measured nearby and we treated missing precipitation measurements (~35.2 % missing values) as no precipitation. Hence, we expected an underestimation of snow depth rather than an overestimation.
Reasons for the overestimation could be failures of the TPS-3100-sensor of Villum Research Station (precipitation data is not quality controlled) or massive influence of snow redistribution due to wind drift.
We can show that Crocus is able to simulate reliable snow depth evolutions for the Arctic when uncertainties of the atmospheric 540 forcing are negligible. Other studies have come to the same conclusion when adjusting the model or the forcing data (Luijting et al., 2018;Jacobi et al., 2010;Sauter and Obleitner, 2015;Barrere et al., 2017). Nevertheless, Luijting et al. (2018) found https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. accumulated over-and underestimations of snow depth of about 1 m over a season. This finding confirms our observations where e.g. ERA5 considerably underestimates snow accumulation during the warm air intrusion in spring 2018 and other precipitation events, which accumulates to a strong underestimated snow depth throughout the season. In our study, also 545 deviations in simulated and measured snow densities e.g. due to an underestimation of new snow densities, lead to differences in snow depth between simulations and measurements. This correspondence with findings from Sauter and Obleitner (2015) and Essery et al. (2016).
We find highest impacts of simulated snow depth on uncertainties in snowfall, specific humidity and air temperature. Sauter and Obleitner (2015) showed that uncertainties in forcing can lead to more than 3 m deviation in snow depth after one year. 550 This again highlights how problematic biases due to missing data, snowdrift and riming on instruments and measurement uncertainties in the measured data or ERA5 data are.
Furthermore, snow redistribution due to strong winds is not parametrised in the model and could be therefore another reason for deviations between observed and simulated snow depth. In the Arctic, snow redistribution is very important. Therefore, adding this process to Crocus could lead to substantial improvements in simulated snow depth (Barrere et al., 2017;Luijting 555 et al., 2018;Libois et al., 2014). The influence on density and increased sublimation due to snowdrift, already parametrised in the model, lead to distinctively lower snow depths in our simulations. However, the quality of the parametrisation of the increased sublimation during snowdrift requires more detailed consideration in future studies.
Crocus reproduces well the observed onset of snow accumulation in our study, while melting starts earlier and is faster than observed. The discrepancies in the timing and development of thawing are a widely observed problem in the application of 560 Crocus in the Arctic (Barrere et al., 2017;Essery et al., 2016;Domine et al., 2019;Luijting et al., 2018;Jacobi et al., 2010).
Reasons for this could be a higher simulated thermal conductivity of the basal snow layer and a lower simulated thermal conductivity of the near-surface layers (Barrere et al., 2017;Domine et al., 2019) as well as the lack of parametrisation of snow redistribution by wind (Luijting et al., 2018).

Snow density and SSA profiles 565
A typical Arctic snow profile shows a basal depth hoar layer with low density (150-200 kg m -3 ) and high density at the top wind slab snow layers (exceeding 400 kg m -3 ) (Barrere et al., 2017). In our SMP measurements, a density peak is visible near the surface. In the density cutter measurements we can observe a density decrease towards the basal layers. In addition, observations during fieldwork also indicate a bottom depth hoar layer in our profiles. However, due to the lack of detailed determination of stratigraphy, we cannot fully confirm the typical snow stratification from our observations. 570 While Essery et al. (2016) found a strong correlation between simulated and observed snow density profiles (r=0.74), we find considerable differences in both CTRL simulations. We see a constant increase in density with depth. This is typical for Alpine snow, where dense basal layers are common (Domine et al., 2019). Simulations overestimate density for near surface snow layers and underestimate density for basal snow layers. This correspondents with results from other studies (Gascon et al., determining the snow stratigraphy in the Arctic is missing in Crocus. Barrere et al. (2017) and Jacobi et al. (2010) considered the main reason for the inverted density profiles in the lack of a parametrisation of upward water vapour transport due to strong temperature gradients occurring in the Arctic. Several other model deficiencies can also lead to an incorrect simulation of the Arctic snow profile (Gascon et al., 2014;Domine et al., 2016b;Barrere et al., 2017).
We found an underestimation of new snow density with Crocus, which confirms the results from other studies in the Arctic 580 (Sauter and Obleitner, 2015;Essery et al., 2016;Lefebre et al., 2005). In the new snow density parametrisation currently used in Crocus, new snow density is set to 50 kg m -3 for temperatures below -30°C, which are common in Arctic winter, and a linear relationship between temperature and density is applied (Vionnet et al., 2012). In contrast, Meister (1985) found increasing new snow densities for temperatures below -15°C as we considered with the introduced new snow density parametrisation K21. Braking through wind impacts and crystallization at low temperatures are reasons for a higher new snow density in the 585 Arctic than in the warmer and more sheltered European Alps (Jordan et al., 1999). Appling K21 results in a denser new snow density in our study and thus higher maximum densities for the uppermost snow layers than with the original formulation.
With an for the Arctic adapted new snow density parametrisation we find that Crocus is able to reproduce the gross density profile and stratification for an Arctic snowpack in line with Sauter and Obleitner (2015).
A realistic simulation of snow density is crucial. Crocus uses a function of density to parametrise important snow properties 590 like SSA or thermal conductivity (Vionnet et al., 2012;Carmagnola et al., 2014). Incorrect simulated snow densities lead to erroneous thermal conductivities, which serve as input for the calculation of the temperature gradient. Therefore, the simulated metamorphism is erroneous (Domine et al., 2019). As density is a key variable influencing other important variables in Crocus, we assume that differences in simulated and measured densities can partly explain differences in simulated and measured SSA.
SSA is parameterised in Crocus as an inverse function of the optimal diameter (Carmagnola et al., 2014). Therefore, e.g. an 595 inverted stratification of snow density leads to a lower SSA and thus to a lower albedo simulated with Crocus (Domine et al., 2019). In reality, differences in SSA of snow profiles can also result from variations in snow depth, density and temperature (Jacobi et al., 2010;Carmagnola et al., 2014). SSA further depends on wind speed and temperature gradient (Domine et al., 2019;Carmagnola et al., 2014). Furthermore, uncertainties in the forcing data can lead to erroneous simulated densities and SSA as visible in deviations in ERA5-CTRL and Insitu-CTRL SSA on 03 April 2018 (Fig. 10). Different timing of a 600 precipitation event in measurements and ERA5 causes these differences. In agreement with Carmagnola et al. (2014) we find discrepancies between simulations of SSA and IceCube measurements. Strong temperature gradients in the Arctic not captures by Crocus (Domine et al., 2019) could partly explains differences in our simulated and measured SSA profiles.
Our density and SSA measurements have weaknesses in capturing the whole spatial variability around Villum Research Station due to limited measurement profiles. However, examined profiles were carried out in representative areas of the vicinity of 605 Villum Research Station. Due to many think ice lenses within the snowpack, the SMP could not penetrate through the entire snowpack. Therefore, the results give just a reasonable approximation for the surface layer.
Another uncertainty of SMP measurements is the calculation of snow density and SSA from penetration force. We applied three parametrisations to our data, with considerably different results (Fig. 9). For the development of the parametrisation, https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License. Proksch et al. (2015) used data from the European Alps and the Arctic, while measurements for the parametrisation after 610 Calonne et al. (2020) were from the Swiss Alps and therefore for warm snow. However, Calonne et al. (2020) developed the parametrisation for the SMP Version 4 we used, while the parametrisation after Proksch et al. (2015) was developed for SMP Version 2 (Calonne et al., 2020). King et al. (2020) found large discrepancies between measured density with a density cutter and SMP after Proksch et al. (2015) for snow on Arctic sea ice. Parametrisation for SSA after Proksch et al. (2015) was developed using the precise MicroCT while Calonne et al. (2020) used measurements from the IceCube. We found that 615 densities for Calonne et al. (2020) and King et al. (2020), whose parametrisation is for Arctic snow on sea ice with the SMP Version 4, differ only slightly. Simulated density is closest to results after King et al. (2020) and Calonne et al. (2020) and especially underestimated compared to densities after Proksch et al. (2015). Simulated SSA agrees better with results after Calonne et al. (2020) and is far off results after Proksch et al. (2015). Nevertheless, a reliable comparison between our simulations and the density and SSA derived from SMP is not possible due to uncertainties in the algorithms used to derive 620 these quantities.
Taken into account uncertainties due to the high spatial variability of snow density and SSA visible in the variety of the measurements, we observe a fair agreement of the ERA5-NSD profiles with the measurements. However, we cannot simulate the strong increase in the measured density in the upper layers. In addition, underestimations of snow density in the upper layers and overestimations in the lower layers remain. We assume that the main reason for this is the lack of parametrisation 625 of the vertical water vapour transport and the resulting mass transport from lower to upper snow layers. With further adaptions of Crocus and inclusion of this missing process relevant for the Arctic, we expect better results.

Conclusion
Working with snow in the Arctic is challenging due to limited temporal and spatial availability of measured snow data. Detailed snow models can help to overcome this problem and can assist to interpret the measurements but require carefully prepared 630 forcing data to ensure a good quality of simulated snow data. Using the snow model Crocus, we can show that already small deviations in solid precipitation, specific humidity and air temperature lead to considerable differences in simulated snow depth.
In situ atmospheric measurements in the Arctic are rare, incomplete, and have to cope with difficult measurement conditions such as riming on the instruments and strong winds, while models need complete time series to drive them. Our study 635 demonstrates that the reanalysis ERA5, except for precipitation, wind speed, and wind direction, agrees well with atmospheric measurements at Villum Research Station in northeast Greenland between October 2014 and August 2018. ERA5 is also physically consistent being an output of a numerical atmospheric model, which does not apply for the observations, explaining partly a higher agreement of simulated and measured snow depth under ERA5 forcing than with in situ forcing, a result not has been expected. 640 https://doi.org/10.5194/tc-2021-100 Preprint. Discussion started: 7 April 2021 c Author(s) 2021. CC BY 4.0 License.
Due to the relatively coarse resolution of ERA5, topography is not resolved in detail. Presumably for that reason, ERA5 cannot reproduce small-scale local precipitation events. Observed and reanalysis precipitation also do not match in time. Differences in topography influence wind direction and wind speed, thus ERA5 underestimates higher wind speeds in particular. Due to differences in elevation between the grid cell and the observation site, there are differences in air temperature and surface air pressure. 645 Overall, our study shows that ERA5 is capable to replace in situ measurements to force snow models where observations are limited, considering the limitations of in situ data and the reanalysis. Our study site is located in an area influenced by orography. We expect higher agreement of ERA5 and observations for flat areas as widely common in the Arctic such as sea ice. ERA5 will be a useful dataset for future studies in the Arctic, which need complete time series of meteorological data. Our first question raised in the introduction if ERA5 data can be used as a surrogate for missing in situ forcing data can be 650 unambiguously affirmed.
Concerning our second research question, whether the Alpine snow model Crocus can reliably simulate snow depth development as well as snow density and SSA profiles, we cannot give a conclusive answer as Crocus can simulate the rough evolution of snow depth, but not of SSA and snow density of an Arctic snowpack. Thereby the ERA5 simulations outperforms in situ simulations. Differences between Alpine and Arctic atmospheric conditions and thus snow properties and processes can 655 explain the deviations between observations and simulations in the Arctic. Crocus strongly underestimates new snow density, leading to an underestimation of snow density over the entire profile. In our study, we introduce a different new snow density parametrisation that results in more representative snow densities for Arctic snow profiles than with the original parametrisation. Using our new snow density parametrisation, Crocus still underestimates snow density in the upper layers and overestimate snow density in the basal layers. We attribute this to a lack of parametrisation of important processes in the 660 Arctic such as vertical water vapour transport and snow redistribution by strong winds.
Our study shows that Crocus has great potential for simulating snow conditions in the Arctic. The model can contribute to complement temporally and spatially limited observed snow depth measurements through representative simulations.
Simulated and observed snow densities and SSA show deviations of which the user needs to be aware. Having that in mind, Crocus can give added value to the evolution and the range of the prevailing density and SSA. This is a clear improvement 665 over the otherwise frequent use of fixed density values for all layers and values from climatologies for e.g. ice-ocean models.
We hope this manuscript is a step forward to achieve these goals.