This paper derives analytical solutions for steady-state one-dimensional (1-D), two-dimensional (2-D), and three-dimensional (3-D) two-phase immiscible subsurface flow for a counter-current problem. Since the governing equations are highly nonlinear, 2-D and 3-D derivations are generally difficult to obtain. The primary purpose for the solutions is to test finite difference/volume/element computer programs for accuracy and scalability using architectures ranging from PCs to parallel high performance computers. This derivation is accomplished by first solving for the saturation of water in terms of a function that is a solution to Laplace's equation to achieve a set of partial differential equations that allows some degree of latitude in the choice of boundary conditions. Separation of variables and Fourier series are used to obtain the final solution. The test problem consists of a rectangular block of soil where specified pressure is applied at the top and bottom of the sample, and no-flow boundary conditions are imposed on the sides. The pressure at the top of the sample is a step function that allows the testing of adaptive meshing or concentration of grid points in action zones. (c) 2008 Published by Elsevier Ltd.