We develop and test a multistate Gaussian model for the distribution of electrostatic solvation energies of a solute in liquid water. The multistate Gaussian model depends on the discovery of simple indicators of structural substates that individually display Gaussian fluctuations of electrostatic interactions. The probability distribution of electrostatic interactions is then modeled as a superposition of Gaussian distributions of electrostatic interactions of the substates. We find that the number of hydrogen bonds to the solute is a suitable substate indicator that eliminates the chief failures of single Gaussian models for the distribution of electrostatic interactions and of quadratic models of the electrostatic contribution to the excess chemical potential. These results should improve calculations of ionic chemical processes in water, i.e., acid-base chemistry, in particular those involving organic acids such as proteins and nucleic acids. The multistate Gaussian approach provides a specific and effective alternative to commonly discussed electrostriction and dielectric saturation modifications of dielectric continuum models. Moreover, the representation of complex energy distributions by a sum of simpler distributions based on structural substates is general and should be applicable in a variety of thermodynamic problems of solution chemistry.