A continuous multiligand model was applied to the study of the systems Cu(II)-,Pb(II)-, Mn(II)-, and Cd(II)-humate. The model is essentially based on a bimodal Gaussian distribution of the binding sites. It allows the calculation of pH-independent metal-binding thermodynamic constants which are more representative than simple empirical curve-fitting parameters. This model was applied to the study of the potentiometric titration data of an humic acid extracted from an Andosol of Macomer (Sardinia, Italy) in order to calculate the formation constants of the complexes of humic acids with four metal ions; Its application to the potentiometric titration experiments showed that different classes of sites are involved in the binding of the metal ions considered. The affinity to the humic matrix followed the order: Cu(II) > Pb(II) > Mn(II) >> Cd(II); phenolic groups play a relevant role in Cu binding while Cd appears to bind exclusively to carboxylic groups.