This study presents a data-driven spatial interpolation algorithm based on physics-informed graph neural networks used to develop a thermal Earth model for the conterminous United States. The model was trained to approximately satisfy Fourier’s Law of conductive heat transfer by simultaneously predicting subsurface temperature, surface heat flow, and rock thermal conductivity. In addition to bottomhole temperature measurements, we incorporated other spatial and physical quantities as model inputs, such as depth, geographic coordinates, elevation, sediment thickness, magnetic anomaly, gravity anomaly, gamma-ray flux of radioactive elements, seismicity, electrical conductivity, and proximity to faults and volcanoes. With a spatial resolution of 18km2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$18 \\ km^2$$\\end{document} per grid cell, we predicted heat flow at surface as well as temperature and rock thermal conductivity across depths of 0-7km\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0-7 \\ km$$\\end{document} at an interval of 1km\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$1 \\ km$$\\end{document}. Our model showed temperature, surface heat flow and thermal conductivity mean absolute errors of 6.4∘C\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$6.4^\\circ C$$\\end{document}, 6.9mW/m2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$6.9 \\ mW/m^2$$\\end{document} and 0.04W/m-K\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$0.04 \\ W/m-K$$\\end{document}, respectively. This thorough modeling of the Earth’s thermal processes is crucial to understanding subsurface phenomena and exploiting natural underground resources. Our thermal Earth model is available as web application at https://stm.stanford.edu, feature layers on ArcGIS at https://arcg.is/nLzzT0, and tabulated data on the Geothermal Data Repository at https://gdr.openei.org/submissions/1592.