During the process of treating soft soil foundations with prefabricated drainage drains (PVD), “soil columns” form around the PVD, and a “weak zone” forms outside the range of the “soil columns.” The difference in properties between the two forms a distinct interface, leading to a gradual decrease in drainage efficiency and obstruction of vertical drainage channels, which in turn causes cracks and lateral displacement in the soil during consolidation. The interfaces between adjacent soil layers are incomplete contact, and the water within the interstices impedes the transfer of heat, manifesting a thermal resistance effect. To address this phenomenon, a synchronous measurement system for the thermal gradient and the heat flux density between the soil interfaces has been developed. Applying Fourier's law of heat conduction, the thermal resistance coefficient has been determined. Based on the theory of thermo‐hydro‐mechanical coupling, a multi‐zone axisymmetric model for saturated soils that considers thermal resistance effect has been proposed. Semi‐analytical solutions were derived and validated through comparison with the custom FEM model and field experiments. The thermal consolidation characteristics of the multi‐zone soils under various thermal contact models have also been discussed, with a comprehensive analysis of the influence of different parameters. Outcomes show that: the generalized incomplete thermal contact model provides a better description of the thermal resistance phenomenon between multi‐zone soils interfaces; ignoring the thermal resistance effect leads to an overestimation of the deformation during the thermal consolidation, and, the thermal resistance effect decreases the influence of the thermo‐osmosis effect on the consolidation characteristics.