A series of high-order surface element discretization schemes for variational boundary element methods are introduced. The surface elements are chosen in accord with angular quadrature rules for integration of spherical harmonics. Surface element interactions are modeled by Coulomb integrals between spherical Gaussian functions with exponents chosen to reproduce the exact variational energy and Gauss’s law for a point charge in a spherical cavity. The present work allows high-order surface element expansions to be made for variational methods such as the conductorlike screening model for solvation and the variational electrostatic projection method for generalized solvent boundary potentials in molecular simulations.