A novel high-fidelity CFD model is proposed for the process optimization and intensification of methane steam reforming in a packed bed reactor. The random packing of spherical particles is constituted by simulating a dumping of particles into the annular reactor by a commercial code based on discrete element method (DEM). The complex spatial geometry between hundreds of particles is created associated with a developed interparticle bridge method, which eliminates the interparticle and particle-wall contact regions by assuming that the small regions are stagnant and would not affect the reactor behavior. The conservation equations of mass, momentum, energy and chemical species are fully solved by employing a commercial CFD code based on finite volume method. A newly coded subroutine is incorporated into the CFD code to evaluate the reaction rates according to the species partial pressures and temperature on each surface cell covering the particles. The values of reaction, equilibrium and adsorption constants are specified based on the experiments in the literature. According to the reaction rates, the mass and heat sources are given to each fluid cell on the surface cell. It is verified that the developed CFD model is capable of predicting the distributions of species and temperature microscopically as well as macroscopically for methane steam reforming. The present DEM-CFD procedure enables one to investigate the behavior of packed bed reactor in detail once the catalyst is characterized in a conventional way. Therefore it would become a powerful tool to optimize and intensify any process in a packed bed reactor where the plug flow assumption is invalid.