The stiffness in the top surface of many biological entities like cornea or articular cartilage, as well as chemically cross-linked synthetic hydrogels, can be significantly lower or more compliant than the bulk. When such a heterogeneous surface comes into contact, the contacting load is distributed differently from typical contact models. The mechanical response under indentation loading of a surface with a gradient of stiffness is a complex, integrated response that necessarily includes the heterogeneity. In this work, we identify empirical contact models between a rigid indenter and gradient elastic surfaces by numerically simulating quasi-static indentation. Three key case studies revealed the specific ways in which (I) continuous gradients, (II) laminate-layer gradients, and (III) alternating gradients generate new contact mechanics at the shallow-depth limit. Validation of the simulation-generated models was done by micro- and nanoindentation experiments on polyacrylamide samples synthesized to have a softer gradient surface layer. The field of stress and stretch in the subsurface as visualized from the simulations also reveals that the gradient layers become confined, which pushes the stretch fields closer to the surface and radially outward. Thus, contact areas are larger than expected, and average contact pressures are lower than predicted by the Hertz model. The overall findings of this work are new contact models and the mechanisms by which they change. These models allow a more accurate interpretation of the plethora of indentation data on surface gradient soft matter (biological and synthetic) as well as a better prediction of the force response to gradient soft surfaces. This work provides examples of how gradient hydrogel surfaces control the subsurface stress distribution and loading response.