Deep-sea mud shows tremendous resource potential for rare earth elements (REEs) and its formation is closely associated with sea-floor hydrothermal activities. Iron (oxyhydr)oxides link the sources and sinks of REEs in sea-floor hydrothermal systems. However, the complexation mechanisms of REEs on iron (oxyhydr)oxides have not been well understood yet. In this study, by using the first principles molecular dynamics technique, we first calculated the pKa’s of surface groups on goethite (110) surface at elevated temperatures relevant to sea-floor hydrothermal systems and then evaluated the complexation structures and free energies of REEs on goethite (110) and (010) surfaces using the method of constraint by taking Sc3+, Y3+, and La3+ as model REE cations. The results show that REE complexation occurs in mildly acidic to neutral conditions. The most thermodynamically stable complexes of REEs are bidentate complexes on two neighboring FeOH sites on goethite (110) surface and tridentate complexes on two neighboring FeOH sites plus one Fe2OH site on goethite (010) surface. Sc3+ complexes match the goethite lattice and can be incorporated into the lattice. The stabilities of REE complexes increase with the distance from hydrothermal vents. Complexation of Y3+ is less favored on goethite compared to other REEs whereas Sc3+ prefers complexation on goethite (010) surface and La3+ exhibits similar stabilities on both (010) and (110) surfaces. The derived atomic level complexation mechanisms would be helpful for the interpretation of experimental data and the prediction of REEs’ behavior in the sea-floor. The findings presented here provide valuable insights into REEs fractionation and enrichment in deep-sea muds.