A New Formulation of the Surface Charge/Surface Potential Relationship in Electrolytes with Valence Less than Three

Summary Surface complexation reactions are frequently used to describe the adsorption of charged ions onto mineral surfaces. However, they can be challenging to implement robustly into numerical software, especially when coupling geochemical equilibrium calculations to transport. To compensate for the buildup of charge at a mineral surface, the composition of the electric diffuse layer next to the surface must balance the net charge of the surface. The mathematical relationship between surface charge and surface potential is obtained from the Poisson-Boltzmann equation. To calculate the composition of the diffusive layer explicitly, a time-consuming numerical integration must be performed inside each grid block at each Newton iteration. In addition to slowing down the simulations, the implementation in popular simulation software such as PHREEQC can sometimes fail to converge. In this paper we show that for electrolytes containing neutral, monovalent, divalent ions and/or complexes with valence less than three (i.e., most electrolytes encountered in practice), we can derive analytical expressions for the electrostatic terms in the governing equations. This allows us to determine the composition of the diffusive layer analytically, without evaluating the integral numerically. To our knowledge, this formulation is novel. For the cases tested so far, the new method greatly increases numerical stability and speed of convergence during Newton-Raphson iterations. Even in cases where diffuse layer concentrations are not required directly, incorrect solutions may be found. The relationship between surface charge and surface potential is then given by the Grahame equation, which has two solutions, only one of which is physical: the solution where surface charge and surface potential have the same sign. Converging to the correct solution during ordinary Newton-Raphson iterations is not guaranteed but depends strongly on the choice of initial guess for the solution to the geochemical system. To avoid this complication, we use our new formulation to derive a version of the Grahame equation having only a single solution, the physically correct one. Finally, we present a model that allows for the preferential accumulation of ions in the diffusive layer. The model, which is implemented mathematically by including ion exchange sites with variable exchange capacity, is consistent with recent data for ion transport in chalk.