Modeling transport of reactive solutes is a challenging problem, necessary for understanding the fate of pollutants and geochemical processes occurring in aquifers, rivers, estuaries, and oceans. Geochemical processes involving multiple reactive species are generally analyzed using advanced numerical codes. The resulting complexity has inhibited the development of analytical solutions for multicomponent heterogeneous reactions such as precipitation/dissolution. We present a procedure to solve groundwater reactive transport in the case of homogeneous and classical heterogeneous equilibrium reactions induced by mixing different waters. The methodology consists of four steps: (1) defining conservative components to decouple the solution of chemical equilibrium equations from species mass balances, (2) solving the transport equations for the conservative components, (3) performing speciation calculations to obtain concentrations of aqueous species, and (4) substituting the latter into the transport equations to evaluate reaction rates. We then obtain the space‐time distribution of concentrations and reaction rates. The key result is that when the equilibrium constant does not vary in space or time, the reaction rate is proportional to the rate of mixing, $\nabla ^{T}uD\nabla u$, where u is the vector of conservative components concentrations and $D$ is the dispersion tensor. The methodology can be used to test numerical codes by setting benchmark problems but also to derive closed‐form analytical solutions whenever steps 2 and 3 are simple, as illustrated by the application to a binary system. This application clearly elucidates that in a three‐dimensional problem both chemical and transport parameters are equally important in controlling the process.