We develop a two-fluid model (TFM) for simulation of thermal transport coupled to particle migration in flows of non-Brownian suspensions. Specifically, we propose a closure relation for the inter-phase heat transfer coefficient of the TFM as a function of the particle volume fraction, particle diameter, magnitude of the particle phase's shear-rate tensor, and the thermal diffusivity of the particles. The effect of shear-induced migration in the particulate phase is captured through the use of state-of-the-art rheological closures. We validate the proposed interphase heat transfer coupling by calibrating it against previous experiments in a Couette cell. We find that, when the shear rate is controlled by the rotation of the inner cylinder, the shear and thermal gradients aid each other to increase the particle migration when temperature difference between the inner and outer walls, ΔT=T_in−T_out<0. Meanwhile, for ΔT>0, the shear and thermal gradients oppose each other, resulting in diminished particle migration, and a more uniform distribution of the particulate phase across the gap. Within the TFM framework, we identify the origin and functional form of a "thermo-rheological" migration force that rationalizes our observations. We also investigate the interplay of shear and thermal gradients in the presence of recirculating regions in an eccentric Couette cell (with offset axis and rotating inner cylinder). Simulations reveal that the system Nusselt number increases with the eccentricity E for ΔT>0, but a maximum occurs for ΔT<0 at E=0.4. This observation is explained by showing that, for E>0.4 and ΔT<0, significant flow recirculation enhances particle inhomogeneity, which in turn reduces heat transfer in the system (compared to ΔT>0).