The rate at which contaminants in groundwater undergo sorption and desorption in routinely described using diffusion models. Such approaches, when incorporated into transport models, lead to large systems of coupled equations, often nonlinear. This has restricted applications of coupled transport and diffusion models to one-dimensional systems. Further, numerical difficulties inherent in many common solution formulations to coupled transport and diffusion problems result in inaccurate and unreliable solutions to problems of common interest. The objective of this work was to develop methods that provide accurate and robust solutions to coupled transport and diffusion problems for single and multicomponent solute systems in both one and two spatial dimensions. Problems involving pore, surface, and combined pore and surface diffusion into spherical particles are considered. A split-operator formulation is proposed in which the reaction operator is separated from the transport operator and solved independently. Models derived from such split-operator formulations to selected coupled transport and diffusion problems are shown to be robust, accurate, and computationally efficient. The developed split-operator approaches are also amenable to solution using parallel, or for some problems massively parallel, processing methods. ¿ American Geophysical Union 1993