The complex interaction between the atmospheric boundary layer and the heterogeneous land surface is typically not resolved in numerical models approximating the turbulent heat exchange processes. In this study, we consider the effect of the land surface heterogeneity on the spatial variability of near-surface air temperature fields and on snow melt processes. For this purpose we calculated turbulent sensible heat fluxes and daily snow depth depletion rates with the physics-based surface process model Alpine3D. To account for the effect of a heterogeneous land surface (such as patchy snow covers) on the local energy balance over snow, Alpine3D is driven by two-dimensional atmospheric fields of air temperature and wind velocity, generated with the non-hydrostatic atmospheric model Advanced Regional Prediction System. The atmospheric model is initialized with a set of snow distributions [snow cover fraction (SCF) and number of snow patches] and atmospheric conditions (wind velocities) for an idealized flat test site. Numerical results show that the feedback of the heterogeneity of the land surface (snow, no snow) on the near-surface variability of the atmospheric fields result in a significant increase in the mean air temperature ΔT_a = 1.8 K (3.7 and 4.9 K) as the SCF is decreased from a continuous snow cover to 55% (25 and 5%). Mean air temperatures over snow heavily increase with increasing initial wind velocities and weakly increase with an increasing number of snow patches. Surface turbulent sensible heat fluxes and daily snow depth depletion rates are strongly correlated to mean air temperatures, leading to 22–40% larger daily snow depth depletion rates for patchy snow covers. Numerical results from the idealized test site are compared with a test site in complex terrain. As slope-induced atmospheric processes (such as the development of katabatic flows) modify turbulent sensible heat fluxes, the variation of the surface energy balance is larger in complex terrain than for an idealized flat test site.