Accurate models to reliably predict sediment and chemical transport in watershed water systems enhance the ability of environmental scientists, engineers and decision makers to analyze the impact of contamination problems and to evaluate the efficacy of alternative remediation techniques and management strategies prior to incurring expense in the field. This dissertation presents the conceptual and mathematical development of a general numerical model simulating (1) sediment and reactive chemical transport in river/stream networks of watershed systems; (2) sediment and reactive chemical transport in overland shallow water of watershed systems; and (3) reactive chemical transport in three-dimensional subsurface systems. Through the decomposition of the system of species transport equations via Gauss-Jordan column reduction of the reaction network, fast reactions and slow reactions are decoupled, which enables robust numerical integrations. Species reactive transport equations are transformed into two sets: nonlinear algebraic equations representing equilibrium reactions and transport equations of kinetic-variables in terms of kinetically controlled reaction rates. As a result, the model uses kinetic-variables instead of biogeochemical species as primary dependent variables, which reduces the number of transport equations and simplifies reaction terms in these equations. For each time step, we first solve the advective-dispersive transport of kinetic-variables. We then solve the reactive chemical system node by node to yield concentrations of all species. In order to obtain accurate, efficient and robust computations, five numerical options are provided to solve the advective-dispersive transport equations; and three coupling strategies are given to deal with the reactive chemistry. Verification examples are compared with analytical solutions to demonstrate the numerical accuracy of the code and to emphasize the need of implementing various numerical options and coupling strategies to deal with different types of problems for different application circumstances. Validation examples are presented to evaluate the ability of the model to replicate behavior observed in real systems. Hypothetical examples with complex reaction networks are employed to demonstrate the design capability of the model to handle field-scale problems involving both kinetic and equilibrium reactions. The deficiency of current practices in the water quality modeling is discussed and potential improvements over current practices using this model are addressed.