Surface complexation models (SCMs) based on Gouy-Chapman theory are often used to describe adsorption of ions onto mineral surfaces. To compensate for the buildup of charge at a solid surface, the composition of the electric diffuse layer next to the surface must balance the surface charge. To calculate the diffuse layer composition, several nonlinear equations and integrals must be solved, usually with an iterative approach. Before convergence, charge balance is typically not fulfilled. One numerical difficulty is that, because of these charge balance errors, the iterative solver may attempt to take the square root of a negative number. Herein, we show that for electrolytes containing only monovalent or divalent ions (i.e., most electrolytes encountered in practice), we can greatly simplify the integrals and eliminate the appearance of complex-valued integrands; it is even possible to derive explicit analytical formulas. Furthermore, using the new method prevents converging to non-physical roots of the Grahame equation, which links surface potential to surface charge. To the best of our knowledge, the presented formulation has not been implemented in geochemical modelling software before, although similar mathematical expressions have been presented in the literature. In Gouy-Chapman theory, ions can only be distinguished by their charges, but this is not consistent with all experimental findings. We present a model that allows for the preferential accumulation of ions in the diffuse layer. The model, which is implemented mathematically by including ion exchange sites with a variable exchange capacity, is flexible and more numerically tractable than the standard models for the diffuse layer.