Reduced order modeling for the purpose of constructing a low dimensional model from high dimensional or infinite dimensional model has important applications in science and engineering such as fast model evaluations and optimization/control. A popular method for constructing reduced-order model is based on finding a suitable low dimensional basis by proper orthogonal decomposition (POD) and forming a model by Galerkin projection of the infinite dimensional model onto the basis. In this paper, we will discuss error estimates for Galerkin proper orthogonal decomposition method for an unsteady nonlinear coupled partial differential equations arising in viscous incompressible flows. A specific finite element in space and finite difference in time discretization scheme will be discussed.