Abstract: Since the shape functions derived from the partition of unity-based meshless method, such as the element-free Galerkin method, are free of the Kronecker delta property, there are great troubles in the exact imposition of the essential boundary condition and boundary continuity of materials. Nevertheless, if adopting the penalty method or the Lagrange multiplier method, problems, like the selection of proper penalty factor and the satisfaction of the inf-sup condition, will occur. This study utilizes the property of partition of unity that once the local solutions satisfy some condition, the global solution will automatically satisfy the same condition. By constructing local approximations in physical patches of different types according to the boundary condition, a new moving least square interpolation-based numerical manifold method（MLS-NMM） is developed. Through the solution of some typical seepage flow problems, it is demonstrated that the proposed procedure is capable to deal with the problems of the singular angular point precisely and may provide an alternative solution for the seepage analysis in engineering.