A Gibbs energy minimization (GEM) approach to calculate chemical equilibria involving surface complexation at the mineral-water interfaces is presented, permitting the whole continuum "aqueous speciation - heterogeneous adsorption surface (co)precipitation - solid solution" to be modeled only in elemental stoichiometry, without additional material balance constraints for the total number of surface sites. The thermodynamic stability of surface-bound species is considered in a way similar to that of the solid-solution endmembers, gases, and aqueous species. This consideration is made possible by introducing: (1) the standard and reference states of surface species involving a single value of reference site density 170; (2) an expression of (electro) chemical potential linking quantities of the surface species with the amount and specific surface area of the sorbent; (3) activity/concentration relationships between aqueous sorbates and surface species, based on the surface activity terms (SAT) as functions of the maximum site densities Gamma(max); (4) the surface type area fractions phi(alpha,t) for describing faces or patches on heterogeneous surface of multi-site-surface sorpaon mineral phases; (5) an elemental stoichiometry and standard partial molal properties of the amphoteric neutral surface functional group =OH and derived surface complexes at ambient conditions, comparable between all (hydr)oxide mineral surfaces. For the =OH group (in 2pK(A) surface complexation models), a O(0.05)Hdegrees stoichiometry is conventionally defined from a reaction 0.5H(2)O(aq) = >O(0.5)Hdegrees, with logK(n) = 1.74436 (at Gamma(o) = 20 mumol . m(-2)) which is independent of temperature, pressure, pH, fO(2), amount and specific surface area of the sorbent, aqueous speciation, and surface speciation. It follows that all the non-reacted >O0.5HOo groups have constant activity in the presence of liquid H2O and are thus macroscopically indistinguishable. The surface heterogeneity is described via standard thermodynamic properties of the surface complexes (reacted groups), such as "surface hydronium" > O0.5H2+ "surface hydroxyl" > O-0.5(-), adsorbed cations or anions, as well as through and the common or individual Gamma(t,max) parameters. Theoretical findings, implemented in the GEM-Selektor code, are illustrated by modeling the literature potentiometric titration and adsorption data for silica and rutile in NaCl solutions with the triple-layer (TLM), the generalized double-layer (DLM). and the non-electrostatic (NEM) surface complexation models (SCM). Some examples were also calculated with the FITEQL3.2 code to check for compatibility. It is shown: (1) how approximate conversions between partial molal properties of surface species and intrinsic adsorption constants can be performed for different total site densities Gamma(C) not equal Gamma(o); (2) what are the regions of "geometrically ideal" and "non-ideal" adsorption behavior for the chosen SCM and how they are controlled by the SAT = f(Gamma(max)) on heterogeneous surfaces; and (3) why different sets of TLM parameters are not thermodynamically equivalent, and how a meaningful "optimal" set can be determined from titration data. Due to the separation of Gamma(max) parameters from Gamma(o) and, hence, from thermodynamic constants, the thermodynamic SCMs (in GEM implementation) can provide more insight into the relationships between surface types, complexes, bulk mineral sorbent, and aqueous sorbates, therefore permitting a more thorough use of microscopic/spectroscopic data in sorption modeling than has ever been possible.