In the case of expensive black-box simulators, efficient surrogate-based methods to perform multi-objective optimization have been proposed. Many rely on Gaussian Process (GP) models to identify optimal compromise solutions in the design and objective spaces, i.e., the Pareto set and Pareto front, respectively. Yet, when it comes to predicting the location of these sets, considering solely the predictive mean of the GP models is insufficient, as it does not propagate uncertainties. To this end, we propose the use of GP conditional simulations of Pareto sets and Pareto fronts. From the latter, based on concepts from random closed set theory, we obtain a continuous approximation of the Pareto front and a measure of the associated variability. We also discuss extensions to enable interactivity, to exploit the Pareto set low-dimensional structure or to cope with high-dimensional inputs. Finally, we show how this methodology naturally adapts to accommodate noise in observations.