Biomass devolatilization and incineration in grate-fired plants are characterized by heterogeneous fuel mixtures, often incompletely mixed, dynamical processes in the fuel bed and on the particle scale, as well as heterogeneous and homogeneous chemistry. This makes modeling using detailed kinetics favorable but computationally expensive. Therefore, a computationally efficient model based on zero-dimensional stochastic reactors and reduced chemistry schemes, consisting of 83 gas-phase species and 18 species for surface reactions, is developed. Each reactor is enabled to account for the three phases: the solid phase, pore gas surrounding the solid, and the bulk gas. The stochastic reactors are connected to build a reactor network that represents the fuel bed in grate-fired furnaces. The use of stochastic reactors allows us to account for incompletely mixed fuel feeds, distributions of local temperature and local equivalence ratio within each reactor and the fuel bed. This allows us to predict the released gases and emission precursors more accurately than if a homogeneous reactor network approach was employed. The model approach is demonstrated by predicting pyrolysis conditions and two fuel beds of grate-fired plants from the literature. The developed approach can predict global operating parameters, such as the fuel bed length, species release to the freeboard, and species distributions within the fuel bed to a high degree of accuracy when compared to experiments. © 2020 American Chemical Society