Abstract: In Getechnical engineering, such as slope, underground tunnel, dam and dike, and soil foundation etc, 3-D unsteady seepage problem often occurs.In performing the numerical analysis of unsteady seepage problem, it is necessary to calculate the integral on free surface boundary.In the paper, a high precision numerical integral method is presented.Based on the 8-node isoparametric element, the equat ion of the free surface boundary in an element is obtained firstly according to zero water pressure condition, and then the integral on the free surface boundary is translated into the integral on the plane ξη of local reference.And the integral region is calculated according to the projections of the intersecting points where free surface and the boundaries of the element cross.Finally the whole integral region is divided into several triangle sub regions, and the Simpson method with changing steps is performed to calculate the integral on each triangle sub-region.The high precision numerical integral method can avoid the assumption that the free surface is a plane within an element.It can improve the calculation accuracy, especially for the case that the free surface within an element varies very sharply The applications in both the sand tank model and an example of homogeneous earth dam are given, it shows that the method is practical.