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.The Reynolds-Averaged Navier-Stokes (RANS) equations combined with turbulence closure models are widely utilized in engineering to predict flows with high Reynolds number. These turbulence closure models are used to obtain an approximate Reynolds stress tensor that is responsible for coupling the mean flow with turbulence. Although many turbulence models exist in the literature, there is no single model that generalizes well to all classes of turbulent flows [1,2]. Specifically, the performance depends on the modeling assumptions and the type of flow used to calibrate the so-called closure coefficients that are needed as inputs to a turbulence model.Since the dominant source of error in the flow prediction comes from the turbulence modeling, a number of approaches have already been developed for the model-form uncertainty quantification (UQ) of RANS simulations, see e.g. [3,4] for recent reviews. The majority of these approaches are based on the perturbation of baseline RANS models. One way to achieve this is by injecting uncertainties in the closure coefficients [5,6,7,8] of turbulence models. Another more general physics-based approaches exists, which typically introduces randomness directly into the modeled Reynolds Stress Tensor (RST), either by perturbing its eigenvalues [9,10,11], tensor invariants [12,13] or the entire RST field [14]. One can also classify these stochastic models in terms of global and local perturbation (in space). For global approaches, such as in [5,6,7,10], the magnitude of perturbations in closure coefficients, eigenvalues of RST, etc. is the same throughout the flow domain. This translates to a low-dimensional UQ problem which can be efficiently handled by deterministic sampling techniques like stochastic collocation or just by simulating flows for Email addresses: pkumar@cwi.nl (Prashant Kumar), m.schmelzer@tudelft.nl (Martin Schmelzer), r.p.dwight@tudelft.nl (Richard P. Dwight) arXiv:1811.00872v1 [physics.comp-ph] 1 Nov 2018