A multi-site adsorption model was employed to determine the supercritical methane adsorption in shales, in which the heterogeneity of surface energy distribution was considered. Additionally, the high-pressure (similar to 20 MPa) methane adsorption behaviors at different temperature conditions (303 K, 333 K and 363 K) were experimentally investigated, and both the traditional single-site model and the multi-site model were employed to analyze the experimental data. Results show that (i) although these two models can achieve acceptable fitting results with the measured adsorption amounts, the adsorption heat interpreted by them shows significant discrepancy. The adsorption heat calculated by the single-site model is always a constant, while it obtained by the multi-site model decreases with an increase in surface coverage. This coverage-dependent heat directly indicates the nature of the surface-energy heterogeneity of our studied shale samples. (ii) The adsorption heat at the zero-coverage limit condition calculated by Henry's law, indicating the actual adsorbate-adsorbent interaction, is significantly larger than that calculated by the single-site model; but the zero-coverage heat can be well matched with the value determined by the multi-site model. Therefore, the present multi-site model can accurately depict not only the adsorption isotherms but also the thermodynamic parameters in energetically heterogeneous shales.