The earth’s surface is characterized by small-scale heterogeneity attributable to variability in land cover, soil characteristics and orography. In atmospheric models, this small-scale variability can be partially accounted for by the so-called mosaic approach, i.e., by computing the land-surface processes on a grid with an explicit higher horizontal resolution than the atmosphere. The mosaic approach does, however, not account for the subgrid-scale variability in the screen-level atmospheric parameters, part of which might be related to land-surface heterogeneity itself. In this study, simulations with the numerical weather prediction model COSMO are shown, employing the mosaic approach together with a spatial disaggregation of the atmospheric forcing by the screen-level variables to the subgrid-scale. The atmospheric model is run with a 2.8 km horizontal grid resolution while the land surface processes are computed on a 400-m horizontal grid. The disaggregation of the driving atmospheric variables at screen-level is achieved by a three-step statistical downscaling with rules learnt from high-resolution fully coupled COSMO simulations, where both, atmosphere and surface, were simulated on a 400-m grid. The steps encompass spline interpolation of the grid scale variables, conditional regression based on the high-resolution runs, and an optional stochastic noise generator which restores the variability of the downscaled variables. Simulations for a number of case studies have been carried out, with or without mosaic surface representation and with or without atmospheric disaggregation, and evaluated with respect to the surface state variables and the turbulent surface exchange fluxes of sensible and latent heat. The results are compared with the high-resolution fully coupled COSMO simulations. The results clearly demonstrate the high importance of accounting for subgrid-scale surface heterogeneity. It is shown that the atmospheric disaggregation leads to clear additional improvements in the structures of the two-dimensional surface state variable fields, but to only marginally impacts on the simulation of the turbulent surface exchange fluxes. A detailed analysis of these results identifies strongly correlated errors in atmospheric and surface variables in the mosaic approach as the main reason for the latter. The effects of these errors largely cancel out in the flux parameterization, and thus explain the comparably good results for the fluxes in the mosaic approach without atmospheric disaggregation despite inferior performance for the surface state variables themselves. Inserting noise in the disaggregation scheme leads to a deterioration of the results.