We present a software tool for simulations of flow and multi-species solute transport in two and three-dimensional domains in combination with comprehensive intra-phase and inter-phase geochemistry. The software uses IPhreeqc, the coupling version of the geochemical modelling framework PHREEQC, as a reaction engine to the multi-purpose finite element solver COMSOL Multiphysics® for flow and transport simulations. Here we used COMSOL to solve Richards’ equation for aqueous phase flow in variably saturated porous media. The coupling procedure presented is in principle applicable to any simulation of aqueous phase flow and solute transport in COMSOL. The coupling with IPhreeqc gives major advantages over COMSOL’s built-in reaction capabilities, i.e., the soil solution is speciated from its element composition according to thermodynamic mass action equations with ion activity corrections. State-of-the-art adsorption models such as surface complexation with diffuse double layer calculations are accessible. In addition, IPhreeqc provides a framework to integrate user-defined kinetic reactions with possible dependencies on solution speciation (i.e., pH, saturation indices, and ion activities), allowing for modelling of microbially mediated reactions. Extensive compilations of geochemical reactions and their parameterization are accessible through associated databases.