A numerical model was developed for the prediction of the density stratification of lakes and reservoirs. It combines a buoyancy-extended k-ε model with a seiche excitation and damping model to predict the diffusivity below the surface mixed layer. The model was applied to predict the seasonal development of temperature stratification and turbulent diffusivity in two medium-sized lakes over time periods ranging from 3 weeks to 2 years. Depending on the type of boundary condition for temperature, two or three model parameters were optimized to calibrate the model. The agreement between the simulated and the observed temperature distributions is excellent, in particular, if lake surface temperatures were prescribed as surface boundary condition instead of temperature gradients derived from heat fluxes. Comparison of different model variants revealed that inclusion of horizontal pressure gradients and/or stability functions is not required to provide good agreement between model results and data. With the aid of uncertainty analysis it is shown that the depth of the mixed surface layer during the stratified period could be predicted accurately within ¿1 m. The sensitivity of the model to several parameters is discussed.