Simulating the influence of heterogeneous surfaces on atmospheric flow using mesoscale models (MSM) remains a challenging task, as the resolution of these models usually prohibits resolving important scales of surface heterogeneity. However, surface heterogeneity impacts fluxes of momentum, heat, or moisture, which act as lower boundary conditions for MSM. Even though several approaches for representing subgrid-scale heterogeneities in MSM exist, many of these approaches rely on Monin–Obukhov similarity theory, preventing those models from resolving all scales of surface heterogeneity. To improve upon these residual heterogeneity scales, a novel heterogeneity parameterization is derived by linking the heterogeneous covariance function in spectral space to an associated homogeneous one. This covariance function approach is subsequently used to derive a parameterization of the aerodynamic resistance to heat transfer of the surface layer. Here, the effect of surface heterogeneity enters as a factor applied to the stability correction functions of the bulk similarity approach. To perform a first comparison of the covariance function approach against the conventional bulk similarity and tile approaches, large-eddy simulations (LESs) of distinct surface heterogeneities are conducted. The aerodynamic resistances from these three parameterizations are subsequently tested against the LES reference by resolving the surface heterogeneities with six different test-MSM grids of varying cell dimension. The results of these comparisons show that the covariance function approach proposed here yields the smallest deviations from the LES reference. In addition, the smallest deviation of the covariance function approach to the reference is observed for the LES with the largest surface heterogeneity, which illustrates the advantage of this novel parameterization.