In this paper we propose a Bayesian method as a numerical way to correct and stabilise projection-based reduced order models (ROM) in computational fluid dynamics problems. The approach is of hybrid type, and consists of the classical proper orthogonal decomposition driven Galerkin projection of the laminar part of the governing equations, and Bayesian identification of the correction term mimicking both the turbulence model and possible ROM-related instabilities given the full order data. In this manner the classical ROM approach is translated to the parameter identification problem on a set of nonlinear ordinary differential equations. Computationally the inverse problem is solved with the help of the Gauss-Markov-Kalman smoother in both ensemble and square-root polynomial chaos expansion forms. To reduce the dimension of the posterior space, a novel global variance based sensitivity analysis is proposed.