Abstract. The definition of partial differential equation (PDE) models usually involves a set of parameters whose values may vary over a wide range. The solution of even a single set of parameter values may be quite expensive. In many cases, e.g., optimization, control, uncertainty quantification, and other settings, solutions are needed for many sets of parameter values. We consider the case of the time-dependent Navier-Stokes equations for which a recently developed ensemble-based method allows for the efficient determination of the multiple solutions corresponding to many parameter sets. The method uses the average of the multiple solutions at any time step to define a linear set of equations that determines the solutions at the next time step. To significantly further reduce the costs of determining multiple solutions of the Navier-Stokes equations, we incorporate a proper orthogonal decomposition (POD) reduced-order model into the ensemble-based method. The stability and convergence results for the ensemble-based method are extended to the ensemble-POD approach. Numerical experiments are provided that illustrate the accuracy and efficiency of computations determined using the new approach.Key words. Ensemble methods, proper orthogonal decomposition, reduced-order models, Navier-Stokes equations.1. Introduction. Computing an ensemble of solutions of fluid flow equations for a set of parameters or initial/boundary conditions for, e.g., quantifying uncertainty or sensitivity analyses or to make predictions, is a common procedure in many engineering and geophysical applications. One common problem faced in these calculations is the excessive cost in terms of both storage and computing time. Thanks to recent rapid advances in parallel computing as well as intensive research in ensemble-based data assimilation, it is now possible, in certain settings, to obtain reliable ensemble predictions using only a small set of realizations. Successful methods that are currently used to generate perturbations in initial conditions include the Bred-vector method, [30], the singular vector method, [5], and the ensemble transform Kalman filter, [4]. Despite all these efforts, the current level of available computing power is still insufficient to perform high-accuracy ensemble computations for applications that deal with large spatial scales such as numerical weather prediction. In such applications, spatial resolution is often sacrificed to reduce the total computational time. For these reasons the development of efficient methods that allow for fast calculation of flow ensembles at a sufficiently fine spatial resolution is of great practical interest and significance.Only recently, a first step was taken in [22,23] where a new algorithm was proposed for computing an ensemble of solutions of the time-dependent Navier-Stokes equations (NSE) with different initial condition and/or body forces. At each time step, the new method employs the same coefficient matrix for all ensemble members. This reduces the problem of solving multi...