A simulation-based interval stochastic bi-level multi-objective programming (SISBLMOP) model was proposed in this research, through integrating the global nutrient export from watersheds model, interval parameter programming and stochastic chance-constrained programming into a general bi-level multi-objective programming framework. The SISBLMOP model can handle multiple uncertainties expressed as discrete intervals and probability density functions in both the simulation and optimization processes. System complexities, including the hierarchy structure of upper- and lower-level decision makers, can also be addressed in the model. The proposed model is applied to a real-world case study of the Xinfengjiang Reservoir Watershed in South China to identify the satisfactory implementation levels of multiple best management practices (BMPs). The model results show that multiple BMP schemes for water quality management can be obtained under different upper- and lower-level decision-making and risk-violation scenarios, reflecting the cooperation and gaming results of the two-level decision makers. Consequently, the corresponding BMP implementation costs are acceptable to both the upper- and lower-level decision makers. The model is widely applicable and can be effectively used for water quality management under multiple uncertainties and complexities.