The size and complexity of 3D unsteady simulations of aeronautical engines have been steadily increasing over the past decades and the exploitation of such simulations has become a challenge. In this paper, Galerkin projection based Reduced Order Models (ROM) s of unsteady turbulent incompressible flows, by the Proper Orthogonal Decomposition (POD) technique, are presented. The Snapshots POD 1, 2 applied to the flow velocity field guarantees the existence of a low dimensional space, which captures the large energetic scales of turbulence. Nevertheless, it is known that the small scales responsible for the dissipation of Turbulent Kinetic Energy (TKE) appear only within the neglected POD modes, because of the reduction process. So, the question which arises naturally is the following: Can we capture the dissipative scales, by a minimal dimensional POD generated space? In what follows, we give an affirmative answer for this question, thanks to an elaborated stable POD-Galerkin formulation of the unsteady turbulent incompressible and filtered Navier-Stokes equations, based on an improved representation of all the turbulent scales by the reduced order model. The key for this, is the enrichment of the flow velocity reduced POD subspace, by POD modes associated with the spatial gradient of the flow velocity. A new numerical approach is proposed to deal with the scaling issue, which results from combining snapshots of the velocity field and its gradient. In what follows, we detail the numerical approach for the stabilized ROM. The performance of this new approach is proved for an unsteady turbulent flow at Re= 80 000 in a complex geometry meshed with 30 million elements.