This paper presents the development of numerical method for unsteady seepage flow problems with phreatic surface. This method is new discrete procedure that applied the concept of hybrid model. Also, present method introduced the general weighting residual procedure with the subsidiary conditions. Compatibility of the potential on the intersection boundary is approximately introduced using penalty function instead of the Lagrange multiplier. In each element, we assume linear potential field with the head and the hydraulic gradient. The element node does not need to have the degree of freedom and is used only for recognition of the domain shape. In this study, we tried to analysis the unsteady seepage flow problems with phreatic surface to compare the experiment.