A critical review of the existing theories of trace metal complexation in aquatic systems is presented. It was demonstrated that original titration curves, the corresponding log-log plot as well as different linearizing plots are not suitable for accurate determination of complexation parameters without additional correction procedures. Models of discrete and continuous distribution of binding strengths were compared. In both cases, the complexing capacity can be determined using the same method from the van den Berg-Ružić-Lee plot. The complexation of trace metals in aquatic systems is equivalent to the adsorption at heterogeneous binding sites (surface complexation). Exact numerical and analytical solutions for surface complexation at heterogeneous binding sites have been proposed for the analysis of titration data, based on the Langmuirian type of local isotherm and the Boltzmann distribution function of binding energies. This exact solution was compared with different semi-empirical isotherms, derived in literature and proposed for the analyses of titration data as well. These isotherms are mainly based on different quasi-Gaussian distribution functions. Disadvantages of these isotherms and their multicomponent versions are discussed. The generalized Langmuir and van Riemsdijk isotherms have a completely different shape (at large additions of trace metals during titration) from the shape of rigorous solution which agrees well with titration experiments of a wide titration window. In addition, the generalized Freundlich isotherm which has a shape similar to the shape of the rigorous solution is shifted from the rigorous solution in order to predict well the complete monolayer. It is concluded that the Boltzmann distribution function of binding energies should be used instead of the quasi-Gaussian type of distribution functions which are not suitable for the interpretation of titration experiments. The rigorous solution can predict the Langmuir, Freundlich, Dubinin-Radushkewich and Henry isotherms observed experimentally at different regions of surface coverages (trace metal concentrations).