Shale gas plays an increasingly important role in the current energy industry. Modeling of gas flow in shale media has become a crucial and useful tool to estimate shale gas production accurately. The second law of thermodynamics provides a theoretical criterion to justify any promising model, but it has been never fully considered in the existing models of shale gas. In this paper, a new mathematical model of gas flow in shale formations is proposed, which uses gas density instead of pressure as the primary variable. A distinctive feature of the model is to employ chemical potential gradient rather than pressure gradient as the primary driving force. This allows to prove that the proposed model obeys an energy dissipation law, and thus, the second law of thermodynamics is satisfied. Moreover, on the basis of energy factorization approach for the Helmholtz free energy density, an efficient, linear, energy stable semi-implicit numerical scheme is proposed for the proposed model. Numerical experiments are also performed to validate the model and numerical method.