By extending Darcy’s law to dry domain above the free surface and specifying a Signorini’s type of boundary condition on the potential seepage boundaries

a new initial-boundary value problem of entire domain is deduced from non-steady seepage flow problems with free surface of in fractured rock masses. A parabolic variational inequality(PVI) formulation equivalent to the partial differential equation(PDE) is proposed for reducing the difficulty in selecting trial functions

while the finite element discretization and iterative algorithm of the PVI formulation are detailed. The validity of the numerical approach is then verified by existing laboratory tests. Finally

as an application to the non-steady seepage flow analysis in rock slope including complex fracture networks

the calculation results well illustrate the variations of free surfaces during drawdown of reservoir water level; characteristics of seepage flow movement in fracture networks and inhomogeneity of flow distribution at seepage surfaces.