This paper presents the development of a 3-D coupled model of subsurface flow, heat transfer, and reactive chemical transport where the local equilibrium assumption is applied. Strong coupling is used for steady-state simulations, while weak coupling may be employed for transient simulations to save computation time. Four types of boundary conditions and a complete suite of geochemical reactions are taken into account. The flow equation is discretized by the Galerkin finite element method. Both the reactive chemical transport and heat transfer equations are solved by either the hybrid Lagrangian-Eulerian or the conventional Eulerian finite element method. Two approaches are considered to compute the transport of aqueous components. The Picard method is used to handle non-linearity in both reactive chemical transport and unsaturated subsurface flow, while the Newton-Raphson method is employed to calculate chemical equilibrium. Two realistic, albeit hypothetical example problems are used to demonstrate the capability of the model. Copyright (C) 1998 Elsevier Science B.V.