Three-dimensional rough surfaces were generated using a modified two-variable Weierstrass-Mandelbrot function with fractal parameters determined from real surface images. The number and size of truncated asperities were assumed to follow power-law relations. A finite element model of a rigid sphere in normal contact with a semi-infinite elastic-plastic homogeneous medium was used to obtain a constitutive relation between the mean contact pressure, real contact area, and corresponding representative strain. The contact model was extended to layered media by modifying the constitutive equation of the homogeneous medium to include the effects of the mechanical properties of the layer and substrate materials and the layer thickness. Finite element simulations of an elastic-plastic layered medium indented by a rigid sphere validated the correctness of the modified contact model. Numerical results for the contact load and real contact area are presented for real surface topographies resembling those of magnetic recording heads and smooth rigid disks. The model yields insight into the evolution of elastic, elastic-plastic, and fully plastic deformation at the contact interface in terms of the maximum local surface interference. The dependence of the contact load and real contact area on the fractal parameters and the carbon overcoat thickness is interpreted in light of simulation results obtained for a tri-pad picoslider in contact with a smooth thin-film hard disk.