First-principles molecular dynamics simulations are carried out to study the stepwise hydrolysis reactions of Y3+-aqua in water. The simulation results reproduce the 8-fold coordinated hydration structure of Y3+. A constrained molecular dynamics method is employed to investigate the hydrolysis reactions, and a thermodynamic integration technique is used to calculate the free-energy profiles. Thus, the stepwise pKa's are derived: pKa1=7.9, pKa2=9.1 and pKa3=10.3, while the fourth hydrolysis is found to be extremely unfavorable. Based on the dissociation constants, the hydrolyzing speciation of Y3+ is constructed. The simulation results also show that in the hydrolysis process, some H2O ligands leave the cation to become part of the bulk solution, which makes the hydration number of Y3+ cations decrease from 8 (Y(H2O)83+) to 7 (Y(OH)(H2O)62+ and Y(OH)2(H2O)5+) and then to 6 (Y(OH)3(H2O)3).