This paper describes a systematic study of the effect of fractal rough surfaces on the contact area with the aim of advancing contact constitutive laws used in the Discrete Element Method. An in-house Boundary Element code is adopted to investigate the mechanical behaviour of computer-generated surface roughness. Surfaces are generated to have methodically controlled root mean square height (Sq), root mean square gradient (Sdq), short wavevector (q0), large wavevector (q1), Hurst exponent (H) and fractal dimension (Df). The effect of each parameter on the contact area is investigated. Two recently proposed analytical solutions in tribology (i.e. Persson-Tosatti and Pastewka-Robbins) are applied to predict the real contact area. A parameter based on Sq, q0, and H is proposed and its sensitivity for real contact area prediction is demonstrated. Surfaces of natural sand are simulated and their mechanical response shows similar trend as the computer-generated surfaces, albeit more complex.