Groundwater management involves a complex decision-making process, often with the need to balance the trade-off between meeting society’s demand for water and environmental protection. Therefore effective management of groundwater resources often involves some form of multi-objective optimization (MOO). Many existing software tools offer simulation model-enabled optimization, including evolutionary algorithms, for solving MOO problems. However, such analyses involve a huge amount of numerical process-based model runs, which require significant computational effort, depending on the nonlinearity and dimensionality of the problem, in order to seek the optimal trade-off function known as the Pareto front. Surrogate modeling, through techniques such as Gaussian Process Regression (GPR), is an emerging approach to significantly reduce the number of these model evaluations thereby speeding up the optimization process. Yet, surrogate model predictive uncertainty remains a profound challenge for MOO, as the current Pareto dominance criteria presumes that model responses are deterministic. Such presumption could mislead surrogate-assisted optimization, which may result in either little computational savings from excessive retraining, or lead to suboptimal and/or infeasible solutions. In this work, we present probabilistic Pareto dominance criteria that considers the uncertainty of GPR emulation during MOO, producing a ”cloudy’” Pareto front which provides an efficient decision space sampling mechanism for retraining the GPR. We then developed a novel acquisition strategy to manage the solution repository from this cloud and generate an ensemble of infill points for retraining. We demonstrate the capabilities of the algorithm through benchmark test functions and a typical density-dependent coastal groundwater management problem.