Turbulent exchange of heat and moisture plays an important role in snow cover dynamics in mountain regions and governs the boundary layer dynamics. While these processes are subject to great spatiotemporal variability, particularly in complex terrain, virtually all measurements of heat, moisture and momentum fluxes are point observations. To quantify the spatial variability, and assess the representativeness of the observations, numerical modelling of the atmosphere and surface is a useful tool. Still, there is substantial uncertainty in the accuracy of how surface models capture this spatial variability, particularly in complex terrain with large spatial variability on small scales. These uncertainties can be partially attributed to (1) the use of Monin-Obukhov similarity theory (MOST) which has limitations in complex terrain due to the role of advection and (2) the errors in representing near-surface atmospheric gradients in the simulations. In this study, we analyse sources of errors in representing energy exchange over snow in mountain regions and look specifically at the spatiotemporal variability during different meteorological events in the region of Davos, Switzerland. To verify common modelling approaches with observations, we use model predictions of turbulent fluxes from CRYOWRF, the atmospheric model WRF coupled to the surface model SNOWPACK. The fluxes at different resolutions are compared to turbulent fluxes measured using the Eddy Covariance method (EC) and calculated with MOST. This model validation is done for different meteorological events representative of the local climate. Preliminary results indicate that the fluxes are highly spatially variable, being an order of magnitude higher on the leeside than on the windward side of a mountain ridge. This indicates that local heat fluxes are not representative of the whole mountain area, which has implications for the calculation of snow melt, sublimation and accumulation across mountainous terrain. The resolution of the model also plays a large role in representing the fluxes as the modelled fluxes differ greatly depending on the resolution. Our results quantify to what extent snow-atmosphere interactions and their spatial variability are correctly represented in state-of-the-art numerical weather and snow models.