Abstract: Statistical prediction is often required in reservoir simulation to quantify production uncertainty or assess potential risks. Most existing uncertainty quantification procedures aim to decompose the input random field to independent random variables, and may suffer from the curse of dimensionality if the correlation scale is small compared to the domain size. In this work, we develop and test a new approach, K-means clustering assisted empirical modeling, for efficiently estimating waterflooding performance for multiple geological realizations. This method performs single-phase flow simulations in a large number of realizations, and uses K-means clustering to select only a few representatives, on which the two-phase flow simulations are implemented. The empirical models are then adopted to describe the relation between the single-phase solutions and the two-phase solutions using these representatives. Finally, the two-phase solutions in all realizations can be predicted using the empirical models readily. The method is applied to both 2D and 3D synthetic models and is shown to perform well in the P10, P50 and P90 of production rates, as well as the probability distributions as illustrated by cumulative density functions. It is able to capture the ensemble statistics of the Monte Carlo simulation results with a large number of realizations, and the computational cost is significantly reduced.