In this paper an extension of the proper orthogonal decomposition to a multi-parameter domain is proposed. The method produces a low-order model after three consequent steps: construction of an optimal basis functions set with respect to model error; computing the model coefficients with the Galerkin approach; calibration of the coefficients to minimize the model error. As a particular example the 2D square cylinder wake flow is used to show the potential of the extended method for different laminar flow regimes. The proposed method shows good interpolation properties of the reduced model with respect to the defined range of Reynolds numbers. (C) 2010 Elsevier Ltd. All rights reserved.