Density functional theory is used to compute the effect of protonation, deprotonation, and dehydroxylation of different reactive sites of a goethite surface modeled as a cluster containing six iron atoms constructed from a slab model of the (1 1 0) goethite surface. Solvent effects were treated at two different levels: (i) by inclusion of up to six water molecules explicitly into the quantum chemical calculation and (ii) by using additionally a continuum solvation model for the long-range interactions. Systematic studies were made in order to test the limit of the fully hydrated cluster surfaces by a monomolecular water layer. The main finding is that from the three different types of surface hydroxyl groups (hydroxo, μ-hydroxo, and μ3-hydroxo), the hydroxo group is most active for protonation whereas μ- and μ3-hydroxo sites undergo deprotonation more easily. Proton affinity constants (pKa values) were computed from appropriate protonation/deprotonation reactions for all sites investigated and compared to results obtained from the multisite complexation model (MUSIC). The approach used was validated for the consecutive deprotonation reactions of the [Fe(H2O)6]3+ complex in solution and good agreement between calculated and experimental pKa values was found. The computed pKa for all sites of the modeled goethite surface were used in the prediction of the pristine point of zero charge, pHPPZN. The obtained value of 9.1 fits well with published experimental values of 7.0-9.5.