To better characterize the intricate coupled thermo-hydro-mechanical dynamic (THMD) response in twodimensional saturated soil and to enrich the research object of Green-Naghdi (G-N) generalized thermoelastic theory, this study innovatively combines the G-N generalized thermoelastic theory and Caputo's fractional order derivative, to obtain the new control equations, and to establish a new fractional order thermoelastic theoretical model. The article is solved by the normal mode analysis (NMA), which can eliminate the integration error and solve the complex fractional order partial differential control equations quickly at the same time. The effects of different boundary conditions of fractional order derivatives, porosity, frequency, and thermal conductivity coefficients on non-dimensional excess pore water pressure, temperature, vertical displacement, and vertical stress are also fully analyzed, and the distribution curves of high precision numerical solutions are given. The results show that the effect of frequency variation on each non-dimensional variable is obvious. The effects of fractional order derivatives, porosity and thermal conductivity coefficients on the non-dimensional variables vary depending on the boundary conditions. The results provide theoretical support for geotechnical and environmental engineering.