A multilevel Monte Carlo (MLMC) method for quantifying model-form uncertainties associated with the Reynolds-Averaged Navier-Stokes (RANS) simulations is presented. Two, high-dimensional, stochastic extensions of the RANS equations are considered to demonstrate the applicability of the MLMC method. The first approach is based on global perturbation of the baseline eddy viscosity field using a lognormal random field. A more general second extension is considered based on the work of [Xiao et al.(2017)], where the entire Reynolds Stress Tensor (RST) is perturbed while maintaining realizability. For two fundamental flows, we show that the MLMC method based on a hierarchy of meshes is asymptotically faster than plain Monte Carlo. Additionally, we demonstrate that for some flows an optimal multilevel estimator can be obtained for which the cost scales with the same order as a single CFD solve on the finest grid level.