The　virtual　element　method　(VEM)　demonstrates　excellent　convergence　when　dealing　with　general　polygonal　mesh　shapes　and　addressing　nonlinear　problems.　Saturated–unsaturated　seepage　represents　a　highly　nonlinear　problem,　and　generating　appropriate　meshes　poses　a　challenge　in　seepage　studies.　This　paper　utilizes　the　VEM　to　discretize　the　Richards　equation　by　employing　the　hydraulic　head　formulation　and　introduces　computational　methods　for　the　corresponding　stiffness　matrix　and　mass　matrix.　Five　numerical　examples　are　employed　to　validate　the　effectiveness　and　feasibility　of　the　algorithm　in　solving　saturated–unsaturated　seepage　problems.　In　the　verification　part,　the　accuracy　of　the　algorithm　and　the　convergence　of　the　mesh　shape　are　verified　by　comparing　the　analytical　solution　of　the　saturation-unsaturated　seepage　problem　in　homogeneous　soil　with　the　classical　experimental　results.　Then,　the　effectiveness　of　the　proposed　algorithm　in　solving　the　saturation-unsaturated　seepage　problem　of　layered　soil　is　verified　by　comparing　with　the　analytical　solution　and　other　research　results.　Finally,　an　engineering　example　is　given　to　illustrate　the　feasibility　of　the　algorithm　for　solving　engineering　problems.　The　simulation　results　demonstrate　that　the　VEM　can　effectively　simulate　saturated–unsaturated　seepage　in　homogeneous　and　layered　soils　and　exhibits　strong　convergence.　©　2024　Elsevier　Ltd