ParCrunchFlow: an efficient, parallel reactive transport simulation tool for physically and chemically heterogeneous saturated subsurface environments. Understanding the interactions between physical, geochemical, and biological processes in the shallow subsurface is integral to the development of effective contamination remediation techniques, or the accurate quantification of nutrient fluxes and biogeochemical cycling. Hydrology is a primary control on the behavior of shallow subsurface environments and must be realistically represented if we hope to accurately model these systems. ParCrunchFlow is a new parallel reactive transport model that was created by coupling a multicomponent geochemical code (CrunchFlow) with a parallel hydrologic model (ParFlow). These models are coupled in an explicit operator-splitting manner. ParCrunchFlow can simulate three-dimensional multicomponent reactive transport in highly resolved, field-scale systems by taking advantage of ParFlow’s efficient parallelism and robust hydrologic abilities, and CrunchFlow’s extensive geochemical abilities. Here, the development of ParCrunchFlow is described and two simple verification simulations are presented. The parallel performance is evaluated and shows that ParCrunchFlow has the ability to simulate very large problems. A series of simulations involving the biologically mediated reduction of nitrate in a floodplain aquifer were conducted. These floodplain simulations show that this code enables us to represent more realistically the variability in chemical concentrations observed in many field-scale systems. The numerical formulation implemented in ParCrunchFlow minimizes numerical dispersion and allows the use of higher-order explicit advection schemes. The effects that numerical dispersion can have on finely resolved, field-scale reactive transport simulations have been evaluated. The smooth gradients produced by a first-order advection scheme create an artificial mixing effect, which decreases the spatial variance in solute concentrations and leads to an increase in overall reaction rates. The work presented here is the first step in a larger effort to couple these models in a transient, variably saturated surface-subsurface framework, with additional geochemical abilities.