Based on the idea of the initial flow method

initial velocity is introduced to eliminate the velocity that does not exist in dry zone. Darcy formula is extended to the entire region and a nonlinear boundary value problem of the entire region is defined. Through categorizing the potential seepage face as a boundary condition of Signorini’s type

an equivalent variational inequality is established. By combining the continuous Heaviside function

a finite element iterative algorithm of the variational inequality based on initial flow method is given. With this algorithm

a program is developed. Two typical examples are employed to demonstrate the efficiency and robustness of the algorithm proposed here in locating the free surface of complex fracture networks subject to underground water.