The　major　difficulty　in　the　analysis　of　unconfined　flow　in　porous　media　is　that　the　free　surface　is　unknown　a　priori,　where　the　nonlinearity　is　even　stronger　than　the　unsaturated　seepage　analysis.　There　is　much　space　for　both　the　adaptive　mesh　methods　and　the　fixed　mesh　methods　to　improve.　In　this　study,　firstly　two　variational　principles　fitted　to　the　numerical　manifold　method　(NMM)　are　formulated,　each　of　which　enforces　the　boundary　conditions　and　the　material　interface　continuity　conditions.　In　the　setting　of　the　NMM　together　with　the　moving　least　squares　(MLS)　interpolation,　then　the　discretization　models　corresponding　to　the　variational　formulations　are　built,　which　are　utilized　to　locate　the　free　surface　and　scrutinize　the　computational　results　respectively.　Meanwhile,　a　novel　approach　is　developed　to　update　the　free　surface　in　iteration.　With　high　accuracy　and　numerical　stability　but　no　need　to　remesh,　the　proposed　procedure　is　able　to　accommodate　complicated　dam　configuration　and　strong　non-homogeneity,　where　internal　seepage　faces　may　develop,　a　seldom　touched　problem　in　the　literature.　(C)　2014　Elsevier　Inc.　All　rights　reserved.