Numerical simulation of fluid-saturated granular materials subject to temporal or spatial variations in temperature is important for a range of industrial applications and to understand a number of natural processes. Amongst the available fluid-coupling options for Discrete Element Method (DEM), pore network models (PNMs) are attractive as the local heterogeneities in the flow field can be captured with a relatively low computational cost. Any coupled DEM-PNM framework that considers thermal behaviour must account for heat transfer between the solid particles and the fluid phase. This contribution develops a novel expression for the Nusselt number for use in coupled DEM-PNM formulations. Previous studies considering Nusselt number formulations for flow in pipes formed the basis for the proposed model. Finite volume simulations of laminar fluid flow through assemblies of spheres were used to assess the model's efficacy and to calibrate the empirical parameters.