Numerical simulation has become a widely practiced and accepted technique for studying flow and transport processes in the vadose zone and other subsurface flow systems. This article discusses a suite of codes, developed primarily at Lawrence Berkeley National Laboratory (LBNL), with the capability to model multiphase flows with phase change. We summarize history and goals in the development of the TOUGH codes, and present the governing equations for multiphase, multicomponent flow. Special emphasis is given to space discretization by means of integral finite differences (IFD). Issues of code implementation and architecture are addressed, as well as code applications, maintenance, and future developments.