Upscaling is often needed in reservoir simulation to coarsen highly detailed geological descriptions. Most existing upscaling procedures aim to reproduce fine-scale results for a particular geological model (realization). In this work we develop and test a new approach, ensemble level upscaling, for efficiently generating upscaled two-phase flow parameters (e.g., upscaled relative permeabilities) for multiple geological realizations. For this purpose, flow-based upscaling calculations are combined with a statistical estimation procedure (cluster analysis). This approach allows us to numerically compute the upscaled two-phase flow functions for only a small portion of the coarse blocks. For the majority of blocks, these functions are estimated statistically based on single-phase velocity information (attributes), determined when the upscaled single-phase parameters are calculated. The procedure is designed to maintain close correspondence between the cumulative distribution functions for the numerically computed and statistically estimated two-phase flow functions. We apply the method to two-dimensional synthetic models of multiple realizations for uncertainty quantification. Models with different geological heterogeneity and fluid mobility ratios are considered. It is shown that the method consistently corrects the biases evident in primitive coarse-scale predictions and can capture the ensemble statistics (e.g., P50, P10, P90) of the fine-scale results almost as accurately as the full flow-based upscaling procedures, but with much less computational effort. The overall approach is flexible and can be used with any combination of upscaling procedures.