Abstract. A general methodology is presented to integrate complex simulation models of hydrological systems into optimization models, as an alternative to scenario-based approaches. A gradient-based hill climbing algorithm is proposed to reach locally optimal solutions from distinct starting points. The gradient of the objective function is estimated numerically with the simulation model. A statistical procedure based on the Weibull distribution is used to build a confidence interval for the global optimum. The methodology is illustrated by an application to a small watershed in Ohio, where the decision variables are related to land-use allocations and the objective is to minimize peak runoff. The results suggest that this specific runoff function is convex in terms of the land-use variables, and that the global optimum has been reached. Modeling extensions and areas for further research are discussed.