Model development in hydrology and geochemistry has been advancing separately with limited integration. We developed a watershed hydrogeochemical code RT‐Flux‐PIHM to understand complex interactions between hydrological processes (PIHM), land‐surface processes (FLUX—Noah Land Surface Model), and multicomponent subsurface reactive transport (RT). The RT module simulates geochemical processes including aqueous complexation, surface complexation, mineral dissolution and precipitation, and cation exchange. The RT module is verified against the widely used reactive transport code CrunchFlow. The code uses semidiscrete finite volume method and irregular gridding and offers data harvesting capabilities from national databases. The application of RT‐Flux‐PIHM is demonstrated in the Susquehanna Shale Hills Critical Zone Observatory (SSHCZO). We aim to understand key processes that govern hydrogeochemical dynamics of the nonreactive chloride and reactive magnesium. Simulation results indicate that watershed characteristics, in particular topography, dictate the spatial distributions of water content and soil dissolution rates. Ion exchange provides buffering capacities and leads to a hysteresis loop of concentration and discharge relationship of magnesium, which differs from the open hysteresis of chloride. RT‐Flux‐PIHM offers physics‐based modeling capabilities to integrate the vast amount of water and chemistry data that have now become available, to differentiate the relative importance of competing processes, and to test hypotheses at the interface of hydrology and geochemistry.