Abstract Real surfaces in contact are subjected to damage such as wear and plastic deformation, which not only affects surface performances but may also lead to material failure. A proper prediction of surface damage is of critical importance in the design of advanced functional materials. However, the prediction is tremendously complicated both experimentally and theoretically by the presence of surface and subsurface imperfections formed during material manufacturing process. Recently, using Eshelby's equivalent inclusion method, we developed a general solution for multiple 3D arbitrarily-shaped inhomogeneous inclusions (which are imperfections characteristic of material dissimilarity and inelastic strain) near surfaces under contact loading. This solution takes into account interactions between all the inhomogeneous inclusions and between them and the loading body as well as considering surface roughness. It provides a detailed knowledge of surface deformation and pressure and subsurface elastic field. Furthermore, a layer of film was modeled as an inhomogeneous inclusion, leading to the successful modeling of elastic-plastic indentation on coated surfaces with imperfections. The inhomogeneous inclusion provides a unified framework to model various surface damage patterns (including chipping wear, gradual wear, and particle pull-out), competition between them, and surface evolution due to them.