This study presents a data-driven spatial interpolation algorithm based on graph neural networks to develop national temperature-at-depth maps for the conterminous US. In addition to bottomhole temperature measurements, we incorporated other physical quantities, such as depth, geographic coordinates, elevation, sediment thickness, magnetic anomaly, gravity anomaly, and gamma-ray flux produced by radioactive element. We constructed subsurface temperature predictions for depths of 0-7 km at an interval of 1 km with spatial resolution of 18 km2 per grid cell. This thorough modeling of subsurface temperature is crucial to understanding subsurface phenomena and exploiting natural underground resources.