In the United States, states are federally mandated to develop watershed management plans to mitigate pollution from increased impervious surfaces due to land development such as buildings, roadways, and parking lots. These plans require a major investment in water retention infrastructure, known as structural Best Management Practices (BMPs). However, the discovery of BMP configurations that simultaneously minimize implementation cost and pollutant load is a complex problem. While not required by law, an additional challenge is to find plans that not only meet current pollutant load targets, but also take into consideration anticipated changes in future precipitation patterns due to climate change. In this dissertation, a multi-scale, multiobjective optimization method is presented to tackle these three objectives. The method is demonstrated on the Bartlett Brook mixed-used impaired watershed in South Burlington, VT. New contributions of this work include: (A) A method for encouraging uniformity of spacing along the non-dominated front in multiobjective evolutionary optimization. This method is implemented in multiobjective differential evolution, is validated on standard benchmark biobjective problems, and is shown to outperform existing methods. (B) A procedure to use GIS data to estimate maximum feasible BMP locations and sizes in subwatersheds. (C) A multi-scale decomposition of the watershed management problem that precalculates the optimal cost BMP configuration across the entire range of possible treatment levels within each subwatershed. This one-time pre-computation greatly reduces computation during the evolutionary optimization and enables formulation of the problem as real-valued biobjective global optimization, thus permitting use of multiobjective differential evolution. (D) Discovery of a computationally efficient surrogate for sediment load. This surrogate is validated on nine real watersheds with different characteristics and is used in the initial stages of the evolutionary optimization to further reduce the computational burden. (E) A lexicographic approach for incorporating the third objective of finding non-dominated solutions that are also robust to climate change. (F) New visualization methods for discovering design principles from dominated solutions. These visualization methods are first demonstrated on simple truss and beam design problems and then used to provide insights into the design of complex watershed management plans. It is shown how applying these visualization methods to sensitivity data can help one discover solutions that are robust to uncertain forcing conditions. In particular, the visualization method is applied to discover new design principles that may make watershed management plans more robust to climate change.