Heterogeneity and a paucity of measurements of key material properties undermine the veracity of quantitative predictions of subsurface flow and transport. For such model forecasts to be useful as a management tool, they must be accompanied by computationally expensive uncertainty quantification, which yields confidence intervals, probability of exceedance, etc. We design and implement novel Multilevel Monte Carlo (MLMC) algorithms that accelerate estimation of the cumulative distribution functions (CDFs) of quantities of interest, e.g., water breakthrough time or oil production rate. Compared to standard non-smoothed MLMC, the new estimators achieve a significant variance reduction at each discretization level by smoothing the indicator function with a Gaussian kernel or replacing standard Monte Carlo (MC) with the recently developed Hierarchical Latinized Stratified Sampling (HLSS). After validating the kernel-smoothed MLMC and HLSS-enhanced MLMC methods on a single-phase flow test bed, we demonstrate that they are orders of magnitude faster than standard MC for estimating the CDF of breakthrough times in multi-phase flow problems.