This paper investigates the interaction between stratification and turbulence by means of turbulence models. The standard and the advanced turbulent kinetic energy-dissipation (k-ϵ) model are derived theoretically, including algebraic stress relations. It is shown that a certain empirical constant in the standard model turns out to be a complicated implicit function in the advanced model, namely, a function of the turbulent shear number, the turbulent buoyancy number, and a wall correction. For a better understanding and physical interpretation of the k-ϵ models, an analysis is carried out for a simplified case where diffusive fluxes are neglected. For this idealization it is shown that (1) the flux Richardson number Rf has a certain lower bound R-f due to the establishment of convection, (2) a steady state flux Richardson number Rstf (which is defined here for this purpose) lebels the borderline between the tendency of turbulence to decrease or collapse (Rf>Rstf) or to increase (Rf<Rstf), and (3) the well-known upper limit for turbulent shear flow,R+f≈0.25, fits our theory. Using the standard model, the advanced model, a modified version of the level-2 model of Mellor and Yamada and a modified version of Kochergin's model, the evolution of thermal stratification in the northern North Sea during the Fladen Ground Experiment (FLEZ'76) is simulated numerically and compared with the measurements. In this specific application, the two k-ϵ models performed best. ¿ Amerian Geophysical Union 1995