Static density functional theory (DFT) calculations with a continuous solvent model as well as classical and Car-Parrinello molecular dynamics (MD) simulations with explicit solvent molecules were performed to study the nature of Zr-monosubstituted monomeric and dimeric polyoxometalates (POMs) in water at different pHs. We have analyzed Zr-aqua, -hydroxo, and -aqua-hydroxo species derived from Linqvist- and Keggin-type anions. Both DFT and Car-Parrinello MD methods suggest that the Zr center tends to have coordination number greater than 6 and can bind up to 3 water molecules. Car-Parrinello MD simulations also show that the Zr atom fluctuates within the oxide POM framework, providing a flexible coordination environment. There is a small thermodynamic preference for the Zr-aqua species over the protonated Zr-hydroxo species; however the prevalence of one or the other species might depend on the pH. Classical MD simulations show that H3O+ interacts mainly with hydroxo ligand, while OH- anions prefer the protons of the H2O ligands. In general, an increase of the acidity favors the formation of Zr-aqua species, explaining why dimer dissociation is promoted at low pH. At basic conditions Zr-hydroxo species are generated, providing the reactive groups to form Zr center dot center dot center dot Zr linkages.