The polynomial PM, integrated quadratic-linear IQLM and factorial analytic FAM solvation models are studied with the aim of estimating and optimizing the quaternary liquid-liquid equilibria. The reliability of existing models is analyzed with respect to three extraction factors by examining three literature quaternary systems, and the general ranking approximates the order PM > FAM > IQLM. A robust hybrid gradient method HGM is proposed for tackling the optimization problem on the basis of the first and second derivative conditions in combination with the interval gradient analysis of the modeled property. Given the nonsmooth nature of the optimization problem, we resort to a derivative solution scheme aimed at simplifying the convergent analysis of finding a global extremum with no a need for the evaluation of any objective function. The model sensitivity to mimic the global optimum of the extraction factors keeps to the order FAM >= IQLM > PM. This paper is completed by the numerical results of three test problems, manifesting the fact that there are significant benefits to be gained from the combination of gradient and interval analysis tools, particularly as the distribution of the solutes is a primary concern to be accounted for optimizing the quaternary equilibria. (C) 2020 Elsevier B.V. All rights reserved.