The determination of the pollutant load distribution among different pollutant sources in a watershed is a critical step in total maximum daily load (TMDL) development. Under current TMDL practices, TMDL allocations are typically determined through a trial-and-error approach of reducing pollutant loadings until a watershed simulation model predicts that water quality standards will be met given a margin of safety. Unfortunately, many feasible combinations of load reductions and significant uncertainties may exist. Therefore it is difficult and time-consuming to compare various allocation scenarios using a trial-and-error approach. A robust optimization model is developed in this study to incorporate the uncertainty of water quality predictions and to minimize pollutant load reductions given various levels of reliability with respect to the water quality standards. The generalized likelihood uncertainty estimation approach is used to explicitly address the uncertainty of a watershed simulation model, the Hydrological Simulation Program--Fortran. The uncertainty is integrated into TMDL allocations using a robust genetic algorithm model linked with a response matrix approach. The developed robust optimization model is demonstrated using a case study based on the Moore's Creek fecal coliform TMDL study. The trade-offs between reliability levels and total load reductions of allocation scenarios are evaluated, and the optimized load reduction scenarios are compared with the scenario generated by a trial-and-error approach and approved by the U.S. Environmental Protection Agency. The results show that the optimized load reduction scenario requires 30% less load reductions than the scenario approved by the U.S. Environmental Protection Agency at the same reliability level.