The continuum theory of elasticity and/or homogeneously discretized finite element models have been commonly used to investigate and analyze subsurface stresses in Hertzian contacts. These approaches, however, do not effectively capture the influence of the random microstructure topology on subsurface stress distributions in Hertzian contacts. In this paper, a finite element model for analyzing subsurface stresses in an elastic half-space subjected to a general Hertzian contact load with explicit consideration of the material microstructure topology is presented. The random internal geometry of polycrystalline microstructures is modeled using a 3D Voronoi tessellation, where each Voronoi cell represents a distinct material grain. The grains are then meshed using finite elements, and an algorithm was developed to eliminate poorly shaped elements resulting from “near degeneracy” in the Voronoi tessellations. Hertzian point and line contacts loads are applied as distributed surface loads, and the model’s response is evaluated with commercial finite element software ABAQUS . Internal stress results obtained from the current model compare well with analytical solutions from theory of elasticity. The influence of the internal microstructure topology on the subsurface stresses is demonstrated by analyzing the model’s response to an over rolling element using a critical plane approach.