Buoyant wakes encountered in the ocean environment are characterized by high Reynolds (Re) and Froude (Fr) numbers, leading to significant space–time resolution requirements for turbulence resolving CFD models (i.e., direct numerical simulations (DNS), large eddy simulations (LES)). Therefore, Reynolds-averaged Navier–Stokes (RANS) based models are attractive for these configurations. The inherently complex dynamics of stratified systems render eddy-viscosity-based modeling inappropriate. RANS second-moment closure (SMC) based modeling is more suitable because it accounts for flow anisotropy by solving the transport equations of important second-moment terms. Accordingly, eleven transport equations are solved at the SMC level, and a range of submodels are implemented for diffusion, pressure strain and scrambling, and dissipation terms. This work studies nonstratified and stratified towed wakes using SMC and DNS. Submodels in the SMC are evaluated in terms of how well their exact Reynolds averaged form impacts the accuracy of the full RANS closure. An ensemble average of 40 and 80–100 DNS realizations are required and conducted for these temporally evolving nonstratified and stratified wakes, respectively, to obtain converged higher-order statistics. SMC over-predicts wake height by over a factor of 2, and under-predicts defect velocity, wake width, and turbulent kinetic and potential energies by factors ranging from 1.3 to 3.5. Also, SMC predicts a near isotropic decay of normal Reynolds stresses , in contrast to the anisotropic decay returned by DNS. The DNS data also provide important insights related to the inaccuracy of the dissipation rate isotropy assumption and the non-negligible contribution of pressure diffusion terms. These results lead to several important recommendations for SMC modeling improvement.