This thesis presents the development of a numerical water quality model using a general paradigm of reaction-based approaches. In a reaction-based approach, all conceptualized biogeochemical processes are transformed into a reaction network. Through the decomposition of species governing equations via Gauss-Jordan column reduction of the reaction network, (1) redundant fast reactions and irrelevant kinetic reactions are removed from the system, which alleviates the problem of unnecessary and erroneous formulation and parameterization of these reactions, and (2) fast reactions and slow reactions are decoupled, which enables robust numerical integrations. The system of species transport equations is transformed to reaction-extent transport equations, which is then approximated with two subsets: algebraic equations and kinetic-variables transport equations. As a result, the model alleviates the needs of using simple partitions for fast reactions. With the diagonalization strategy, it makes the inclusion of arbitrary number of fast and kinetic reactions relatively easy, and, more importantly, it enables the formulation and parameterization of kinetic reactions one by one. To demonstrate the general paradigm, QAUL2E was recasted in the mode of a reaction network. The model then was applied to the Loxahatchee estuary to study its response to a hypothetical biogeochemical loading from its surrounding drainage. Preliminary results indicated that the model can simulate four interacting biogeochemical processes: algae kinetics, nitrogen cycle, phosphorus cycle, and dissolved oxygen balance.