Abstract <p>Consolidation equations of unsaturated soil can achieve a closed-form solution when the definite conditions (boundary and initial conditions) are given. The boundary conditions are also one of the factors that affect the consolidation behavior of unsaturated soils. To accurately predict the consolidation behavior of double-layered unsaturated soils under three types of boundary conditions, this paper proposes an efficient and stable implicit difference method. The proposed solution is based on the implicit difference method and virtual grid method to discretize the consolidation equation and boundary conditions, respectively. The stability of the numerical solutions is proved using the von Neumann criterion. Furthermore, the accuracy of the difference method is validated by comparing the results of the numerical solutions and analytical solutions. According to the validated numerical solutions, the consolidation behavior of double-layer unsaturated soil under five representative boundary conditions has been explored. The results show that the boundary conditions will significantly affect the dissipation rate of excess pore pressures and further change the ground settlement.</p>