Abstract. We consider settings for which one needs to perform multiple flow simulations based on the Navier-Stokes equations, each having different values for the physical parameters and/or different initial condition data, boundary conditions data, and/or forcing functions. For such settings, we propose a second-order time accurate ensemble-based method that to simulate the whole set of solutions, requires, at each time step, the solution of only a single linear system with multiple righthand-side vectors. Rigorous analyses are given proving the conditional stability and error estimates for the proposed algorithm. Numerical experiments are provided that illustrate the analyses.Key words. Navier-Stokes equations, parameterized flow, ensemble method MSC (2010). 65M60, 76D051. Introduction. Many computational fluid dynamics applications require multiple simulations of a flow under different input conditions. For example, the ensemble Kalman filter approach used in data assimilation first simulates a forward model a large number of times by perturbing either the initial condition data, boundary condition data, or uncertain parameters, then corrects the model based on the model forecasts and observational data. A second example is the construction of low-dimensional surrogates for solutions of partial differential equation (PDE) such as sparse-grid interpolants or proper orthogonal decomposition approximations for which one has to first obtain expensive approximations of solutions corresponding to several parameter samples. Another example is sensitivity analyses of solutions for which one often has to determine approximate solutions for a number of perturbed inputs such as the values of certain physical parameters. In this paper, we consider such applications and develop a second-order time-stepping scheme for efficiently simulating an ensemble of flows. In particular, we consider the setting in which one wishes to determine the PDE solutions for several different values of the physical parameters and with several different choices of initial condition and boundary condition data and forcing functions appearing in the PDE model.The ensemble algorithm we use was first developed in [11] to find a set of J solutions of the Navier-Stokes equations (NSE) subject to different initial condition and forcing functions. The main idea is that, based on the introduction of an ensemble average and a special semi-implicit time discretization, the discrete systems for the multiple flow simulations share a common coefficient matrix. Thus, instead of solving J linear system with J right-hand sides (RHS), one only need solve one linear system with J RHS. This leads to great computational saving in linear solvers when either the LU factorization (for small-scale systems) or a block iterative algorithm (for large-