In this paper, we study the contact stiffness of a fractal rough surface compressed by a rigid flat plane. A numerical model based on the analysis of flat punch indentation is proposed for simulated hierarchical surfaces, which are generated using statistical and fractal descriptors collected by surface profilometry. The contact stiffness of surfaces under increasing normal load is determined on the basis of the total truncated area at varying heights. The results are compared with experimental data from nanoindentation on four types of treated rough surfaces, showing good agreement with experimental observations below a certain truncation depth. Furthermore, the limits of the model's validity are discussed by focusing on surface geometries and deformation of contacting asperities. With this proposed truncation method, we present a parametric analysis to establish a correlation between contact stiffness and surface roughness descriptors. The contact stiffness shows a unified power-law scaling with respect to the applied load over a wide range for simulated surfaces with distinct sets of roughness descriptors. The exponent of the power-law relationship is found to correlate positively to the fractal dimension while its amplitude is inversely correlated to the surface roughness amplitude. This study provides an easily implemented and computationally efficient method to connect mechanical behaviour with multi-scale surface structure, which can be utilized in design and optimization of engineering applications involving rough contacts.