The design of industrial wastewater treatment networks is a challenging problem due to the presence of non-convex bilinear terms in the mathematical formulation. Success in finding the global optimal solution is highly dependent on the quality of the relaxation that is used to compute the lower bound for the objective of minimizing the total flowrate going through the treatment units. We derive rigorous relaxations of increasing complexity for the two alternative nonlinear formulations that have been previously proposed. The main focus is on mixed-integer linear programming relaxations from multiparametric disaggregation and on finding out how to select the most convenient set of variables to discretize. Each such relaxation can be used as the basis of a two-stage algorithm for global optimization that iterates over an accuracy parameter that consecutively increases the tightness as well as the size of the underlying problem. Through the solution of sixteen benchmark problems from the literature we show that the computational performance of the algorithm for the best relaxation compares favorably with commercial solvers BARON and GloMIQO. (C) 2015 The Institution of Chemical Engineers. Published by Elsevier B.V. All rights reserved.