In this paper we describe the finite element codes for two three-dimensional models of reactive transport in groundwater. The transport models are linear and include advection, dispersion, decay, and sorption. Sorption under both equilibrium (LEA3D code) and nonequilibrium (NONLEA3D) conditions is considered. In the latter case a "dual porosity" concept is used that subdivides the porous medium into five distinct interacting regions. The coupled model for nonequilibrium sorption is solved using an integro- differential approach. The large sparse systems of equations generated by the finite element discretization are nonsymmetric and are solved using efficient preconditioned conjugate-like methods. Tetrahedral elements and linear basis functions are used for the discretization in space, and a weighted finite difference formula is used for the discretization in time. The code handles: temporally and spatially variable boundary conditions of Dirichlet, Neumann, or Cauchy type; heterogeneous material and solute properties, including dispersivities,porosities,distribution coefficients, mass transfer coefficients, and decay constant; and both saturated and unsaturated flow regimes. The velocity and water saturation values need as input can be defined by the user or calculated as output from the FLOW3D flow code, wich is a companion paper.