An Eulerian-Lagrangian localized adjoint method (ELLAM) numerical solution is developed for the multiphase contaminant transport equations in two dimensions. The ELLAM uses finite volume test functions in the space-time domain defined by the characteristics of the hyperbolic part of the governing equation. The use of the characteristics results in an approximation that allows large time steps while still maintaining accurate solutions. This greatly reduces the computational effort required to find solutions to the advection-dispersion equation. Combination of the finite volume test functions with the conservative form of the governing equation results in a local conservation of mass property. For problems with a constant saturation such as that of contaminant transport in the saturated zone this property is highly advantageous as the resulting numerical approximation conserves mass both globally and locally. For problems with a variable saturation the requirement of local conservation of mass is too stringent, particularly around first type boundaries where oscillations can occur because of errors inherent in the numerical determination of fluid velocities and in the backtracking routine. A combined conservative/nonconservative ELLAM is developed with an ELLAM formulation based on the nonconservative form of the governing equation being applied to subdomains intersecting first-type boundaries and a conservative ELLAM being used for all other subdomains. The combined conservative/nonconservative ELLAM is compared to a Galerkin finite element scheme and is found to have greatly superior performance, requiring far fewer time steps to obtain a solution of equivalent accuracy. ¿ American Geophysical Union 1996.