Complexation of rare earth elements (REEs) with humic substances (HSs) at various REE loading levels was studied by solvent extraction and equilibrium dialysis methods coupled with inductively coupled plasma-mass spectrometry (ICP-MS). Apparent stability constants (β_M-HS) of REE ions with standard Suwannee River humic acid (SRHA) were determined simultaneously for all REEs except for Pm. β_M-HS depends on the REE loading level to HSs, where β_M-HS decreases with an increase in the REE loading level at a lower REE loading region, whereas β_M-HS is almost constant beyond a certain REE loading level. The variation in β_M-HS with the REE loading level suggests the existence of several types of binding sites of REEs in HSs. Next, β_M-HS was converted into K_MA (which is the REE-HS interaction parameter for Tipping's “Humic Ion-Binding Model VI”) to be compared with the data by previous studies on REE-HS complexation at various REE loading levels in terms of relative differences among REEs (i.e., REE pattern). In the case of solvent extraction, the REE pattern of logK_MA changes with variation in the REE loading level; a peak around Sm and Eu was observed in the REE pattern of logK_MA at a higher REE loading region, whereas logK_MA increases with an increase in the atomic number of REEs at a lower REE loading region. On the other hand, such variation in the REE pattern of logK_MA was not observed in the results of equilibrium dialysis. This difference in the REE pattern of logK_MA can be explained by the total metal loading level, including the contributions of Fe³⁺ and Al³⁺ originally contained in HS. The main binding sites of REEs in HSs were determined by comparing the REE pattern of logK_MA experimentally obtained for HSs with those estimated based on the linear free-energy relationship (LFER) for carboxylates. Similarity of the REE pattern suggested that simple carboxylates are the main binding sites of REEs at a higher metal loading region, whereas chelating sites play a key role at a lower metal loading region. Our results show that relative stabilities of HS complexes among various REEs can change as a function of the metal loading level.