Simple first-order closure models for covariances of concentration fluctuations, for use in modeling turbulent flow reactors, were tested by direct numerical simulations. Concentration covariances and other statistical functions were evaluated for a series-parallel reaction scheme in decaying, homogeneous turbulent flow. The simulations involve solving the unsteady Navier-Stokes and mass conservation equations by a pseudo-spectral method in a 64(3) wavenumber domain, with initially segregated reactants, for an initial turbulence Reynolds number of 29.9. Simulation results show that covariances of concentration fluctuations normalized with respect to mean concentration values are almost constant and that the time dependence of concentration covariances can be estimated if the mean concentrations are known at any one time after the initial time. Predictions of the first-order closure models of Bourne and Toor, Brodkey and Lewalle, Li and Toor, and Dutta and Tarbell were compared to simulation results. While none of these closures are satisfactory for all the conditions tested, the Brodkey-Lewalle closure agrees best with the simulations.