It is of practical interest to investigate the natural evaporation of volatile organic compounds (VOCs) after the removal of a leaking tank situated on the top of the soil. This study aims to develop a mathematical model to predict mole fraction distributions and migration of evaporation front for two VOCs emanating from residual non-aqueous phase liquid (NAPL) due to the leak from the tank in a homogeneous soil. Considering the location of the front and the regions above and below the front, a numerical model for the diffusive transport of VOCs in unsaturated soils was developed using the finite difference method with a moving grid approach. The model was further simplified to the case of single VOC and solved analytically by Boltzmann's transformation with a moving boundary. Analytical expressions for the depth and moving speed of the front for a single VOC were then obtained for practical use. Finally, the developed model was used to predict the concentration distributions of VOCs below the land surface and examine the factors affecting the location and moving speed of the evaporation front.