The interplay of surface roughness and stable stratification is investigated by direct numerical simulation of Ekman flow. Our setup is well within the turbulent regime, reaching a friction Reynolds number of Reτ ≈ 2700. Further, we reach the verge of the fully rough regime under neutral conditions with a non-dimensional obstacle height H+ ≈ 40, corresponding to a z-nought parameter in viscous units z+ 0 ≈ 2. Stability is imposed via a gradual decrease of surface buoyancy from neutral (no stratification) to very strong stratification. The reduced Reynolds number (Reτ ) in comparison to atmospheric problems warrants consideration of viscous effects on our results, and we demonstrate a correction method that consistently incorporates viscous effects, thus reducing the spread of data from our numerical results. The weakly stable regime is maintained at higher stability due to efficient production of turbulence kinetic energy which counteracts buoyant restoring forces in the presence of roughness. When scaled according to Monin–Obukhov similarity theory (MOST) our results for weak stability compares excellent to known formulations based on atmospheric observations. The coefficients of the stability correction functions formomentum and heat are estimated as βm = 3.45, βh = 5.21 respectively, and we observe a slight but significant increase of the turbulent Prandtl number with stability. In the very stable regime, global flow properties (e.g. friction velocity, Obukhov length) oscillate with a decaying amplitude and global intermittency, i.e. the co-occurrence of turbulent/laminar fluid at large scale, is observed in the presence of roughness. In such very stable conditions, a strong veering of the surface wind with respect to the large-scale forcing (< 90◦) is observed.