Despite one fundamental issue in the adsorption theory of coalbed methane, little is known about the thermodynamic properties of gas adsorption in a porous coal matrix. In this work, considering the heterogeneity of pore structure and the exothermic characteristics of gas adsorption, a fractional heat conduction model with an unsteady volumetric heat source is proposed to study the heat transfer process induced by gas adsorption in a heterogeneous coal matrix. The heat conduction equation with a fractional time derivative is discretized by using an implicit numerical method based on the generalization of a standard finite-difference scheme. First, to validate the fractional heat conduction model, gas adsorption experiments on a microcalorimeter were carried out on 5 g coal samples of 0.3 mm diameter at 25 °C. The experimental heat flux with initial adsorption pressures of 3.23 bar, 5.83 bar and 9.77 bar increases rapidly from zero to peak values of 7.17 mW, 12.05 mW and 16.81 mW in less than 7 min (i.e., fast thermal diffusion stage) and then decreases slowly to zero again in approximately 2 h (i.e., slow thermal diffusion stage). It is revealed that for all tested gas pressures the fractional heat conduction model with a fractional order α=0.86 can reproduce the experimental process of heat flux with better accuracy than the Fourier law-based model (i.e., α=1), suggesting that anomalous thermal diffusion is the governing heat transfer process of gas adsorption in the coal matrix. Second, the spatial distribution and temporal evolution of temperature patterns with different model parameters are numerically simulated. It is found that the time to reach the peak temperature decreases from 760 s at the center of the coal particles to 490 s at the boundary. Finally, the parametric sensitivity of the thermodynamic properties of gas adsorption such as temperature, heat flux and integral adsorption heat is discussed in detail. Particularly, it is shown that as one of the most important thermodynamic parameters, the integral heat is very sensitive to the fractional order α. In the case of 3.23 bar, if α increases from0.75 to 1, while other model parameters remain unchanged, the integral heat could be enhanced from 1.1 J/g to 8.5 J/g.