Buoyant shear layers encountered in many engineering and environmental applications have been studied by researchers for decades. Often, these flows have high Reynolds and Richardson numbers, which leads to significant/intractable space–time resolution requirements for direct numerical simulation (DNS) or large eddy simulation (LES). On the other hand, many of the important physical mechanisms, such as stress anisotropy, wake stabilization, and regime transition, inherently render eddy viscosity-based Reynolds-averaged Navier–Stokes (RANS) modeling inappropriate. Accordingly, we pursue second-moment closure (SMC), i.e., full Reynolds stress/flux/variance modeling, for moderate Reynolds number nonstratified, and stratified shear layers for which DNS is possible. A range of submodel complexity is pursued for the diffusion of stresses, density fluxes and variance, pressure strain and scrambling, and dissipation. These submodels are evaluated in terms of how well they are represented by DNS in comparison to the exact Reynolds-averaged terms, and how well they impact the accuracy of full RANS closure. For the nonstratified case, SMC model predicts the shear layer growth rate and Reynolds shear stress profiles accurately. Stress anisotropy and budgets are captured only qualitatively. Comparing DNS of exact and modeled terms, inconsistencies in model performance and assumptions are observed, including inaccurate prediction of individual statistics, non-negligible pressure diffusion, and dissipation anisotropy. For the stratified case, shear layer and gradient Richardson number growth rates, and stress, flux and variance decay rates, are captured with less accuracy than corresponding flow parameters in the nonstratified case. These studies lead to several recommendations for model improvement.