We quantify the effect of uncertainties on quantities of interest related to contact mechanics of rough surfaces. Specifically, we consider the problem of frictionless non adhesive normal contact between two semi infinite 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 field whose covariance function encodes the surface’s roughness, which is experimentally measurable. Then, we introduce the multilevel Monte Carlo method which is a computationally efficient sampling method for the computation of the expectation and higher statistical moments of uncertain system outputs, such as those derived from contact simulations. In particular, we consider two different quantities of interest, namely the contact area and the number of contact clusters, and show via numerical experiments that the multilevel Monte Carlo method offers significant computational gains compared to an approximation via a classic Monte Carlo sampling.