The high-Reynolds-number turbulent flows of technological importance contain such a wide range of excited length and time scales that the application of direct or large-eddy simulations is all but impossible for the foreseeable future. Reynolds stress models remain the only viable means for the solution of these complex turbulent flows. It is widely believed that Reynolds stress models are completely ad hoc, having no formal connection with solutions of the full Navier-Stokes equations for turbulent flows. While this belief is largely warranted for the older eddy viscosity models of turbulence, it constitutes a far too pessimistic assessment of the current generation of Reynolds stress closures. It will be shown how secondorder closure models and two-equation models with an anisotropic eddy viscosity can be systematically derived from the Navier-Stokes equations when one overriding assumption is made: the turbulence is locally homogeneous and in equilibrium. A brief review of zero equation models and one equation models based on the Boussinesq eddy viscosity hypothesis will first be provided in order to gain a perspective on the earlier approaches to Reynolds stress modeling. It will, however, be argued that since turbulent flows contain length and time scales that change dramatically from one flow configuration to the next, two-equation models constitute the minimum level of closure that is physically acceptable. Typically, modeled transport equations are solved for the turbulent kinetic energy and dissipation rate from which the turbulent length and time scales are built up; this obviates the need to specify these scales in an ad hoc fashion. While two-equation models represent the minimum acceptable closure, second-order closure models constitute the most complex level of closure that is currently feasible from a computational standpoint. It will be shown how the former models follow from the latter in the equilibrium limit of homogeneous turbulence. However, the two-equation models that are formally consistent with second-order closures have an anisotropic eddy viscosity with strain-dependent coefficients - a feature that most of the commonly used models do not possess.