A software tool for the simulation of one-dimensional unsaturated flow and solute transport together with biogeochemical reactions in the vadose zone was developed by integrating a numerical solution of Richards' equation into the geochemical modelling framework PHREEQC. Unlike existing software, the new simulation tool is fully based on existing capabilities of PHREEQC without source code modifications or coupling to external software packages. Because of the direct integration, the full set of PHREEQC's geochemical reactions with connected databases for reaction constants is immediately accessible. Thus, unsaturated flow and solute transport can be simulated together with complex solution speciation, equilibrium and kinetic mineral reactions, redox reactions, ion exchange reactions and surface adsorption including diffuse double layer calculations. Liquid phase flow is solved as a result of element advection, where the Darcy flux is calculated according to Richards' equation. For accurate representation of advection-dominated transport, a total- variation-diminishing scheme was implemented. Dispersion was simulated with a standard approach; however, PHREEQC's multicomponent transport capabilities can be used to simulate species-dependent diffusion. Since liquid phase saturation is recalculated after every reaction step, biogeochemical processes that modify liquid phase saturation, such as dissolution or precipitation of hydrated minerals are considered a priori. Implementation of the numerical schemes for flow and solute transport have been described, along with examples of the extensive code verification, before illustrating more advanced application examples. These include (i) the simulation of infiltration with saturation-dependent cation exchange capacity, (ii) changes in hydraulic properties due to mineral reactions and (iii) transport of mobile organic substances with complexation of heavy metals in varying geochemical conditions.