Iron oxyhydroxides are among the most important colloids that control rare earth elements (REE) concentrations and transport in natural hydrosystems. In this study, REE surface complexation to iron oxyhydroxides (Fe(OH)3(a)) was described by using the Donnan diffuse layer model and a two-site (i.e. ≡FesOH and ≡FewOH) model. The specific surface area and pH of zero charge were fixed as 100000 m2/mol and 8.0, respectively. The surface site density for weak and strong binding sites were fixed at 0.1 moL/mol Fe (≡FewOH) and 0.001 moL/mol Fe (≡FesOH) respectively. The two site types were used with pKa1int = 7.29 and pKa2int = 8.93. Using linear free energy relationship, the estimated equilibrium surface complexation constants (log K) increased from light REE (LREE) to heavy REE (HREE). Results of REE modeling calculation using the determined log K revealed a good fit of experimental data, showing an order of sorption on iron oxyhydroxides: HREE > MREE > LREE and preferential sorption of HREE at a lower pH. However, sorption edges only showed a slight change with ionic strength (0.1–0.7 moL/L) for the whole REE series. The generalized model was subsequently used to evaluate the impact of iron oxyhydroxides on REE speciation in groundwater. Application of the model to “model groundwater” showed that iron oxyhydroxide complexes of REE were significant in near neutral and weakly alkaline pH. This study contributes to putting forward a comprehensive database which would be useful for the application of surface complexation model to describe REE sorption by amorphous ferric hydroxides in nature.