The reliability of water distribution networks is susceptible to large impacts resulting from uncertainties. Several factors may contribute to such uncertainties. These include nodal demands, reservoir/tank levels, and availability of system components. Additional sources of uncertainties arise when water-energy nexus is considered. For instance, when water desalination systems that are coupled with power plants are incorporated in the water distribution networks, uncertainties in energy demand influence the design of the water network system. In this work, a multi-objective optimization framework is developed for the design of water networks that involve dual-purpose power plants for simultaneous water desalination and power generation. The multi-objective model accounts for sustainability goals through environmental, economic, and societal metrics such as greenhouse gas emissions, net profit, and number of jobs created, respectively. Fossil fuels, as well as renewable biofuels and solar power, are considered as candidate sources of energy. Water resources include multiple systems such as dams, lakes, rivers, aquifers, and artificial storage lakes. Demands from agricultural, industrial as well as domestic needs in the region are considered. In order to incorporate uncertainty in the proposed approach, a stochastic optimization model is formulated with a probabilistic objective function and constraints. The network has multiple sources as well as demands with varied distributions. In order to find an overall distribution for probability density estimation, a decomposition scheme is applied where each scenario is solved separately as a sub-problem. The solution from sub-problem is used to find the probability using kernel density estimator. The modified model equations lead to a stochastic MINLP problem. The proposed framework and optimization formulations are applied to a case study in the Mexican state of Sonora. From the demand and natural resource distribution, 12 scenarios have been identified and the model is weighted by their probability of occurrence. The results are compared with the deterministic solution to obtain the cost associated with uncertainty.