The dependence of hydrogen coverage on the bulk doping concentration is investigated for the polar O-terminated (000) ZnO surface. We use hybrid density-functional theory in combination with ab initio thermodynamics to determine a doping-dependent phase diagram of this surface. For hydrogen coverages lower than 50% dangling oxygen bonds remain at the surface, where they subsequently become charged by bulk electrons. For such charged surfaces, a computational first-principles approach is presented, with which long-range band bending can now be included in first-principles supercell calculations. In this work, we use a modified and extended version of the recently introduced charge-reservoir electrostatic sheet technique (Sinai et al 2015 Phys. Rev. B 91 075311) to incorporate band bending effects directly into our first-principles calculations. This allows us to investigate the effect of space charge layers and the resulting band bending on the hydrogen coverage of the ZnO (000) surface. After introducing a new implementation of CREST, we show that the structure and stability of polar ZnO surfaces are indeed sensitive to the amount of free charge carriers in the bulk. At low doping concentrations our results corroborate the previously reported (2 × 1) hydrogen phase, at higher doping concentrations the hydrogen coverage diminishes notably.