Natural shale is mainly characterized by various mineral compositions and multi-scale pore sizes, leading to complex adsorbate-adsorbent interactions, also named the energetic heterogeneity of shale rocks. In this study, a multi-site Langmuir model was proposed to address the energetic heterogeneity, where the interaction energy for each adsorption site was related to the multi-scale pore sizes and pore surface property. Besides, methane adsorptions on shale rocks were experimentally measured, and the proposed multi-site model was used to describe the measured adsorption isotherms. Compared to the traditional single-site Langmuir model, the multi-site model can not only accurately predict the adsorption isotherms but also describe the adsorption energy distributions for shale rocks. Furthermore, the proposed multi-site model was used to estimate the total Gas-In-Place (GIP) and investigate the gas contributed by individual pores. This estimation achieves more accurate results compared with the traditional method.