The increasing number of physics parametrization schemes adopted in numerical weather forecasting models has resulted in a proliferation of intercomparison studies in recent years. Many of these studies concentrated on determining which parametrization yields results closest to observations rather than analyzing the reasons underlying the differences. In this work, we study the performance of two 1.5-order boundary layer parameterizations, the quasi-normal scale elimination (QNSE) and Mellor–Yamada–Janjić (MYJ) schemes, in the weather research and forecasting model. Our objectives are to isolate the effect of stability functions on the near-surface values and vertical profiles of virtual temperature, mixing ratio and wind speed. The results demonstrate that the QNSE stability functions yield better error statistics for 2 m virtual temperature but higher up the errors related to QNSE are slightly larger for virtual temperature and mixing ratio. A surprising finding is the sensitivity of the model results to the choice of the turbulent Prandtl number for neutral stratification (Prt0): in the Monin–Obukhov similarity function for heat, the choice of Prt0 is sometimes more important than the functional form of the similarity function itself. There is a stability-related dependence to this sensitivity: with increasing near-surface stability, the relative importance of the functional form increases. In near-neutral conditions, QNSE exhibits too strong vertical mixing attributed to the applied turbulent kinetic energy subroutine and the stability functions, including the effect of Prt0.