Subsurface contamination problems of metals and radionuclides are ubiquitous. Metals and radionuclides may exist in the solute phase or may be bound to soil particles and interstitial portions of the geologic matrix. Accurate tools to reliably predict the migration and transformation of these metals and radionuclides in the subsurface environment enhance the ability of environmental scientists, engineers, and decision makers to analyze their impact and to evaluate the efficacy of alternative remediation techniques prior to incurring expense in the field. A mechanistic-based numerical model could provide such a tool. This paper communicates the development and verification of a mechanistically coupled fluid-flow thermal-reactive biogeochemical-transport model where both fast and slow reactions occur in porous and fractured media. Theoretical bases, numerical implementations, and numerical experiments using the model are described. A definition of the "rates" of fast/equilibrium reactions is presented to come up with a consistent set of governing equations. Two example problems are presented. The first one is a reactive transport problem which elucidates the non-isothermal effects on heterogeneous reactions. It also demonstrates that the rates of fast/equilibrium reactions are not necessarily greater than that of slow/kinetic reactions in the context of reactive transport. The second example focuses on a complicated but realistic advective-dispersive-reactive transport problem. This example exemplifies the need for innovative numerical algorithms to solve problems involving stiff geochemical reactions. It also demonstrates that rates of all fast/equilibrium reactions are finite and definite. Furthermore, it is noted that a species-versus-time curve cannot be used to characterize the rate of homogeneous fast/equilibrium reaction in a reactive transport system even if one and only one such reaction is responsible for the production of this species.