Adsorption of protons on a heterogeneous solid surface is modeled using the Monte Carlo (MC) simulation method. The surface of an oxide is assumed to consist of adsorption sites with pK assigned according to a quasi-Gaussian distribution. The influence of the electrostatic interactions combined with the energetic heterogeneity of the surface is examined, and the MC results are compared with the predictions of the mean field theory (MFT). It is demonstrated that the heterogeneity affects strongly the shape of the isotherms while it does not change the location of the common intersection point of the isotherms. On the other hand, introduction of repulsive interactions into the system is found to shift the CIP toward lower values of pH. It is also shown that the MFT, in general, describes correctly the behavior of the system. On the contrary the condensation approximation, used to derive relatively simple expressions for the adsorption isotherms, introduces serious errors unless the surface is strongly heterogeneous. Some practical remarks how to eliminate the errors associated both with the MC simulations and with the theory are also presented.