Water splitting using an iron-based looping processes is a promising method to produce high purity hydrogen. The cyclical, heterogeneous oxidation and reduction reactions are carried out within high surface area stable porous structures as the solid reactant. The random walk method is used to simulate chemical reaction and species transport, as it is capable of handling stiff reaction kinetics and varying hydrodynamic dispersion tensors caused by pore-level velocity fluctuations. The original random walk formulation is extended to account for bulk density variations and source terms due to the chemical reaction. The species transport equation is recast in the form of the Fokker-Planck equation, and the trajectories of fluid particles are obtained by solving an appropriate Langevin equation that has additional drift terms as a result of spatial variations in bulk density and dispersion tensor. The source term is accounted for by changing the composition of fluid particles based on the reaction kinetics. The extended approach for each new term is validated against exact solutions or highly resolved finite difference solutions. Finally, a new coupled model for bulk fluid flow, species transport, and chemical reaction in porous media is developed and applied to simulate a bench-scale water splitting reactor with a porous iron-silica fixed bed structure. Excellent agreement with the measured hydrogen production rate at different operating conditions is obtained. Copyright (C) 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.