This paper proposes the solution of seepage problems with free surface. It is obtained through meshless Galerkin method. The corresponding basic equations and integration schemes are derived. Applications of penalty function method to the processing procedure of seepage boundary condition and free surface boundary are also put forward. Two-dimensional meshless method software LIDAREFM is developed with IDL language. Research on seepage field characteristics and distribution of free surface under different discontinuity combinations are carried out, taking the fractured rock slope at Qiaozi Town, Huairou County in Beijing City as an example. The results show that the meshless method has certain advantages in solving seepage field of intensive fracture problem. The layout independence of gridding nodes on integration points is achieved. The re-processing of free surface and sub-domain with cracks crossing is avoided. The process of solving the seepage problem is simplified.