Arctic permafrost thaw holds the potential to drastically alter the Earth’s surface in Northern high latitudes. We utilize high-resolution Large Eddy Simulations to investigate the impact of the changing surfaces onto the neutrally stratified Atmospheric Boundary Layer (ABL). A stochastic surface model based on Gaussian Random Fields modeling typical permafrost landscapes is established in terms of two land cover classes: grass land and open water bodies, which exhibit different surface roughness length and surface sensible heat flux. A set of experiments is conducted where two parameters, the lake areal fraction and the surface correlation length, are varied to study the sensitivity of the boundary layer with respect to surface heterogeneity. Our key findings from the simulations are the following: The lake areal fraction has a substantial impact on the aggregated sensible heat flux at the blending height. The larger the lake areal fraction, the smaller the sensible heat flux. This result gives rise to a potential feedback mechanism. When the Arctic dries due to climate heating, the interaction with the ABL may accelerate permafrost thaw. Furthermore, the blending height shows significant dependency on the correlation length of the surface features. A longer surface correlation length causes an increased blending height. This finding is of relevance for land surface models concerned with Arctic permafrost as they usually do not consider a heterogeneity metric comparable to the surface correlation length.