Abstract: To solve non-steady seepage flow problems with free surface in fractured rock masses, Darcy's law is extended to the entire domain. A parabolic variational inequality (PVI) formulation, in which a Signorini's type of boundary condition enforced on the potential seepage surface, is established for transforming the flux condition on the potential seepage surface into natural boundary conditions, and then proved to be equivalent to the partial differential equation (PDE) formulation, and then the difficulty in solving this problem is reduced. Finite element numerical solution of the PVI formulation is proposed, and the validity of the numerical approach is verified by comparison of theoretical and calculated results for cross fracture model. Finally, the proposed approach is applied for non-steady seepage flow behaviors in a complex fractured rock slope, and the calculation results validate the reliability and robustness of this method well.