Seepage problems with complex drainage systems are commonly encountered in civil engineering, with strong non-linearity. A numerical solution based on the Finite Element Method combining the substructure technique with a variational inequality formulation of Signorini's type is proposed to solve these problems. The aims of this work are to accurately characterize the boundary conditions of the drainage systems, to reduce the difficulty in mesh generation resulting from the drainage holes with small radius and dense spacing, and to eliminate the singularity at the seepage points and the resultant mesh dependency. Numerical stability and robustness of the proposed method are guaranteed by an adaptive procedure for progressively relaxing the penalized Heaviside function associated with the formulation of the discrete variational inequality. Two challenging numerical examples are presented to validate the effectiveness and robustness of the proposed method. (c) 2007 Elsevier Ltd. All rights reserved.