We developed a Monte-Carlo-based radiative heat transfer model capable of simulating solar exposure and subsequent warming of rough snow and ice surfaces on ice-covered airless solar system bodies. The model accounts for wavelength-dependent internal light scattering and heat conduction in the snow interior down to meter-scale depths. We validated the model against analytical and experimental test cases with relevant applications to Europa, one of Jupiter's moons. We examined differential heating across the surface, from the centimeter to meter scale, to reveal potential patterns of preferential sublimation that could lead to rough ice morphologies, such as penitentes. An exploration of parameters such as penitente height-width ratios, shape, size, snow grain size, and thermal properties revealed that taller, thinner, larger penitentes with sharper peaks, coarser snow grain sizes, and lower thermal inertias are more likely to grow in Europa's environment near the equator.