Multilevel Monte Carlo is an efficient variance reduction technique to quantify uncertainties in simulations of subsurface flows with highly varying, rough coefficients, e.g. in radioactive waste disposal. We study model elliptic problems of single phase flow in random media described by correlated lognormal distributions, discretized in space via (standard and mixed) finite elements. We will demonstrate (theoretically and numerically) significant gains with respect to standard Monte Carlo, leading (in the best case) to an asymptotic cost that is proportional to solving one deterministic PDE.