We quantify the eect of uncertainties on quantities of interest related to contact mechanics of rough surfaces. Specically, we consider the problem of frictionless non adhesive normal contact between two semi innite linear elastic solids subject to uncertainties. These uncertainties may for example originate from an incomplete surface description. To account for surface uncertainties, we model a rough surface as a suitable Gaussian random eld whose covariance function encodes the surface's roughness, which is experimentally measurable. Then, we introduce the multilevel Monte Carlo method which is a computationally ecient sampling method for the computation of the expectation and higher statistical moments of uncertain system output's, such as those derived from contact simulations. In particular, we consider two dierent quantities of interest, namely the contact area and the number of contact clusters, and show via numerical experiments that the multilevel Monte Carlo method oers signicant computational gains compared to an approximation via a classic Monte Carlo sampling.