This dissertation presents the design of an integrated watershed model, WASH123D version 3.0, a first principle, physics-based watershed-scale model of integrated hydrology/hydraulics and water quality transport. This numerical model is comprised of three modules: (1) a one-dimensional (1-D) simulation module that is capable of simulating separated and coupled fluid flow, sediment transport and reaction-based water quality transport in river/stream/canal networks and through control structures; (2) a two-dimensional (2-D) simulation module, capable of simulating separated and coupled fluid flow, sediment transport, and reactive biogeochemical transport and transformation in two-dimensional overland flow systems; and (3) a three-dimensional (3-D) simulation module, capable of simulating separated and coupled fluid flow and reactive geochemical transport and transformation in three-dimensional variably saturated subsurface systems. The Saint Venant equation and its simplified versions, diffusion wave and kinematic wave forms, are employed for surface fluid flow simulations and the modified Richards equation is applied for subsurface flow simulation. The reaction-based advection-dispersion equation is used as the governing equation for water quality transport. Several physically and mathematically based numerical options are provided to solve these governing equations for different application purposes. The surface-subsurface water interactions are considered in the flow module and simulated on the basis of continuity of interface. In the transport simulations, fast/equilibrium reactions are decoupled from slow/kinetic reactions by the decomposition of reaction networks; this enables robust numerical integrations of the governing equation. Kinetic variables are adopted as primary dependent variables rather than biogeochemical species to reduce the number of transport equations and simplify the reaction terms. In each time step, hydrologic/hydraulic variables are solved in the flow module; kinetic variables are then solved in the transport module. This is followed by solving the reactive chemical system node by node to yield concentrations of all species. Application examples are presented to demonstrate the design capability of the model. This model may be of interest to environmental scientists, engineers and decision makers as a comprehensive assessment tool to reliably predict the fluid flow as well as sediment and contaminant transport on watershed scales so as to evaluate the efficacy and impact of alternative watershed management and remediation techniques prior to incurring expense in the field.