An analytical method is proposed to solve the contaminant transport in a double-layered contaminated soil with imperfect diffusion boundaries. By virtue of the separation of variables scheme, the governing equation for contaminant diffusion is split into two ordinary differential equations, and the corresponding general solutions are obtained. By utilizing the initial, continuity and boundary conditions of the system, and the orthogonality of trigonometric functions, the analytical solution for contaminant concentration in the double-layered contaminated soil is derived, which can be further used to describe the average degree of diffusion. The reliability and accuracy of the developed solution is verified by comparing with the existing analytical solution and numerical simulation results. Selected numerical examples are further presented to analyze the influence of the imperfect diffusion boundaries on the spatial/time distribution of contaminant concentration, average degree of diffusion and mass flux. The results show that greater imperfect diffusion capacity coefficients lead to lower contaminant concentration distribution within the entire depth range and higher average degree of diffusion for a given time. Particularly, it only takes 23.32 years to complete the entire diffusion process when imperfect diffusion capacity coefficients are infinite. The developed solution can provide useful guidance for engineering practice.