In this paper, we seek to optimize closed loop supply chain (CLSC) network design by determining the location of distribution and collection facilities and the assignment of forward and reverse flows that maximize the system profit and minimize CO2 emissions. To account for the difficulty of predicting and controlling recoverable items in practice, we incorporate into the proposed bi-objective model three types of uncertainties: return quantity, return quality and remanufacturing costs. We use the ε-constraint method to solve the resulting bi-objective scenario-based two-stage stochastic program and generate a set of Pareto-optimal solutions. Our numerical experiments show that both network configuration and performance levels are sensitive to variations of the quality and quantity of returns, in particular when the penalty for not processing returns is high. They also reveal the importance of stochastic modelling since the value of the stochastic solution can be significant in many cases. However, by increasing the size of the problem, the runtime of the ε-constrained stochastic model can be dramatically increased and memory shortage issues may occur. Therefore, we investigate the computational performance of a metaheuristic method based on a non-dominated sorting genetic algorithm (NSGAII) and a linear programming (LP) relaxation. Our method is capable of approximating the Pareto frontier on instances where the ε-constraint method cannot provide any feasible solution within the same time limit.