The detailed solvation structure at the (110) surface of rutile (alpha-TiO2) in contact with bulk liquid water has been obtained primarily from experimentally verified classical molecular dynamics (CMD) simulations of the ab initio-optimized surface in contact with SPC/E water. The results are used to explicitly quantify H-bonding interactions, which are then used within the refined MUSIC model framework to predict surface oxygen protonation constants. Quantum mechanical molecular dynamics (QMD) simulations in the presence of freely dissociable water molecules produced H-bond distributions around deprotonated surface oxygens very similar to those obtained by CMD with nondissociable SPC/E water, thereby confirming that the less computationally intensive CMD simulations provide accurate H-bond information. Utilizing this H-bond information within the refined MUSIC model, along with manually adjusted Ti-O surface bond lengths that are nonetheless within 0.05 angstrom of those obtained from static density functional theory (DFT) calculations and measured in X-ray reflectivity experiments (as well as bulk crystal values), give surface protonation constants that result in a calculated zero net proton charge pH value (pH(znpc)) at 25 degrees C that agrees quantitatively with the experimentally determined value (5.4 +/- 0.2) for a specific rutile powder dominated by the (110) crystal face. Moreover, the predicted pH,,, values agree to within 0.1 pH unit with those measured at all temperatures between 10 and 250 degrees C. A slightly smaller manual adjustment of the DFT-derived Ti-O surface bond lengths was sufficient to bring the predicted pH,,,value of the rutile (110) surface at 25 degrees C into quantitative agreement with the experimental value (4.8 +/- 0.3) obtained from a polished and annealed rutile (110) single crystal surface in contact with dilute sodium nitrate solutions using second harmonic generation (SHG) intensity measurements as a function of ionic strength. Additionally, the H-bond interactions between protolyzable surface oxygen groups and water were found to be stronger than those between bulk water molecules at all temperatures investigated in our CMD simulations (25, 150 and 250 degrees C). Comparison with the protonation scheme previously determined for the (110) surface of isostructural cassiterite (alpha-SnO2) reveals that the greater extent of H-bonding on the latter surface, and in particular between water and the terminal hydroxyl group (Sn-OH) results in the predicted protonation constant for that group being lower than for the bridged oxygen (Sn-O-Sn), while the reverse is true for the rutile (110) surface. These results demonstrate the importance of H-bond structure in dictating surface protonation behavior, and that explicit use of this solvation structure within the refined MUSIC model framework results in predicted surface protonation constants that are also consistent with a variety of other experimental and computational data.