In this paper, we present a new method to quantify the uncertainty introduced by the drastic dimensionality reduction commonly practiced in the field of computational fluid dynamics, the ultimate goal being to simulate accurate priors for real-time data assimilation. Our key ingredient is a stochastic Navier--Stokes closure mechanism that arises by assuming random unresolved flow components. This decomposition is carried out through Galerkin projection with a proper orthogonal decomposition (POD-Galerkin) basis. The residual velocity fields, model structure, and evolution of coefficients of the reduced order's solutions are used to compute the resulting multiplicative and additive noise's correlations. The low computational cost of these consistent correlation estimators makes them applicable to the study of complex fluid flows. This stochastic POD-reduced order model (POD-ROM) is applied to 2-dimensional and 3-dimensional direct numerical simulations of wake flows at Reynolds 100 and 300, respectively, with uncertainty quantification and forecasting outside the learning interval being the main focus. The proposed stochastic POD-ROM approach is shown to stabilize the unstable temporal coefficients and to maintain their variability under control, while exhibiting an impressively accurate predictive capability.