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.,DNS, LES). Therefore, 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 sub-models are implemented for diffusion, pressure strain and scrambling, and dissipation terms. This work studies non-stratified and stratified towed wakes using SMC and DNS. Sub-models 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 conducted for the non-stratified and stratified wakes, respectively, to obtain converged higher-order statistics for direct comparison. SMC over-predicts wake height and under-predicts defect velocity, wake width, and turbulent kinetic and potential energies. 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.