This chapter presents the development of the latest version of HYDROGEOCHEM a multidimensional numerical model of coupled fluid flow, thermal transport, hydrologic transport, and biogeochemical kinetic/equilibrium reactions in saturated/unsaturated media. It iteratively solves the Richards equation for fluid flow, the thermal transport equation for temperature fields, and reactive biogeochemical transport equations for concentration distributions. For the latter, the advectivedispersive-reactive transport equations are solved for mobile components and kinetic variables. The biogeochemical reaction equations along with the component-and kinetic-variable equations are solved for concentration distributions of all species. This version of HYDROGEOCHEM is designed for generic applications to reactive transport problems under non-isothermal conditions in subsurface media. It considers all types of reactions (aqueous complexation, adsorption-desorption, precipitationdissolution, ion-exchange, hydrolysis, and abiotic and biotic-mediated redox) as both fast/equilibrium and/or slow/kinetic processes, using a consistent definition of fast/equilibrium reaction rates. Input to the program includes the finite element numerical representation of the system, the properties of the media, reaction network, and initial and boundary conditions. Output includes the spatial distributions of pressure and total heads, velocity fields, moisture contents, temperature, and biogeochemical concentrations at user specified times and locations (finite element nodes). Six examples are employed to demonstrate the design capabilities of HYDROGEOCHEM, to illustrate the calculations of fast/equilibrium reaction rates, and to highlight the non-intuitive notion that the rates of slow/kinetic reactions are not necessarily smaller than those of fast/equilibrium reactions.