We present, by means of a simple example, a comprehensive step-by-step procedure to consistently derive a pH model of aquatic systems. As pH modeling is inherently complex, we make every step of the model generation process explicit, thus ensuring conceptual, mathematical, and chemical correctness. Summed quantities, such as total inorganic carbon and total alkalinity, and the influences of modeled processes on them are consistently derived. The model is subsequently reformulated until numerically and computationally simple dynamical solutions, like a variation of the operator splitting approach (OSA) and the direct substitution approach (DSA), are obtained. As several solution methods are pointed out, connections between previous pH modelling approaches are established. The final reformulation of the system according to the DSA allows for quantification of the influences of kinetic processes on the rate of change of proton concentration in models containing multiple biogeochemical processes. These influences are calculated including the effect of re-equilibration of the system due to a set of acid-base reactions in local equilibrium. This possibility of quantifying influences of modeled processes on the pH makes the end-product of the described model generation procedure a powerful tool for understanding the internal pH dynamics of aquatic systems.