In non-isothermal flows through porous media, the rate of heat transport across the domain and the distribution of temperature profile are strongly affected by the thermal properties of the fluid and solid phases, flow conditions, and the topological geometry of the pore structure. The combined effect of these parameters on the dynamics of heat transfer is generally characterized by the thermal dispersion coefficient. In this study, we model thermal dispersion in granular porous media using pore-scale numerical simulations. The flow and heat transport equation are solved in a porous domain at the pore-level to obtain the distribution of temperature inside the fluid phase and solid grains. The effluent temperature profile at the outlet is then fitted to the analytical solution of the macroscale heat advection-dispersion equation to determine the longitudinal component of the effective thermal dispersion tensor. Several simulation runs are carried out to investigate the effect of Peclet number, fluid to solid thermal conductivity and heat capacity ratios, and the temperature-dependent viscosity on the variation of the longitudinal thermal dispersion coefficient. The results indicate the existence of different thermal dispersion regimes separated by a threshold Peclet number. Below the threshold, thermal dispersion is dominated by conduction and longitudinal dispersion coefficient is equal to the effective thermal conductivity of the porous medium. Above the threshold, there is a transition region in which both conduction and advection are responsible for the heat ttansport. At higher Peclet numbers, dispersion regime is advection-dominated and normalized dispersion coefficient increases with the Peclet number and thermal conductivity ratio whereas it decreases with the heat capacity ratio. In this regime, the relationship between dispersion coefficient and Peclet number follows a power-law functionality with the exponent decreasing with the heat capacity ratio. To summarize the results, a characterization diagram is constructed for thermal dispersion regime based on the values of Peclet number and thermal conductivity ratio. Finally, the presence of viscous fingering due to the monotonically increasing viscosity along the flow is observed to further increase the longitudinal thermal dispersion. The outcomes of this study are of great importance in the design, implementation, and performance optimization of thermal processes in granular porous media.