We firstly employ a proper orthogonal decomposition (POD) method, Crank-Nicolson (CN) technique, and two local Gaussian integrals to establish a PODbased reduced-order stabilized CN mixed finite element (SCNMFE) formulation with very few degrees of freedom for non-stationary parabolized Navier-Stokes equations. Then, the error estimates of the reduced-order SCNMFE solutions, which are acted as a suggestion for choosing number of POD basis and a criterion for updating POD basis, and the algorithm implementation for the POD-based reduced-order SCN-MFE formulation are provided, respectively. Finally, some numerical experiments are presented to illustrate that the numerical results are consistent with theoretical conclusions. Moreover, it is shown that the reduced-order SCNMFE formulation is feasible and efficient for finding numerical solutions of the non-stationary parabolized Navier-Stokes equations.Keywords: proper orthogonal decomposition method, reduced-order stabilized Crank-Nicolson mixed finite element formulation, non-stationary parabolized Navier-Stokes equations, error estimate.A POD-Based Reduced-Order SCNMFE Formulation 347 Problem I. Seek u = (u 1 , u 2 ) T and p satisfying for T ∞ > 0,where Ω ⊂ R 2 is a bounded and connected domain, ∂ t = ∂/∂t represents the first order partial derivative with respect to time t, ∂ 2 yy = ∂ 2 /∂y 2 the second order partial derivative with respect to y, u = (u 1 , u 2 ) T the velocity vector, p the pressure, T ∞ the total time, ν = (RePr ) −1 , Re the Reynolds number, Pr the Prandtl number, and ϕ(x, y, t) and u 0 (x, y) are given vector functions. As a matter of convenience and without loss of generality, we might assume that ϕ(x, y, t) = 0 in the following theoretical analysis.A lot of numerical examples have shown that, if the fluid mainstream direction does not appear on a wide range of separation zone, the numerical results for the simplified parabolized Navier-Stokes equations are very close to those for the full Navier-Stokes equations. Especially, for fluid flow with high Reynolds number, the numerical viscosity for the full Navier-Stokes equations tends to hide some of the real physical viscosity. In other word, for fluid flow problems with high Reynolds number, the numerical solutions obtained from the parabolized Navier-Stokes equations are closer to real physical solutions than those obtained from the full Navier-Stokes equations (see [2]). Although, it seems, in principle, unreasonable that the fluid flow of the computational domain appearing separation zone is described by the simplified parabolized Navier-Stokes equations, a lot of computational examples show that, for high Reynolds number fluid flows, which feature a local small separation zone along the main stream direction (for example, the front or backward facing step flow, the separation bubble flow, the compression corner flow, the air intake channel flow field), the numerical solutions obtained by the simplified parabolized Navier-Stokes equations are very close to those obtained from the full Nav...