A two-dimensional flow and hydrogeochemical transport model is designed for estimating the movement and attenuation of chemicals after they are released into subsurface systems. This model combines the flow of groundwater and the transport of multispecies-multicomponent chemical systems. The flow of groundwater is calculated basically according to the updated FEMWATER model, which solves the Richards' equation by the finite element method. The transport of chemicals is determined by modifying the LEHGC code which solves the coupled geochemical equilibrium and hydrological transport equations using the Lagrangian-Eulerian approach. The modification is made by casting the multicomponent transport equations in a modified implicit iteration form, and the retarded velocities are used in the Lagrangian step rather than the fluid velocities used in the LEHGC code. The model is verified by several examples.