In this article, we present a model order reduction method based on time-scale analysis for chemical reaction systems. The method can be applied to any reaction system exhibiting multiple time scales and described by the set of differential equations dc/dt = f(c), where c (dimension n) is the vector of chemical species and f is the operator describing the kinetics. From the Jacobian of the kinetic operator and its eigenvalues, regions which exhibit different timescale behavior are identified. Within each region, the set of fast variables (dimension nf) is identified and these are linearly lumped into a smaller set of pseudo species. The fast and slow time scales can be separated, and the concentrations of the fast species can then be approximated by explicit algebraic expressions. Thus, the dynamics of the reaction system can be simulated by a smaller set of variables (dimension n-n(f)) characteristic of each region. The generated reduced models are non-stiff, computationally efficient, and valid over a range of initial conditions. The use of the method is illustrated by simplifying a description of cyclohexane oxidation.