The speciation and thermodynamic properties of Pb(II) complexes with different ligands in aqueous solutions are critical for understanding Pb mineralization in the crust and Pb pollution in soil and groundwater. Although HS − is an important ligand for Pb(II) in hydrothermal solutions, little is known about the structure and thermodynamic properties of Pb(II)−HS complexes, especially at temperatures higher than 350 °C. In this study, we applied ab initio molecular dynamics (AIMD) simulations to investigate the speciation and geometrical properties of different Pb(II)−HS complexes in sulfur-rich solutions from ambient to magmatic-hydrothermal conditions (25−600 °C, 1−2000 bar). The AIMD simulations show that Pb(II)−HS complexes are coordinated with 2 or 3 H 2 O molecules to form a 4-fold saddle-like structure at 25 °C and then transform into a 3-fold pyramidal structure at temperatures ≥200 °C. The complexes Pb(HS) 3 − , Pb(HS) 2 (H 2 O), and Pb(HS) 2 (OH) − are predominant in sulfur-rich solutions. The mixed−ligand complex Pb(HS) 2 (OH) − has not been described in previous studies, but we found that it is stable at high temperatures. In addition, the formation constants (log K) of different Pb(II)−HS complexes were calculated by thermodynamic integration simulations at 200−600 °C and 20−2000 bar. The calculated log K values are in good agreement with the theoretical and experimental data. Furthermore, the log K values were extrapolated to a wide range of temperatures and pressures based on the modified Ryzhenko−Bryzgalin model. The newly obtained thermodynamic data were used for thermodynamic modeling. The results demonstrate that Pb(HS) 2 (OH) − is dominant in high-temperature hydrothermal solutions with low sulfur concentration and high pH (>7). These results provide new insight into the complexation behavior of Pb(II) in sulfur-rich solutions.