A heterogeneous complexation model approach has been developed to describe the adsorption of large organic molecules by goethite taking the full speciation of the adsorbed molecules into account. The essence of the model is the calculation of the mean mode of an adsorbed organic molecule, defined by an overall affinity, charge distribution, and reaction stoichiometry. The calculation of the mean mode of the adsorbed organic molecule is based on the pH-dependent interaction of a functional group for binding a proton or forming an inner or outer sphere complex with specific surface sites. The distribution of the organic molecules over the solid and solution phases is calculated using the CD-MUSIC model with parameters obtained from the mean mode of the adsorbed molecule. The mean mode of the adsorbed molecule is calculated using the competitive Langmuir or NICCA equation in combination with an electrostatic model comprising charge distribution, taking the competitive binding of functional groups with protons or surface sites into account. The new model is applied to the adsorption data of a series of benzenecarboxylic acids. The model approach can describe simultaneously the concentration, pH, and salt dependency of the benzenecarboxylate adsorption. Furthermore, the model is applied to describe data obtained from IR spectroscopy. The new approach predicts the data reasonably well using only three adjustable parameters in case of the NICCA-CD-MUSIC approach. The average stoichiometry of the overall reaction is not an input parameter of the model but is an output parameter. In future work the new model will be applied to the adsorption of larger organic acids, e.g., fulvic acid.