A first principles study of the MgO(001) surface hydroxylation is reported. Using ab initio Hartree-Fock embedded cluster calculations, pair potentials for the species H+-O-surf, H+-Mg-surf, HO--O-surf, and HO--Mg-surf pairs have been obtained. These potentials are used to simulate, by means of molecular dynamics, the hydroxylation of the MgO(001) surface. The hydroxyl groups are found to remain at the surface, bound to Mg cations, inducing a noticeable relaxation. The protons are found to leave the surface, stabilizing beyond the third layer. This proton penetration agrees with recent elastic recoil detection analysis in which the protonation is observed to at least 5000 Angstrom. The structure of the surface, as well as the coordination environment of inner protons, is discussed in light of radial distribution functions and spectra of density of states.