This paper describes the development of an optimal biochemical oxygen demand (BOD) load model for the Beargrass Creek watershed in Louisville, Kentucky. A conceptual, macro-level BOD load model is first developed based on the steady state Streeter-Phelps equations for dissolved oxygen deficit for the three forks of the watershed. The model is used to compute the effective total ultimate BOD concentration for the stream system based on observed dissolved oxygen deficits and stream flows measured in the field. A mass balance of flows and corresponding BOD concentrations is performed to compute the corresponding loads at the junction points of the streams. The computed loads are disaggregated into three distinct types namely, point source, non-point source, and an unknown source attributed to other suspect sources of pollution. The BOD loading model for the watershed is used as a simulation model in a disaggregated nonlinear constrained optimization framework to obtain optimal load reductions for all three types of pollution sources for the stream system. The optimization framework uses two different types of optimization algorithms namely genetic algorithms and the box complex method of constrained optimization. The optimization framework can serve as a useful management tool for watershed decision makers to formulate and evaluate different management strategies leading to effective capital improvement projects for the watershed. Such a framework can be effectively used to determine total maximum daily load (TMDL) for streams that are impaired for low dissolved oxygen and/or nutrient enrichment.