In this paper, a numerical technique-the thermal particle dynamics method (TPD)-is extended to study heat conduction in granular media in the presence of stagnant interstitial fluids. The method, which generates a multi-particle simulation by explicitly modeling many two-particle interactions, allows bed heterogeneities to be directly included and dynamic temperature distributions to be obtained at the particle-level. Comparison with experimental data shows that TPD yields quantitatively accurate values of the effective thermal conductivity without introducing new adjustable parameters for a wide range of stagnant interstitial media. The model not only sheds light on fundamental issues in heat conduction in particulate materials, but also provides a valuable test bed for existing continuous theories.