Revealing the scattering behaviour of gas molecules on porous surfaces is essential to develop accurate boundary conditions for kinetic transport models that describe the gas dynamics in shale reservoirs. Here, we use high-fidelity molecular dynamics simulations to resolve the gas–surface interactions between methane molecules and realistic organic kerogen surfaces, and to assess the applicability of the widely used scattering kernels. Our results show that the tight matrix porosities have a negligible effect on the timescale and lengthscale of the scattering process, which can be considered instantaneous in time and local in space. Although reflected velocity distributions reveal that the common Maxwell, Cercignani–Lampis and Yamamoto scattering models fail to fully capture the scattering details of methane on kerogen, especially when the incident molecular speeds are high, the Maxwell model predicts best the reflected angular beam pattern and the overall reflected velocity distribution for rough kerogen surfaces. However, for low-speed impingement, more characteristic of shale applications, all scattering models give similar velocity distributions, which are driven by the high degree of gas–surface accommodation observed. We find that a Maxwell model with a calibrated tangential momentum accommodation coefficient, which approaches unity as the surface roughness increases to 2 , is enough to reproduce comparable velocity profiles and mass flow rates inside moderately confined kerogen mesopores. Deviations between the Maxwell model and our molecular simulations are only observed for highly rarefied transport problems, but this rarefaction lies beyond the realm of shale reservoir applications. This paper, therefore, reports the first scattering study on porous and rough kerogen surfaces, and demonstrates the applicability of the Maxwell model, which can be readily incorporated into gas kinetic solvers to predict the apparent permeability of shale with mesopore and macropore networks.