Abstract: The discontinuities such as cracks, joints and faults can have the varying degrees of developments. The permeability of fractured rock can show the characteristics of anisotropy and discontinuity. The traditional finite element method encounters difficulty in the solution of the seepage field for intensive fractured rocks. 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 as an example. The seepage field features under different degrees of opening and connectivity are discussed. 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.