In this work, we present a stable proper orthogonal decomposition-Galerkin approximation for parametrized steady incompressible Navier-Stokes equations with low Reynolds number. Supremizer solutions are added to the reduced velocity space in order to obtain a stable reduced-order system, considering in particular the fulfillment of an inf-sup condition. The stability analysis is first carried out from a theoretical standpoint and then confirmed by numerical tests performed on a parametrized two-dimensional backward-facing step flow. 1137 ological developments were provided in more recent years [15][16][17][18], including the application to new unconventional fields [19,20].In recent years, growing attention has been devoted to the combination of Galerkin strategies with POD [21-23] and to stabilization techniques both for POD [9,[24][25][26] and RB methods [27][28][29][30].The aim of this work is to exploit some features of the RB framework in its state-of-the-art formulation for the stable and accurate approximation of parametrized flows in a POD setting, in order to improve the performance of the latter. In particular, we aim at investigating possible sources of pressure instabilities, in order to avoid spurious pressure modes in the POD approximation of parametrized flows. In fact, it is very well known that the FE approximation of Navier-Stokes (NS) equations presents two main difficulties that can make their numerical approximations meaningless: (i) the indefinite nature of the system and (ii) the stability loss due to convection dominant regimes for high Reynolds numbers. The first issue can be cured by using suitable velocitypressure FE spaces (also known as inf-sup stable elements), satisfying a discrete version of the so-called Ladyzhenskaya-Babuška-Brezzi (LBB) condition (see Section 4). Instead, in order to overcome the velocity stability loss for convection dominant regimes, we can rely on FE stabilization techniques, such as the streamline upwind Petrov-Galerkin (SUPG) method. We underline that stabilization techniques can also be employed to overcome pressure instability, for instance, through the pressure stabilizing Petrov-Galerkin (PSPG) method. Concerning pressure instability, a possible way to cure it in the RB context is the enrichment of velocity spaces through the so-called supremizers.Here, we want to analyze this supremizer stabilization technique to exploit it within a POD context and to develop a stability analysis based on the introduction of an LBB inf-sup condition at the reduced level, too. In this way, POD could benefit from a previously developed robust framework for the RB stabilization of viscous flows that can also manage with a correct pressure recovery. In this work, we test such a method on nonlinear steady viscous flows modeled by NS equations and characterized by both physical and geometrical parametrization.In the offline stage of the resulting strategy, several NS truth solutions are computed, and a POD is performed to extract a low-dimensional representation of both velo...