Large-eddy simulation (LES) models are widely used to study atmospheric turbulence. The effects of small-scale motions that cannot be resolved need to be modeled by a subfilter-scale (SFS) model. The SFS contribution to the turbulent fluxes is typically significant in the surface layer. This study presents analytical solutions of the classical Smagorinsky SFS turbulent kinetic energy (TKE) model including a buoyancy flux contribution. Both a constant length scale and a stability-dependent one as proposed by Deardorff are considered. Analytical expressions for the mixing functions are derived and Monin-Obukhov similarity relations that are implicitly imposed by the SFS TKE model are diagnosed. For neutral and weakly stable conditions, observations indicate that the turbulent Prandtl number (PrT) is close to unity. However, based on observations in the convective boundary layer, a lower value for PrT is often applied in LES models. As a lower Prandtl number promotes a stronger mixing of heat, this may cause excessive mixing, which is quantified from a direct comparison of the mixing function as imposed by the SFS TKE model with empirical fits from field observations. For a strong stability, the diagnosed mixing functions for both momentum and heat are larger than observed. The problem of excessive mixing will be enhanced for anisotropic grids. The findings are also relevant for high-resolution numerical weather prediction models that use a Smagorinsky-type TKE closure.