Abstract: This paper presents a wave functions method to calculate the three dimensional dynamic interaction between a circular tunnel and the surrounding unsaturated foundation soil. The tunnel lining is conceptualized as an infinite Flügge cylindrical shell, and the unsaturated foundation soil is conceptualized as a three-phase medium consisting of fluid, solid and gas. The motion equations of infinite Flügge cylindrical shell are solved by the separation of variable method, and the Helmholtz decomposition theorem is used to solve the unsaturated soil governing equations. Based on the boundary conditions of displacement, stress and pore fluid pressure at the unsaturated soil-tunnel interface and the ground surface, the dynamic interaction between the unsaturated soil and circular tunnel is solved when the tunnel invert subjected to a unit harmonic load. The transformation properties of plane wave functions and cylindrical wave functions are used to apply the boundary conditions expressed in both the rectangular and cylindrical coordinate systems before the coupling. The proposed method is verified by comparing with the results obtained by the existing 2.5 dimensional coupled FE-BE method for single phase elastic medium, 2.5 dimensional coupled FE-BE method for two phase saturated porous medium and Pip in Pip semi analytical method for three-phase unsaturated medium. Finally, by using the proposed method, the dynamic response of the unsaturated soil-tunnel system under different soil saturations is analyzed by a case study. The results show that saturation has a great effect on the dynamic response amplitude of soil displacement and excess pore water pressure. Because the three-phase medium simulating unsaturated soil can be reduced to two-phase medium or single-phase medium, this method can also be used to calculate and analyze the dynamic response of a circular tunnel buried in two-phase saturated soil or in single-phase elastic soil after the degeneration of the calculation parameters of unsaturated soil.