[1]Mohammad J. Aljubran

[1]Energy Science & Engineering, Stanford University, Building 367, Stanford,
94305, California, United States

# Thermal Earth Model for the Conterminous United States Using an
Interpolative Physics-Informed Graph Neural Network (InterPIGNN)

<EMAIL_ADDRESS>Roland N. Horne<EMAIL_ADDRESS>*

###### Abstract

This study presents a data-driven spatial interpolation algorithm based on
physics-informed graph neural networks used to develop national temperature-
at-depth maps for the conterminous United States. The model was trained to
approximately satisfy the three-dimensional heat conduction law by
simultaneously predicting subsurface temperature, surface heat flow, and rock
thermal conductivity. In addition to bottomhole temperature measurements, we
incorporated other physical quantities as model inputs, such as depth,
geographic coordinates, elevation, sediment thickness, magnetic anomaly,
gravity anomaly, gamma-ray flux of radioactive elements, seismicity, and
electric conductivity. We constructed surface heat flow, and temperature and
thermal conductivity predictions for depths of $0$-$7\ km$ at an interval of
$1\ km$ with spatial resolution of $18\ km^{2}$ per grid cell. Our model
showed superior temperature, surface heat flow and thermal conductivity mean
absolute errors of $4.8^{\circ}C$, $5.817\ mW/m^{2}$ and $0.022\ W/(C\cdot
m)$, respectively. The predictions were visualized in two-dimensional spatial
maps across the modeled depths. This thorough modeling of the Earth’s thermal
processes is crucial to understanding subsurface phenomena and exploiting
natural underground resources. The thermal Earth model is made available as
feature layers on ArcGIS at https://arcg.is/nLzzT0.

###### keywords:

Temperature-at-Depth; Heat Flow; Rock Thermal Conductivity; InterPIGNN;
Physics-Informed; Graph Neural Networks

##  1 Introduction

Understanding temperature variations across depths and geographies holds
significance to many domains. Temperature is fundamental to modeling heat
distribution, thermal anomalies, groundwater flow and management, volcanic
eruptions, stresses and earthquakes, hydrocarbon maturity, carbon
sequestration, subsurface energy storage, amongst others [1, 2, 3, 4]. In
particular, accurate modeling of temperatures at significant depths within the
Earth’s crust is central to the exploration and development of geothermal
energy. A conventional geothermal resource requires three ingredients: (1)
heat source, (2) fluid for heat transport, and (3) permeable/fractured rock
for fluid flow. Other unconventional geothermal systems can also be viable
based on site-specific circumstances [5, 6]. Specifically, enhanced geothermal
systems (EGS) only require sufficiently high-temperature resources as heat
transport and fluid flow can be achieved through water injection and rock
stimulation, respectively [7]. Many recent efforts were focused on assessing
the EGS resource potential across the United States [8, 9, 10, 11] using
either regional or national temperature-at-depth maps.

Regional temperature-at-depth maps were developed for various locations. One
study by researchers at Southern Methodist University (SMU) Geothermal
Laboratory modeled heat flow of the Appalachian Basin using two-dimensional
finite difference heat conduction, basin cross sections, equilibrium
temperature, and oil and gas bottom-hole temperature (BHT) data to quantify
heat flow at the surface and at the base of the sedimentary basin. Other
efforts by researchers at Cornell University used heat flow estimates based on
BHT measurements to estimate temperature-at-depth at the Appalachian Basin
using a one-dimensional heat conduction model [12] coupled with Monte Carlo
analysis of various physical quantities for uncertainty quantification. SMU
researchers also developed a regional temperature-at-depth model for the
Cascades using a three-dimensional finite difference heat conduction model
[13]. Other thermal models were focused on the Great Basin and led by multiple
institutes, e.g., Nevada Bureau of Mines and Geology, US Geological Survey
(USGS), and INGENIOUS project [14, 15, 16]. Others estimated temperature-at-
depth and subsurface thermal properties at the Snake River Plain using field
measurements alongside a conductive and radiogenic heat transfer and
production model. [17, 18]. Whereas such regional models are useful, they are
spatially limited and vary in input quantities, assumptions, BHT accuracy,
spatial and depth-wise resolutions, inconsistency if integrated, amongst other
aspects. This prompts the need for a national temperature-at-depth models to
guide conventional and non-conventional geothermal exploration and
development.

National temperature-at-depth maps for the conterminous US were developed by
the SMU Geothermal Laboratory as part of a 2006 study on The Future of
Geothermal Energy, mainly focused on assessing the potential of EGS across the
US [19, 11, 8]. However, the SMU maps have several limitations:

  * •

They were developed at a gridding interval of 5 minutes (or,
$0.08333^{\circ}$) of spatial resolution, which translates to an average
spacing of about 8 km representing an area of about $64\ km^{2}$ per grid
cell. A typical 250 MWe EGS plant might require about $10-20$ $km^{2}$ of
reservoir planar area to accommodate the thermal resource needed, assuming
that heat removal occurs in a 0.5 km-thick region of hot rock at depth [11].
With such a large areal extent, a $64\ km^{2}$ grid cell is likely to filter
out many local heat anomalies. These maps also ignore EGS resources shallower
than 3 km entirely, and only model depth intervals between 3-10 km [9].

  * •

The analytical approach used in developing these temperature-at-depth maps
requires knowledge of heat flow across the conterminous US. Note that all
hydrothermal system-influenced data (very high values, i.e., generally greater
than $120\ mW/m^{2}$ were conveniently excluded in this process to avoid the
difficulty involved in analytically modeling regions with hot fluid upflow
without overestimating temperatures at neighbouring sites. This inconveniently
leads to smoothing many local heat anomalies which developers would rather
like to capture for economic EGS development.

  * •

In developing such heat flow maps, a minimum curvature algorithm was used.
However, such interpolation algorithms (1) are less effective with irregularly
spaced data points, and (2) tend to generate smooth fits that dismiss
potential heat anomalies in undersampled regions.

  * •

SMU adopted a two-layer model where the sediment and basement contributions
were captured separately. This requires knowledge of either
unavailable/assumed or interpolated spatial properties across layers, such as
rock thermal conductivity, measured sediment heat flow, basement heat flow,
radioactive depth variable, and radioactive heat generation.

Interested in estimating shallow, low-temperature geothermal resources across
the conterminous US, researchers at the National Renewable Energy Laboratory
developed a temperature-at-depth model to capture depths shallower than $3\
km$ [20]. Approximately 300,000 BHT measurements were used in that study. To
reduce the nugget effect associated with the variogram of the kriging
algorithm used in that study, the BHT observations were calibrated to the
nearest $500\ m$ depth interval using geothermal gradients obtained from the
$3.5$-$km$ temperatures originally predicted by the SMU temperature-at-depth
maps. Using a grid of $16\ km^{2}$, BHT measurements were sampled down to
approximately 71,000 records per grid cell using the median statistic. Then,
they fit a linear regression model to predict BHT using (1) depth, and (2) the
$3.5$-$km$ temperatures originally predicted by the SMU temperature-at-depth
maps. Next, they applied ordinary kriging to the linear regression residual
which were used to improve local fit. In addition to the aforementioned
caveats associated with the SMU temperature-at-depth maps, these shallower
maps lack spatial and depth-wise resolution due to downsampling real BHT
measurements. Our investigation shows that subtracting the $3$-$km$
temperature-at-depth map developed in this study from the $3.5$-$km$
temperature-at-depth map developed by SMU results in an average difference of
over $-50^{\circ}\ C$.

This study presents an Interpolative Physics-Informed Graph Neural Network
(InterPIGNN) used to develop national temperature-at-depth maps for the
conterminous US. In addition to bottomhole temperature measurements, we
incorporated other physical quantities as model inputs, such as depth,
geographic coordinates, elevation, sediment thickness, magnetic anomaly,
gravity anomaly, gamma-ray flux of radioactive elements, seismicity, and
electric conductivity. Data curation, preprocessing, and normalization are
significant steps in this process. We trained, validated, and tested our model
on unseen BHT measurements to evaluate generalizability. We then used it to
predict temperature-at-depth for depths of $0$-$7\ km$ at an interval of $1\
km$ with spatial resolution of $18\ km^{2}$ per grid cell that is sufficient
for typical EGS projects.

##  2 Data Collection

Various physical quantities and data sources are required toward an accurate
estimate of temperature-at-depth across the conterminous US. Measurements span
various geographic locations (Easting/Northing) and/or depth. We collected
measurements for three thermal quantities: BHT, heat flow, and rock thermal
conductivity. Additionally, we allocated measurements with sufficiently high
spatial resolution of other physical quantities which can be explicitly or
implicitly correlated with temperature-at-depth, i.e., average surface
temperature, elevation, sediment thickness, magnetic anomaly, gravity anomaly,
Uranium radiation, Thorium radiation, Potassium radiation, seismicity, and
electric conductivity. Note that these measurements have various spatial
resolutions, hence we downsampled/upsampled them to populate the continental
United States with the target resolution of $18\ km^{2}$ using inverse
distance weighing [21]. Note that all visualizations reported in this article
were generated by the authors accordingly. In the following, we briefly
describe the source and distribution associated with these quantities.

###  2.1 Thermal Quantities

BHT measurements are acquired through direct wellbore temperature measurements
at one or more depths. The collected temperature measurements could be
associated with geothermal, oil & gas, groundwater, observation, and
monitoring wells. After exhaustive search and curation of BHT measurements, we
identified nine different raw and aggregated data sources which we integrated
into a single database. These sources are SMU, Association of American State
Geologists (ASSG), USGS, Utah Geological Survey (UGS), Colorado Geological
Survey (CGS), Maine Geological Survey (MGS), Washington State Department of
Natural Resources (WSDNR), Geothermal Information Layer for Oregon (GTILO),
and Great Basin Center for Geothermal Energy (GBCGE).

Due to the temperature differences between the rock and the wellbore fluid,
those BHT measurements have to be corrected [22]. Due to mud circulation and
borehole temperature disequilibrium, this thermal disturbance is not equally
distributed along depth. Consequently, compared to the adjacent rock, deeper
and shallower sections indicate lower and higher temperatures, respectively
[23]. To correct BHT measurements, we used the Harrison correction which is
based on a second-order polynomial equation developed by the Oklahoma
Geological Survey [24]. We found this method to be more suitable to this study
as it can be used with only temperature and depth measurements available.
Meanwhile, duplicated BHT measurements across the different data sources were
detected and filtered out. Note that records are considered to be duplicates
if they have the same depth, Easting, and Northing when rounded to the nearest
0.1 m, 0.1 km, and 0.1 km, respectively. As shown in Fig 1, our database
included a total of 400,134 unique BHT measurements with maximum and minimum
corrected BHT measurements of $0.16$ and $533.71^{\circ}\ C$, respectively
(going forward, we will drop the term “corrected” for brevity). Figs 2 and 3
show all data points projected on the conterminous US and colored based on
source and BHT, respectively.

Figure 1: Log-scale plot of BHT record count per data source. Figure 2: BHT
records projected on the conterminous US, colored based on data source. Figure
3: BHT records projected on the conterminous US, colored based on magnitude.
Note that we set 250${}^{\circ}\ C$ as an upper threshold for visual
convenience.

The SMU dataset also aggregated data for heat flow and rock thermal
conductivity [19]. As seen in Fig 4, this dataset indicates heat flow
estimates as low as nearly zero and as high as $23,012\ mW/m^{2}$ at
Yellowstone, with median of $59\ mW/m^{2}$. Generally, we observed several
anomalous heat flow estimates around calderas, such as Yellowstone, Newberry,
and Valles. Meanwhile, records indicate rock thermal conductivity values
ranging between $0.33$ and $7.24\ W/m-K$, with median of $2.13\ W/m-K$, as
seen in Fig 5.

Figure 4: Heat flow records for the conterminous US, colored based on
magnitude. Note that we set $300\ mW/m^{2}$ as an upper threshold for visual
convenience. Figure 5: Thermal conductivity records for the conterminous US,
colored based on magnitude.

Near-surface temperature data is based on average 40-inch soil temperatures
reported by the National Weather Service. As seen in Fig 6, the average
surface temperature ranges between $-7$ and $25^{\circ}\ C$, with mean and
standard deviation of $6.5\pm 7.6^{\circ}\ C$.

Figure 6: Surface temperature records for the conterminous US, colored based
on magnitude.

###  2.2 Elevation and Sediment Thickness

Additionally, we gathered high-resolution maps depicting elevation and
sediment thickness across the conterminous US. Elevation data is queried from
the National Map’s Elevation Point Query Service offered by the USGS, which
has an overall resolution of 0.53 meters [25]. As seen in Fig 7, elevation
ranges between $-84$ and $4,192\ m$, with mean and standard deviation of
$791\pm 726\ m$. Sediment thickness data were acquired through the Oak Ridge
National Laboratory Distributed Active Archive Center. This dataset provides
estimates of the thickness of the permeable layers above bedrock with spatial
resolution of $1\ km^{2}$ based on input data for topography, climate, and
geology. These data are modeled to represent estimated thicknesses by landform
type for the geological present [26]. Seen in Fig 9, this dataset has
resolution up to $50\ m$, where values of $50\ m$ indicate sediment thickness
of $\geq 50\ m$ while lower values are indicative of exact value estimated by
the data source. To emphasize the basin-and-ridge scale structural features,
we created detrended maps for both elevation and sediment thickness, seen in
Figs 8 and 10. We followed a simple procedure using a $30\ km\ \times\ 30\ km$
moving window to compute a trend map for each quantity which is then
subtracted from the corresponding raw map to finally arrive at the detrended
map.

Figure 7: Elevation records for the conterminous US, colored based on
magnitude. Figure 8: Detrended elevation map for the conterminous US, colored
based on magnitude. Figure 9: Sediment thickness records for the conterminous
US, colored based on magnitude. Values of $50m$ indicate sediment thickness of
$\geq 50m$ while lower values are indicative of exact value estimated by the
data source Figure 10: Detrended sediment thickness map for the conterminous
US, colored based on magnitude.

###  2.3 Magnetic and Gravity Anomalies

Magnetic and gravity anomalies were acquired through the USGS Mineral
Resources Online Spatial Data. This data is a result of a joint effort between
the USGS, the Geological Survey of Canada, and Consejo de Recursos Minerales
of Mexico [27]. The airborne measurement of the earth’s magnetic field over
North America describes the magnetic anomaly caused by variations in earth
materials and structure. We accessed the corrected aeromagnetic grid data with
a resolution of $1\ km^{2}$ per grid cell. Seen in Fig 11, magnetic anomalies
for the conterminous US range between $-1,513$ to $2,614\ nanoteslas$, with
mean and standard deviation of $41\pm 206\ nanoteslas$. Bouguer gravity
anomaly data was based on National Information Mapping Agency data files which
were gathered and curated by the USGS with a spatial resolution of nearly $10\
km^{2}$ per grid cell. Gravity anomalies are produced by density variations
within the rocks of the Earth’s crust and upper mantle. The free-air
correction reduces the measurement to sea level by assuming there is no
intervening mass as a uniform slab of constant density, and the Bouguer
correction includes the effects of constant density topography within 166.7 km
of the measurement location [28]. Seen in Fig 12, Bouguer gravity anomaly
ranges between $-224$ and $105\ milligal$, with mean and standard deviation of
$-8.9\pm 19.8\ milligal$.

Figure 11: Magnetic anomaly for the conterminous US, colored based on
magnitude. Note that we set $-1000$ and $1000\ nanoteslas$ as an lower and
upper thresholds, respectively, for visual convenience. Figure 12: Bouguer
gravity anomaly for the conterminous US, colored based on magnitude.

###  2.4 Radioactive Elements

Aerial gamma-ray surveys measure the gamma-ray flux produced by radioactive
decay of the naturally occurring elements K-40, U-238, and Th-232 in the top
few centimeters of rock or soil. With a spatial resolution of nearly $2.5\
km^{2}$, the USGS aggregated and curated aeroradiometric data for the
Conterminous United States which originate to the National Uranium Resource
Evaluation Program of the US Department of Energy. [29, 30]. Radioactive decay
of K-40, U-238, and Th-232 data is reported in the units of percent potassium
(% K), parts per million equivalent thorium (ppm eTh), and parts per million
equivalent uranium (ppm eU), respectively. Seen in Figs 13, 14, and 15,
radioactive decay of K-40, U-238, and Th-232 shows mean and standard deviation
of $1.17\pm 1.53\ \%\ K$, $1.72\pm 1.65\ ppm\ eU$, and $6.19\pm 3.64\ ppm\
eTh$, respectively.

Figure 13: K-40 radioactive decay for the conterminous US, colored based on
magnitude. Note that we set $5\ \%K$ as an upper threshold for visual
convenience. Figure 14: U-238 radioactive decay for the conterminous US,
colored based on magnitude. Note that we set $5\ ppmeU$ as an upper threshold
for visual convenience. Figure 15: Th-232 radioactive decay for the
conterminous US, colored based on magnitude. Note that we set $30\ ppmeTh$ as
an upper threshold for visual convenience.

###  2.5 Seismicity and Electric Conductivity

Seismic wave velocity and electrical conductivity are highly affected by the
subsurface rock and fluid properties, and in-situ conditions (e.g., porosity,
pore pressure, temperature, elastic moduli). Generally, with increasing in-
situ temperatures, seismic wave velocity and thermal conductivity decrease
[31, 32]. We acquired seismic and electric conductivity measurements nation-
wide using the USArray project, which deployed transportable seismic and
magnetic stations across the United States over years to passively capture the
Earth’s seismic activities and magnetic field. We particularly used inverted
and processed maps generated by various studies based on data captured by this
project, including the average geothermal gradient at the mantle [33], crustal
thickness [34], compressional-shear wave velocity ratio [35], velocity
perturbation [36], and electric conductivity [37]. Figs 16, 17, and 18
demonstrate some of these properties after we upsampled them to our target
spatial resolution of $18\ km^{2}$. Where properties vary with depth (e.g.,
seismic velocities, velocity ratio, electric conductivity), we additionally
normalized each quantity along depth to eliminate their naturally increasing
or decreasing trends with depth.

Figure 16: Average geothermal gradient at a mantle depth of $60\ km$ across
the continental United States as estimated using seismic measurements. Figure
17: Crustal thickness at a depth of $60\ km$ across the continental United
States as estimated using seismic measurements. Figure 18: Compressional wave
velocity at a depth of $2\ km$ across the continental United States as
estimated using seismic measurements.

##  3 Methods

Three-dimensional point clouds are typically unordered, sparse, and unevenly
distributed over an irregular grid. They are often associated with point-level
features but with various spatial resolutions (e.g., sensor data where
measurements are taken at different three-dimensional locations). We are
interested in interpolation tasks where a quantity of interest is to be
predicted for some desired irregular three-dimensional grid. For instance, our
study was interested in regressing temperature-at-depth with accuracy along
depth intervals and across the conterminous US.

Popular three-dimensional spatial interpolation techniques include nearest
neighbours [38], inverse distance weighing [21], kriging [39], radial basis
function kernels [40, 41], and spline methods [42, 43]. Each of these methods
suffer one or more caveats, e.g., discrete and abrupt surfaces in the presence
of dense spatial measurements, smoothing, dismissal of spatial structure
information and correlation, linearity, lack of sufficient non-linearity, use
of feature engineering and heuristics, and local trends [44, 45]. Multiple
deep learning architectures were proposed to treat point cloud inputs on
irregular grids for various tasks, including classification, regression,
segmentation, amongst others. They can be lumped into three categories: (1)
point-based methods (e.g., PointNet, PointNet++, PointWeb) [46, 47, 48], (2)
convolution-based methods (e.g., Relation-shape CNN, DensePoint, Pointwise-
CNN) [49, 50, 51], (3) registration methods (e.g., PPFNet33, FoldNet,
PointNetLK) [52, 53, 54, 55, 56, 57], and (4) graph-based methods (e.g., Edge-
Conditioned Convolution, EdgeConv, 3DGCN) [58, 59, 60]. Whereas these
architectures can input irregular three-dimensional point clouds, they are not
particularly designed for interpolation where the target quantity can
explicitly inform predictions at neighbouring grid cells. Additionally, the
presence of extensive sparsity in point cloud data makes it challenging to
generalize spatially. We treated these challenges using a novel Interpolative
Physics-Informed Graph Neural Network (InterPIGNN). The interpolative aspect
of this architecture incorporates a graph-based operator which, unlike
convolutional neural networks, does not require input data to lie on regular
grids. Graph-based architectures are capable of extracting point cloud
features without the need for heuristically projecting them on a regular grid.
Meanwhile, the physics-informed aspect involves satisfying physical laws
related to the problem of interest. In this section, we describe the
InterPIGNN module along with the data preprocessing and modeling steps. We
conveniently use the three-dimensional temperature-at-depth mapping task as an
example interpolation task to further motivate the InterPIGNN module.

###  3.1 Interpolative Convolutional Module

In three-dimensional point cloud interpolation tasks across space of interest,
$\Omega$, a target quantity $y$ is measured at a sparsely distributed finite
set of points, which we denote as the set of anchor nodes $V_{a}$ with
$|V_{a}|=n_{a}$. The objective is to predict or interpolate the target
quantity at a finite set of desired points (regular or irregular) sampled from
$\Omega$, which we denote as the set of grid nodes $V_{g}$ with with
$|V_{g}|=n_{g}$. We denote the union set of the anchor and grid nodes as
$V=V_{a}\bigcup V_{g}$ with $|V|=n_{a}+n_{g}=n$. Each node $v\in V$ is
associated with a position vector $\boldsymbol{p_{v}}\in\Omega$ (in three-
dimensional point clouds, $\Omega=\mathbb{R}^{3}$). We consider two sets of
features: (1) high-resolution features $F_{hr}$ which are known at all nodes
$v\in\mathcal{V}$, and (2) low-resolution features $F_{lr}$ which are only
known at anchor nodes $v\in V_{a}$. For a node $v\in V$, the high-resolution
and low-resolution feature vectors are denoted as $\boldsymbol{x}_{v,hr}$ and
$\boldsymbol{x}_{v,lr}$, respectively, where the latter is only defined when
$v\in V_{a}$.

We convert the point cloud into a directed graph representation $G=(V,E)$ by
constructing a set of directed edges $E\subseteq V_{a}\times V$ with $|E|=m$.
A directed edge $e_{uv}$ is constructed between nodes $u\in V_{a}$ and $v\in
V$ if $\boldsymbol{p}_{u}$ belongs to (1) the $k$-nearest neighbours of
$\boldsymbol{p}_{v}$, and/or (2) a sphere centered at $\boldsymbol{p}_{v}$
with radius $r$, where $k$ and $r$ are hyperparameters. An edge $e_{uv}$ has
edge weight $w_{uv}$ defined as the inverse euclidean distance and normalized
such that $\sum_{u\in\mathcal{N}(v)}w_{uv}=1$, where $\mathcal{N}(v)$ is the
set of neighbours of node $v$.

Following the GNN framework of message passing, aggregation, and combination,
we define our InterPIGNN convolutional module which can be used flexibly as a
nonlinear preprocesser or transformer in any deep learning algorithm inputting
point clouds. Consider an GNN architecture based on InterPIGNN with $L$ neural
layers. At layer $l=0$, the hidden state of each node $v\in V$ is initialized
using the high-resolution feature vector,
$\boldsymbol{h}_{v}^{(0)}=\boldsymbol{x}_{v,hr}$. We also assign each edge
$e_{uv}$ an attribute
$\boldsymbol{a}_{uv}^{(l)}=\boldsymbol{h}_{u}^{(l)}-\boldsymbol{h}_{v}^{(l)}$
that is updated across convolutional layers. For a target node $v\in V$, the
InterPIGNN module, seen in Equation 1, consists of four neural networks: (1)
root network $f_{\theta_{r}}$ to convolve root high-resolution features, (2)
neighbour network $f_{\theta_{n}}$ to transform each individual message passed
by neighbouring nodes including features, positions, edge attributes, and
labels, and (3) global network $f_{\theta_{g}}$ to transform aggregated
messages, and (4) post-processing network $f_{\theta_{p}}$ to jointly
transform the outputs of the root and global networks. Note that $\parallel$
represents the vector concatenation operator. Also, note that permutation-
invariant aggregation is achieved in the max pooling operation. The final
hidden $\boldsymbol{h}_{v}^{(L)}$ for each node $v\in V$ can then be used for
the desired task.

$\boldsymbol{h}_{v}^{(l+1)}=f_{\theta_{p}}\Bigg{[}f_{\theta_{r}}(\boldsymbol{h}_{v}^{(l)})+f_{\theta_{g}}\Bigg{(}\frac{1}{|\mathcal{N}(v)|}\sum_{u\in\mathcal{N}(v)}\Big{[}w_{uv}\cdot
f_{\theta_{n}}(\boldsymbol{x}_{u}\parallel\boldsymbol{p}_{u}\parallel\boldsymbol{a}_{uv}^{(l)}\parallel
y_{u})\Big{]}\Bigg{)}\Bigg{]}$ (1)

###  3.2 Three-Dimensional Heat Conduction

Temperature within the Earth’s subsurface is governed by a complex interplay
of physical processes, mainly thermal conduction, which drives heat transfer
through the subsurface materials. A significant contribution also comes from
the radioactive decay of isotopes such as uranium, thorium, and potassium,
which generates heat. In addition, convective movements of magma and
groundwater can redistribute heat vertically and horizontally through
advection, altering temperature profiles. The latent heat associated with
phase changes of water and tectonic activities, such as volcanic upwelling,
further influence subsurface temperatures [61]. Collectively, these factors
establish the thermal state of the subsurface, varying significantly with
geographic location, depth, and local geophysical and hydrological conditions.

Using the collected data for this study, we could satisfy the steady-state
three-dimensional heat conduction in the Earth’s subsurface which is
represented by a partial differential equation (PDE) as seen in Eqs. 2-3 [62,
63]. Subsurface temperature $T$, rock thermal conductivity $K$, crustal heat
flow $Q_{c}$, and local heat generation $\dot{g}$ vary with spatially
according to easting $x$, northing $y$ and depth $z$. Meanwhile, surface heat
flow $Q$ is only a function of $x$ and $y$.

$\frac{\partial}{\partial x}\Big{(}K\frac{\partial T}{\partial
x}\Big{)}+\frac{\partial}{\partial y}\Big{(}K\frac{\partial T}{\partial
y}\Big{)}+\frac{\partial}{\partial z}\Big{(}K\frac{\partial T}{\partial
z}\Big{)}+\dot{g}=0$ (2)

$Q_{c}(z)=Q-\int_{0}^{z}\dot{g}\ dz^{\prime}$ (3)

###  3.3 Physics-Informed Graph Neural Networks

Considering the unique BHT measurements, seen in Fig 3, the temperature-at-
depth interpolation task is a three-dimensional setting with $n_{a}=400,134$
anchor nodes. We aim to map temperature-at-depth from surface and down to $7\
km$ for the conterminous US with $18\ km^{2}$ and $1\ km$ spatial and depth-
wise resolutions, hence we have a total of $n_{g}=4,279,536$ grid points.
Edges were constructed following the $k$-nearest neighbours strategy with
$k=5$. The three-dimensional position vector $\boldsymbol{p}$ was defined
using depth, Easting, and Northing. Node labels $y$ were set as the BHT
measurement at anchor nodes, but they are unknown at grid nodes. Low-
resolution features $F_{lr}$ are BHT, heat flow, and thermal conductivity,
depth, northing, easting, sediment thickness, magnetic anomaly, Bouguer
gravity anomaly, K-40, U-238, and Th-232 radioactive decay measurements, and
seismic quantities. High-resolution features $F_{hr}$ are elevation, sediment
thickness, near-surface temperature, magnetic anomaly, Bouguer gravity
anomaly, and K-40, U-238, and Th-232 radioactive decay measurements.

Anchor nodes were split into $80\%$, $10\%$, and $10\%$ for training,
validation, and testing, respectively. To decorrelate spatial features and
ensure proper model evaluation, we split data using three-dimensional cubic
blocks, each with length of $10/km$. We used Feature vectors, both
$\boldsymbol{x}_{v,lr}$ and $\boldsymbol{x}_{v,hr}$, were normalized using the
training split mean $\boldsymbol{\mu}$ and standard deviation
$\boldsymbol{\sigma}$ statistics, seen in Equation 4. Toward satisfying the
three-dimensional conductive heat transfer laws, we created three graph neural
networks (GNNs) $T_{\theta}$, $Q_{\theta}$, and $K_{\theta}$ to predict
subsurface temperature, heat flow, and thermal conductivity simultaneously.
Each of these networks utilized the interpolative convolutional module where
$f_{\theta_{r}}$, $f_{\theta_{n}}$, $f_{\theta_{g}}$, and $f_{\theta_{p}}$
were each chosen as two-layer perceptrons with 128 and 64 neurons. To ensure
proper first and second gradient flow when training InterPIGNN, we used $Mish$
as an activation function. It has an unbounded positive domain which
alleviates issues of gradient saturation, bounded negative domain which
enhances generalization, and a continuously differentiable range [64].
Moreover, we incorporated two types of dropouts: (1) traditional network
neuron dropout of $10\%$ to encourage generalization across input samples
[65], and (2) graph edge dropout of $1\%$ to encourage generalization across
neighbouring nodes [66].

Whereas $Q_{\theta}$ and $K_{\theta}$ were each designed to output a quantity
each, i.e., $\hat{Q}$ and $\hat{K}$, respectively, the $T_{\theta}$
architecture involved a final node-level feedforward neural layer to output
two quantities: average surface-to-depth geothermal gradient $\hat{g}$, and
near-surface temperature correction $\Delta\hat{T}_{0}$. Consequently, the
predicted subsurface temperature at depth $z$ was computed as
$\hat{T}(z)=T_{0}+\Delta\hat{T}_{0}+\hat{g}\cdot z$, where $T_{0}$ is the
surface temperature seen in Fig. 6. We also enforced positivity using the
absolute value of all network outputs such that $\hat{T},\hat{Q},\hat{K}\geq
0$.

$\boldsymbol{\overline{x}}=\frac{\boldsymbol{x}-\boldsymbol{\mu}}{\boldsymbol{\sigma}}$
(4)

In this physics-informed supervised regression task, we incorporated multiple
loss terms. Temperature loss $L_{T}$, heat flow loss $L_{Q}$, and thermal
conductivity loss $L_{K}$ were used to learn the measured data of thermal
quantities, seen in Eqs. 5, 6, and 7, where $V_{a}^{(train)}$ represents the
training anchor nodes. Considering a finite and homogenous element (i.e.,
fixed average thermal conductivity $\overline{K}$) in three-dimensional space
at depth $z$, we combined Eqs 2 and 3 to quantify the deviation from the
three-dimensional conductive heat transfer as a loss term $L_{pde}$, seen in
Eq. 8. Furthermore, we noted that deep subsurface temperatures generally
increase along depth given intervals of hundreds of meters (except for rare
scenarios, such as subduction zones and deep aquifers undergoing convection to
shallower depths [67]), hence we introduced an additional physics-informed
loss $L_{c}$, seen in Eq. 9, to constrain our model accordingly. Note that
derivatives were computed using automatic differentiation of neural networks
[68]. Finally, the total loss $L$ was computed as the weighted sum using
hyperparameters
$\\{\lambda_{T},\lambda_{Q},\lambda_{K},\lambda_{pde},\lambda_{c}\\}$, seen in
Eq. 10.

$L_{T}=\mathbb{E}_{v\sim
V_{a}^{(train)}}\Big{[}(\hat{T}_{v}-T_{v})^{2}\Big{]}$ (5)

$L_{Q}=\mathbb{E}_{v\sim
V_{a}^{(train)}}\Big{[}(\hat{Q}_{v}-Q_{v})^{2}\Big{]}$ (6)

$L_{K}=\mathbb{E}_{v\sim
V_{a}^{(train)}}\Big{[}(\hat{K}_{v}-K_{v})^{2}\Big{]}$ (7)

$\begin{split}&L_{pde}=\mathbb{E}_{v\sim
V_{g}}\Bigg{[}\Bigg{(}\overline{K_{v}}\Big{(}\frac{\partial\hat{T}_{v}}{\partial
x}+\frac{\partial\hat{T}_{v}}{\partial y}+\frac{\partial\hat{T}_{v}}{\partial
z}\Big{)}+\hat{Q}_{v}+\\\ &\int_{0}^{z_{v}}\Bigg{(}\frac{\partial}{\partial
x}\Big{(}\hat{K}_{v}\frac{\partial\hat{T}_{v}}{\partial
x}\Big{)}+\frac{\partial}{\partial
y}\Big{(}\hat{K}_{v}\frac{\partial\hat{T}_{v}}{\partial
y}\Big{)}+\frac{\partial}{\partial
z^{\prime}}\Big{(}\hat{K}_{v}\frac{\partial\hat{T}_{v}}{\partial
z^{\prime}}\Big{)}\Bigg{)}dz^{\prime}\Bigg{)}^{2}\Bigg{]}\\\ \end{split}$ (8)

$L_{c}=\mathbb{E}_{v\sim
V_{g}}\Big{[}\max\Big{(}-\frac{\partial\hat{T}_{v}}{\partial
z},0\Big{)}^{2}\Big{]}$ (9)

$L=\lambda_{T}L_{T}+\lambda_{Q}L_{Q}+\lambda_{K}L_{K}+\lambda_{pde}L_{pde}+\lambda_{c}L_{c}$
(10)

Minimizing such multi-objective loss is an inherently challenging aspect of
training physics-informed neural networks. In this work, we alleviated this
challenge using a self-adaptive loss balancing scheme called Relative Loss
Balancing with Random Lookback (ReLoBRaLo), seen in Eq. 11 [69], where
$\alpha$ is the exponential decay rate, $\rho$ is a Bernoulli random variable
with $\mathbb{E}[\rho]\simeq 1$, $t$ is the current iteration, $t-1$ is the
previous iteration, $m$ is the number of loss terms, and $\tau$ is temperature
term. We used a uniform prior over loss weights at step zero and set
$\alpha=0.9$, $\tau=0.1$, and $\mathbb{E}[\rho]=1$ as recommended in the
original ReLoBRaLo paper [69].

$\begin{split}&\forall\ \ i\in\mathcal{L}=\\{T,Q,K,pde,c\\}\ and\
m=|\mathcal{L}|\\\
&\lambda^{bal}_{i}(t,t^{\prime})=m\cdot\frac{exp\Big{(}\frac{L_{i}(t)}{\tau
L_{i}(t^{\prime})}\Big{)}}{\sum_{j=1}^{m}exp\Big{(}\frac{L_{j}(t)}{\tau
L_{j}(t^{\prime})})}\\\
&\lambda^{hist}_{i}(t)=\rho\lambda_{i}(t-1)+(1-\rho)\lambda^{bal}_{i}(t,0)\\\
&\lambda_{i}(t)=\alpha\lambda^{hist}_{i}(t)+(1-\alpha)\lambda^{bal}_{i}(t,t-1)\end{split}$
(11)

In this study, we minimized $L$ using the Adam optimizer [70] for a total of
$400$ epochs, with initial learning rate of $0.001$, which was reduced when
the validation loss stopped to improve. Unlike uncorrelated data structures,
records (i.e., nodes) are interdependent in graph settings where the
prediction for a given node is affected by its neighbourhood. In sampling a
batch of nodes, we created a node-induced subgraph. However, especially in
large and sparse graphs, it is very likely that the neighbourhood of the
sampled nodes would be misrepresented in this subgraph. Consequently, the
InterPIGNN model would not be able to properly leverage neighbourhood
information during training, which would hinder the learning process. Our
graph is fairly large with $|V|=n=4,679,670$ nodes and $|E|=23,398,350$ edges,
so it is extremely sparse (i.e., graph density is defined as
$|E|/(|V|(|V|-1))$). To resolve this challenge, we further included the
neighbourhood of those already randomly sampled batch nodes.

###  3.4 Uncertainty Quantification

Quantifying uncertainty in subsurface engineering is crucial for decision
making and risk minimization, budgeting and investment allocation, resource
management, modeling and simulation, regulatory compliance, amongst other
considerations. Generally, there are two types of uncertainties: aleatoric or
data uncertainty, and epistemic or model uncertainty [71, 72, 73]. In this
work, we focused on estimating the epistemic uncertainty associated with the
trained model using the Monte Carlo Dropout [74]. This approach provides a
theoretical framework showing that an arbitrary neural network with dropout
applied before every weight layer, is mathematically equivalent to an
approximation to the probabilistic deep Gaussian process [75].

We can directly apply Monte Carlo Dropout using our trained model with no
changes since dropout layers are already included prior to each weight layer.
To estimate the predictive mean and predictive uncertainty we simply collected
the results of stochastic forward passes through the model. Note that these
forward passes are done concurrently, resulting in constant running time
identical to that of standard dropout. We note that the size of the variance
of the posterior distribution generated by Monte Carlo Dropout mainly depends
on (1) dropout rate, (2) the size and architecture of the neural network, and
(3) magnitude of the output quantity. Nevertheless, the posterior variance is
independent of the input data size and variance [76].

###  3.5 Prediction Explainability

Beyond achieving accurate predictions of subsurface temperature, we were also
interested in understanding which and how input quantities (i.e., features)
are most important to model predictions. Model explainability can be local or
global, where the former motivate the prediction of a given instance while the
latter provides insights into model predictions across the entire input
domain. The machine learning literature shows a variety of explainability
techniques, such as Gradients, Grad-CAM, Guided BackPropagation, DeepLIFT,
Deconvolution, amongst others [77, 78, 79, 80, 81, 82, 83]. Other approaches
provide axiomatic attribution that goes beyond using feature values and
gradients at the input vector. Integrated Gradients is an attribution
algorithm developed on axiomatic basis [79]. Particularly, this algorithm
satisfies six axioms:

  1. 1.

Completeness: given an input vector, the sum of feature attributions is equal
to the difference between the model outputs for that input vector and a
baseline vector.

  2. 2.

Implementation Invariance: if two neural networks compute identical functions
for all inputs, they receive identical attributions.

  3. 3.

Sensitivity: if an input vector differs from a baseline vector in a single
variable resulting in a distinct output then that variable should be
attributed accordingly

  4. 4.

Insensitivity: no attribution is assigned to a variable that has no impact on
the output..

  5. 5.

Linearity Preservation: the attribution of a linear combination of input
vectors is equal to the linear combination of the attributions of the
individual input vectors.

  6. 6.

Symmetry: symmetric variables with identical values receive equal
attributions.

Theoretically, the Integrated Gradients method represents the unique path
integral from a baseline vector to an input vector of interest satisfying the
aforementioned axioms. Historically, the Integrated Gradients algorithm is the
Aumann-Shapley method from cooperative game theory, which holds similar
characteristics [84]. We apply Integrated Gradients on the model prediction of
each of the $4,279,536$ grid nodes to capture local feature attributions. The
magnitudes of these attributions are then aggregated to understand the
cumulatively most influential features.

##  4 Results

We first examined the performance of our InterPIGNN model compared to the
combined SMU/NREL temperature-at-depth model, linear regression, feedforward
neural network, and EdgeConv point-based GNN. Except for the combined SMU/NREL
model, all of these models make use of the physical quantities described in
this study. Whereas these algorithms have inherently different architectures,
we ensured that they were trained and tested using the same preprocessing and
postprocessing steps. Seen in Fig. 19, we compared the different models using
the mean absolute error based on the $10\%$ hold-out test set and found that
the combined SMU/NREL mode, linear regression, feedforward, EdgeConv, and
InterPIGNN achieved $49.5$, $13.0$, $6.5$, $6.3$, and $4.8^{\circ}C$,
respectively. We first note that our models are superior to the combined
SMU/NREL model, which indicates that the considered physical quantities are
overall useful predictors of BHT. Whereas using GNN algorithms was
advantageous, the proposed InterPIGNN model was found to be the most accurate.
Comparing the InterPIGNN and EdgeConv architectures, we can theoretically
attribute this improved InterPIGNN performance to the inclusion of the target
quantities (e.g., temperature, surface heat flow, and rock thermal
conductivity) as neighbouring features in the message-passing computation.
This was further confirmed empirically through feature explainability, as
described later. Additionally, InterPIGNN simultaneously generates predictions
of surface heat flow and rock thermal conductivity which were found to have
mean absolute errors of $5.817\ mW/m^{2}$ and $0.022\ W/(C\cdot m)$,
respectively. Meanwhile, our model closely satisfied the the three-dimensional
conductive heat transfer loss with mean absolute error of $0.9\ mW/m^{2}$, and
resulted in non-decreasing temperatures with depth almost everywhere.

Figure 19: Comparison of the mean absolute error of the predicted temperature-
at-depth on the $10\%$ hold-out test set for different algorithms (smaller
values are better).

Fig. 20 shows the predicted surface heat flow across the continental United
States. Heat flow values were limited to a range of $0-300\ mW/m^{2}$ for
visual convenience. Whereas the SMU maps conveniently excluded heat flow
measurements greater than $120\ mW/m^{2}$, our model overcomes this limitation
by incorporating all heat flow measurements. The minimum, mean, and maximum
predicted surface heat flow values were $17\ mW/m^{2}$, $68\ mW/m^{2}$, and
$3,754\ mW/m^{2}$, respectively. Seen in Fig. 20, our predictions show
elevated surface heat flow across the Western United States largely due to the
presence of the Pacific Ring of Fire where tectonic movements favor volcanic
activities and geothermal heat transfer. Particularly, multiple locations
stand out such as Yellowstone in Wyoming, Basin and Range Province extending
mainly across Nevada and neighbouring states with thinner crustal thickness,
Cascade Range starting from northern California and reaching to British
Columbia, Canada and including active volcanoes like Mount St. Helens and
Mount Rainier, Coso Volcanic Field in California, Rio Grande Rift stretching
from New Mexico to Colorado, and Imperial Valley in California partly due to
the San Andreas Fault system. Away from the Western flank, elevated surface
heat flow values were also observed at the Appalachian Mountains and along the
Gulf Coast.

Figure 20: Predicted surface heat flow map.

The minimum, mean, and maximum predicted rock thermal conductivity values were
$0.099\ W/(C\cdot m)$, $2.093\ W/(C\cdot m)$, and $5.512\ W/(C\cdot m)$,
respectively. Unlike the diffusive heat flow and temperature quantities,
thermal conductivity is a rock property with compartmentalized distribution as
can be observed in the raw measurements and model predictions. It is governed
by various parameters such as rock composition, porosity, fluid saturation,
temperature, and pressure. It tends to be lower in dry and unconsolidated
rocks and varies with rock composition [85]. Seen in Fig. 21, our predictions
show a clear trend of decreasing rock thermal conductivity along depth. This
is attributed to the increase in temperature with depth which, in most rock
types, result in decreasing rock thermal conductivity [86].

Figure 21: Subsurface thermal conductivity profiles across grid nodes $V_{g}$
(gray lines) and on average (red line).

Figs. 22, 23, 24, 25, 26, 27, and 28 show the predicted temperature-at-depth
across the conterminous United States for depths of 1-7 km, respectively.
Temperatures were limited to a range of $25-250^{\circ}C$ for visual
convenience. The greatest temperature predicted at a depth of 7 km was at the
Yellowstone Caldera with a magnitude of $381^{\circ}C$. As anticipated, other
high-temperature spots were also observed around volcanic and thermally active
areas, such as Valles Caldera in New Mexico, La Garita Caldera in Colorado,
Newberry Volcano in Oregon, The Great Basin, The Geysers in California, Coso
Volcanic Field in California, Salton Buttes in California, amongst others.

Figure 22: Predicted temperature-at-depth map at a depth of 1 km. Figure 23:
Predicted temperature-at-depth map at a depth of 2 km. Figure 24: Predicted
temperature-at-depth map at a depth of 3 km. Figure 25: Predicted temperature-
at-depth map at a depth of 4 km. Figure 26: Predicted temperature-at-depth map
at a depth of 5 km. Figure 27: Predicted temperature-at-depth map at a depth
of 6 km. Figure 28: Predicted temperature-at-depth map at a depth of 7 km.

As seen in Eqs. 8-9, the InterPIGNN model was constrained to respect three-
dimensional conductive heat transfer and the fact that temperature generally
increases with depth. This resulted in physical predictions as seen in Fig. 29
which shows the predicted subsurface temperature profiles across various
locations (gray lines) as well as the average profiles (red line). We observe
the desired behavior of strictly increasing temperatures with depth with
declining geothermal gradients. Seen in Figs. 30-31, the average and local
geothermal gradient profiles can be visualized more closely across locations
(gray lines) as well as on average (red line). We generally observe decreasing
average geothermal gradients with depth, which is aligned with the established
finding that radiogenic heat generation in the Earth’s crust exponentially
declines with depth [87].

Figure 29: Subsurface temperature profiles across grid nodes $V_{g}$ (gray
lines) and on average (red line). Figure 30: Average geothermal gradient
profiles across grid nodes $V_{g}$ (gray lines) and on average (red line).
Figure 31: Local geothermal gradient profiles across grid nodes $V_{g}$ (gray
lines) and on average (red line).

With the uncertainties surrounding subsurface physical quantities, we used
Monte Carlo Dropout to estimate the model epistemic uncertainty associated
with our predictions of temperature, surface heat flow, and rock thermal
conductivity. Fig 32 shows the distribution of uncertainty across our target
thermal quantities. The median uncertainty in temperature, surface heat flow,
and rock thermal conductivity was $10.3\%$, $9.7\%$, and $5.2\%$,
respectively.

Figure 32: Temperature-at-depth uncertainty distribution as modeled using
Monte Carlo Dropout.

To verify the model performance and demonstrate its utility, we evaluated its
predictions against measured temperature log data. This involved 15 wells from
multiple geothermal projects: Utah FORGE wells 58-32, 78-32, OH1, OH4, and
Acord 1-26, Fervo Energy well Frisco-1, Cornell University Borehole
Observatory (CUBO) well, Fallon Forge 88-24, Snake River Wells WO-2 and
Camas-1, HOTSPOT Kimberly well, Paynton, Texas well, California well 5190016,
Coso Forge well 83-11, and Brady’s Hot Springs well SP-2 [88, 89, 90, 91, 92,
93]. In the case of temperature logs which were run soon after well drilling
and completion, we applied the Harrison correction [94] to account for the
disequilibrated thermal conditions in the wellbore as a result of drilling and
completion fluid circulation. It is important to note that these corrections
are correlations derived based on data measured in the field before and after
thermal equilibrium across wellbores. Thus, the log temperature data might be
a few degrees Celsius off compared to the actual rock temperatures at those
locations. We also included the NREL and SMU predictions to demonstrate how
they differ from our model. Seen in Fig. 33, the InterPIGNN model predictions
closely matched most of the temperature log data. Due to the increasing
spatial sparsity of BHT measurements along depth, our model showed reasonable
uncertainty magnitudes that were generally increasing with depth. We note that
our model was superior to the NREL model, where the latter often overestimates
subsurface temperatures. Additionally, unlike the SMU model whose predictions
tend to frequently show static temperatures along deep intervals, our model
shows increasing temperatures along depth.

Figure 33: Verification of the proposed InterPIGNN model using temperature log
data across different locations and projects. Whereas the solid line
represents the InterPIGNN model prediction, the shaded area represents the
uncertainty band. the NREL and SMU models were also visualized in dashed
lines.

##  5 Conclusion

We aggregated and processed BHT measurements for the conterminous US from
multiple sources. To construct a thermal Earth model using these records, we
introduced a novel interpolative physics-informed graph neural network module,
named InterPIGNN, that is suitable for the interpolation of point cloud data
structures. We constructed surface heat flow, and temperature and thermal
conductivity predictions for depths of $0$-$7\ km$ at an interval of $1\ km$
with spatial resolution of $18\ km^{2}$ per grid cell. Our model showed
superior temperature, surface heat flow and thermal conductivity mean absolute
errors of $4.8^{\circ}C$, $5.817\ mW/m^{2}$ and $0.022\ W/(C\cdot m)$,
respectively. We also found that elevation, crustal thickness, depth, electric
conductivity and seismic features were most important to the model
predictions. We additionally visualized local feature attributions to provide
further understanding of the model predictions across three-dimensional space.

Data Availability It is also available as feature layers on ArcGIS at
https://arcg.is/nLzzT0.

Conflict of Interest The authors have no conflicts of interest to declare that
are relevant to the content of this article.

Ethical Statement This material is the authors’ own original work, which has
not been previously published elsewhere. This paper is not currently being
considered for publication elsewhere. The paper reflects the authors’ own
research and analysis in a truthful and complete manner. The paper properly
credits the meaningful contributions of co-authors and co-researchers. The
results are appropriately placed in the context of prior and existing
research. All sources used are properly disclosed. All authors have been
personally and actively involved in substantial work leading to the paper, and
will take public responsibility for its content.

Funding Declaration No funding was received to assist with the preparation of
this manuscript.

Author Contribution Mohammad Aljubran contributed to conceptualization,
methodology, implementation and analysis, and writing first draft. Roland
Horne contributed to conceptualization, research supervision, and reviewing
and editing.

## References

  *   * Kukkonen et al. [1994] Kukkonen, I., Čermák, V., Šafanda, J.: Subsurface temperature—depth profiles, anomalies due to climatic ground surface temperature changes or groundwater flow effects. Global and Planetary Change 9(3-4), 221–232 (1994)
  * Chen et al. [2016] Chen, S., Liu, P., Liu, L., Ma, J.: Bedrock temperature as a potential method for monitoring change in crustal stress: Theory, in situ measurement, and a case history. Journal of Asian Earth Sciences 123, 22–33 (2016)
  * Vilarrasa and Rutqvist [2017] Vilarrasa, V., Rutqvist, J.: Thermal effects on geologic carbon storage. Earth-science reviews 165, 245–256 (2017)
  * Head et al. [2003] Head, I.M., Jones, D.M., Larter, S.R.: Biological activity in the deep subsurface and the origin of heavy oil. Nature 426(6964), 344–352 (2003)
  * DiPippo [1991] DiPippo, R.: Geothermal energy: Electricity production and environmental impact, a worldwide perspective. MIT PRESS, CAMBRIDGE, MA(USA)., 741–754 (1991)
  * Mock et al. [1997] Mock, J.E., Tester, J.W., Wright, P.M.: Geothermal energy from the earth: its potential impact as an environmentally sustainable resource. Annual review of Energy and the Environment 22(1), 305–356 (1997)
  * Kana et al. [2015] Kana, J.D., Djongyang, N., Raïdandi, D., Nouck, P.N., Dadjé, A.: A review of geophysical methods for geothermal exploration. Renewable and Sustainable Energy Reviews 44, 87–95 (2015)
  * Blackwell et al. [2006] Blackwell, D.D., Negraru, P.T., Richards, M.C.: Assessment of the enhanced geothermal system resource base of the united states. Natural Resources Research 15, 283–308 (2006)
  * Augustine [2016] Augustine, C.: Update to enhanced geothermal system resource potential estimate. Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States) (2016)
  * Augustine et al. [2019] Augustine, C.R., Ho, J.L., Blair, N.J.: Geovision analysis supporting task force report: Electric sector potential to penetration. Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States) (2019)
  * Tester et al. [2006] Tester, J.W., Anderson, B.J., Batchelor, A.S., Blackwell, D.D., DiPippo, R., Drake, E.M., Garnish, J., Livesay, B., Moore, M., Nichols, K., et al.: The future of geothermal energy. Massachusetts Institute of Technology 358, 1–3 (2006)
  * Smith [2019] Smith, J.: Appalachian basin temperature-depth maps and structured data in support of feasibility study of direct district heating for the cornell campus utilizing deep geothermal energy. Technical report, USDOE Geothermal Data Repository (United States); Cornell Univ., Ithaca, NY (United States) (2019)
  * Frone et al. [2015] Frone, Z., Richards, M., Blackwell, D., Augustine, C.: Shallow egs resource potential maps of the cascades. In: Proceedings, 40th Stanford Workshop on Geothermal Reservoir Engineering (2015)
  * Coolbaugh et al. [2005] Coolbaugh, M.F., Zehner, R.E., Kreemer, C.: Geothermal potential map of the great basin region, western united states. Technical report, Nevada Bureau of Mines and Geology (2005)
  * DeAngelo et al. [2023] DeAngelo, J., Burns, E.R., Gentry, E., Batir, J.F., Lindsey, C.R., Mordensky, S.P.: New maps of conductive heat flow in the great basin, usa: Separating conductive and convective influences. In: Proceedings, 48th Stanford Workshop on Geothermal Reservoir Engineering (2023)
  * Ayling et al. [2022] Ayling, B., Faulds, J., Morales Rivera, A., Koehler, R., Kreemer, C., Mlawsky, E., Coolbaugh, M., Micander, R., dePolo, C., Kraal, K., et al.: Ingenious-great basin regional dataset compilation. Technical report, USDOE Geothermal Data Repository (United States); GBCGE, NBMG, UNR (2022)
  * Batir et al. [2020] Batir, J., Richards, M., Hornbach, M., Blackwell, D., Kolker, A., Waibel, A.: Shallow geothermal potential of the snake river plain. Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States) (2020)
  * Nielson et al. [2012] Nielson, D.L., Delahunty, C., Shervais, J.W.: Geothermal systems in the snake river plain, idaho, characterized by the hotspot project. Geothermal Resources Council Transactions 36, 727–730 (2012)
  * Blackwell et al. [2011] Blackwell, D., Richards, M., Frone, Z., Batir, J., Ruzo, A., Dingwall, R., Williams, M.: Temperature-at-depth maps for the conterminous us and geothermal resource estimates. Technical report, Southern Methodist University Geothermal Laboratory, Dallas, TX (United States) (2011)
  * Mullane et al. [2016] Mullane, M., Gleason, M., McCabe, K., Mooney, M., Reber, T., Young, K.R.: An estimate of shallow, low-temperature geothermal resources of the united states. Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States) (2016)
  * Lu and Wong [2008] Lu, G.Y., Wong, D.W.: An adaptive inverse-distance weighting spatial interpolation technique. Computers & geosciences 34(9), 1044–1055 (2008)
  * Jessop [1990] Jessop, A.: Developments in solid earth geophysics. Thermal Geophysics 17 (1990)
  * Schumacher and Moeck [2020] Schumacher, S., Moeck, I.: A new method for correcting temperature log profiles in low-enthalpy plays. Geothermal Energy 8(1), 27 (2020)
  * Harrison et al. [1982] Harrison, W.E., Luza, K.V., Prater, M.L., Cheung, P.K., Ruscetta, C.: Geothermal resource assessment in oklahoma. Technical report (1982)
  * Arundel et al. [2018] Arundel, S., Bulen, A., Adkins, K., Brown, R., Lowe, A., Mantey, K., Phillips, L.: Assimilation of the national elevation dataset and launch of the 3d elevation program through the usgs spatial data infrastructure. International Journal of Cartography 4(2), 129–150 (2018)
  * Pelletier et al. [2016] Pelletier, J., Broxton, P., Hazenberg, P., Zeng, X., Troch, P., Niu, G., Williams, Z., Brunke, M., Gochis, D.: Global 1-km gridded thickness of soil, regolith, and sedimentary deposit layers. ORNL DAAC (2016)
  * Bankey et al. [2002] Bankey, V., Cuevas, A., Daniels, D., Finn, C., Hernandez, I., Hill, P., Kucks, R., Miles, W., Pilkington, M., Roberts, C., et al.: Digital data grids for the magnetic anomaly map of North America: USGS Open-File Report 02-414 (2002)
  * Kucks [1999] Kucks, R.P.: Bouguer gravity anomaly data grid for the conterminous us (1999)
  * Duval et al. [2005] Duval, J.S., Carson, J.M., Holman, P.B., Darnley, A.G.: Terrestrial radioactivity and gamma-ray exposure in the united states and canada. US geological survey open-file report 1413, 2005 (2005)
  * Hill et al. [2009] Hill, P.L., Kucks, R.P., Ravat, D.: Aeromagnetic and aeroradiometric data for the conterminous united states and alaska from the national uranium resource evaluation (nure) program of the us department of energy. Technical report, US Geological Survey (2009)
  * Poletto et al. [2018] Poletto, F., Farina, B., Carcione, J.M.: Sensitivity of seismic properties to temperature variations in a geothermal reservoir. Geothermics 76, 149–163 (2018)
  * Yoshino [2011] Yoshino, T.: Electrical properties of rocks. Encyclopedia of Earth Sciences Series, 270–276 (2011)
  * Shinevar et al. [2023] Shinevar, W.J., Golos, E.M., Jagoutz, O., Behn, M.D., Van Der Hilst, R.D.: Mantle thermochemical variations beneath the continental united states through petrologic interpretation of seismic tomography. Earth and Planetary Science Letters 602, 117965 (2023)
  * Schmandt et al. [2015] Schmandt, B., Lin, F.-C., Karlstrom, K.E.: Distinct crustal isostasy trends east and west of the rocky mountain front. Geophysical Research Letters 42(23), 10–290 (2015)
  * Shen and Ritzwoller [2016] Shen, W., Ritzwoller, M.H.: Crustal and uppermost mantle structure beneath the united states. Journal of Geophysical Research: Solid Earth 121(6), 4306–4342 (2016)
  * Golos et al. [2020] Golos, E.M., Fang, H., Hilst, R.D.: Variations in seismic wave speed and vp/vs ratio in the north american lithosphere. Journal of Geophysical Research: Solid Earth 125(12), 2020–020574 (2020)
  * Murphy et al. [2023] Murphy, B.S., Bedrosian, P.A., Kelbert, A.: Geoelectric constraints on the precambrian assembly and architecture of southern laurentia (2023)
  * Rukundo and Cao [2012] Rukundo, O., Cao, H.: Nearest neighbor value interpolation. arXiv preprint arXiv:1211.1768 (2012)
  * Ryu et al. [2002] Ryu, J.-S., Kim, M.-S., Cha, K.-J., Lee, T.H., Choi, D.-H.: Kriging interpolation methods in geostatistics and dace model. KSME International Journal 16, 619–632 (2002)
  * Bhatia and Arora [2016] Bhatia, G.S., Arora, G.: Radial basis function methods for solving partial differential equations-a review. Indian Journal of Science and Technology 9(45), 1–18 (2016)
  * Mo et al. [2022] Mo, J., Shou, H., Chen, W.: Implicit surface reconstruction via rbf interpolation: A review. Recent Patents on Engineering 16(5), 49–66 (2022)
  * Hutchinson [1995] Hutchinson, M.F.: Interpolating mean rainfall using thin plate smoothing splines. International journal of geographical information systems 9(4), 385–403 (1995)
  * Wahba [1990] Wahba, G.: Spline models for observational data. Technical report, SIAM (1990)
  * Qiao et al. [2019] Qiao, P., Li, P., Cheng, Y., Wei, W., Yang, S., Lei, M., Chen, T.: Comparison of common spatial interpolation methods for analyzing pollutant spatial distributions at contaminated sites. Environmental geochemistry and health 41, 2709–2730 (2019)
  * Li and Heap [2014] Li, J., Heap, A.D.: Spatial interpolation methods applied in the environmental sciences: A review. Environmental Modelling & Software 53, 173–189 (2014)
  * Qi et al. [2017a] Qi, C.R., Su, H., Mo, K., Guibas, L.J.: Pointnet: Deep learning on point sets for 3d classification and segmentation. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 652–660 (2017)
  * Qi et al. [2017b] Qi, C.R., Yi, L., Su, H., Guibas, L.J.: Pointnet++: Deep hierarchical feature learning on point sets in a metric space. Advances in neural information processing systems 30 (2017)
  * Zhao et al. [2019] Zhao, H., Jiang, L., Fu, C.-W., Jia, J.: Pointweb: Enhancing local neighborhood features for point cloud processing. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 5565–5573 (2019)
  * Liu et al. [2019a] Liu, Y., Fan, B., Xiang, S., Pan, C.: Relation-shape convolutional neural network for point cloud analysis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 8895–8904 (2019)
  * Liu et al. [2019b] Liu, Y., Fan, B., Meng, G., Lu, J., Xiang, S., Pan, C.: Densepoint: Learning densely contextual representation for efficient point cloud processing. In: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 5239–5248 (2019)
  * Hua et al. [2018] Hua, B.-S., Tran, M.-K., Yeung, S.-K.: Pointwise convolutional neural networks. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 984–993 (2018)
  * Besl and McKay [1992] Besl, P.J., McKay, N.D.: Method for registration of 3-d shapes. In: Sensor Fusion IV: Control Paradigms and Data Structures, vol. 1611, pp. 586–606 (1992). Spie
  * Deng et al. [2018a] Deng, H., Birdal, T., Ilic, S.: Ppfnet: Global context aware local features for robust 3d point matching. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 195–205 (2018)
  * Deng et al. [2018b] Deng, H., Birdal, T., Ilic, S.: Ppf-foldnet: Unsupervised learning of rotation invariant 3d local descriptors. In: Proceedings of the European Conference on Computer Vision (ECCV), pp. 602–618 (2018)
  * Gojcic et al. [2019] Gojcic, Z., Zhou, C., Wegner, J.D., Wieser, A.: The perfect match: 3d point cloud matching with smoothed densities. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 5545–5554 (2019)
  * Choy et al. [2019] Choy, C., Park, J., Koltun, V.: Fully convolutional geometric features. In: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 8958–8966 (2019)
  * Aoki et al. [2019] Aoki, Y., Goforth, H., Srivatsan, R.A., Lucey, S.: Pointnetlk: Robust & efficient point cloud registration using pointnet. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 7163–7172 (2019)
  * Simonovsky and Komodakis [2017] Simonovsky, M., Komodakis, N.: Dynamic edge-conditioned filters in convolutional neural networks on graphs. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 3693–3702 (2017)
  * Wang et al. [2019] Wang, Y., Sun, Y., Liu, Z., Sarma, S.E., Bronstein, M.M., Solomon, J.M.: Dynamic graph cnn for learning on point clouds. ACM Transactions on Graphics (tog) 38(5), 1–12 (2019)
  * Lin et al. [2020] Lin, Z.-H., Huang, S.-Y., Wang, Y.-C.F.: Convolution in the cloud: Learning deformable kernels in 3d graph convolution networks for point cloud analysis. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 1800–1809 (2020)
  * Jaupart et al. [2007] Jaupart, C., Labrosse, S., Lucazeau, F., Mareschal, J.: 7.06-temperatures, heat and energy in the mantle of the earth. Treatise on geophysics 7, 223–270 (2007)
  * Cengel et al. [2012] Cengel, Y., Cimbala, J., Turner, R.: EBOOK: Fundamentals of Thermal-Fluid Sciences (SI Units). McGraw Hill, Singapore (2012)
  * Mareschal and Jaupart [2013] Mareschal, J.-C., Jaupart, C.: Radiogenic heat production, thermal regime and evolution of continental crust. Tectonophysics 609, 524–534 (2013)
  * Misra [2019] Misra, D.: Mish: A self regularized non-monotonic activation function. arXiv preprint arXiv:1908.08681 (2019)
  * Srivastava et al. [2014] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: Dropout: a simple way to prevent neural networks from overfitting. The journal of machine learning research 15(1), 1929–1958 (2014)
  * Papp et al. [2021] Papp, P.A., Martinkus, K., Faber, L., Wattenhofer, R.: Dropgnn: Random dropouts increase the expressiveness of graph neural networks. Advances in Neural Information Processing Systems 34, 21997–22009 (2021)
  * Ziagos and Blackwell [1986] Ziagos, J., Blackwell, D.: A model for the transient temperature effects of horizontal fluid flow in geothermal systems. Journal of Volcanology and Geothermal Research 27(3-4), 371–397 (1986)
  * Baydin et al. [2018] Baydin, A.G., Pearlmutter, B.A., Radul, A.A., Siskind, J.M.: Automatic differentiation in machine learning: a survey. Journal of machine learning research 18(153), 1–43 (2018)
  * Bischof and Kraus [2021] Bischof, R., Kraus, M.: Multi-objective loss balancing for physics-informed deep learning. arXiv preprint arXiv:2110.09813 (2021)
  * Kingma and Ba [2014] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * Gal et al. [2016] Gal, Y., et al.: Uncertainty in deep learning (2016)
  * Hüllermeier and Waegeman [2021] Hüllermeier, E., Waegeman, W.: Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods. Machine Learning 110, 457–506 (2021)
  * Kendall and Gal [2017] Kendall, A., Gal, Y.: What uncertainties do we need in bayesian deep learning for computer vision? Advances in neural information processing systems 30 (2017)
  * Gal and Ghahramani [2016] Gal, Y., Ghahramani, Z.: Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In: International Conference on Machine Learning, pp. 1050–1059 (2016). PMLR
  * Damianou and Lawrence [2013] Damianou, A., Lawrence, N.D.: Deep gaussian processes. In: Artificial Intelligence and Statistics, pp. 207–215 (2013). PMLR
  * Verdoja and Kyrki [2020] Verdoja, F., Kyrki, V.: Notes on the behavior of mc dropout. arXiv preprint arXiv:2008.02627 (2020)
  * Linardatos et al. [2020] Linardatos, P., Papastefanopoulos, V., Kotsiantis, S.: Explainable ai: A review of machine learning interpretability methods. Entropy 23(1), 18 (2020)
  * Simonyan et al. [2013] Simonyan, K., Vedaldi, A., Zisserman, A.: Deep inside convolutional networks: Visualising image classification models and saliency maps. arXiv preprint arXiv:1312.6034 (2013)
  * Sundararajan et al. [2017] Sundararajan, M., Taly, A., Yan, Q.: Axiomatic attribution for deep networks. In: International Conference on Machine Learning, pp. 3319–3328 (2017). PMLR
  * Shrikumar et al. [2017] Shrikumar, A., Greenside, P., Kundaje, A.: Learning important features through propagating activation differences. In: International Conference on Machine Learning, pp. 3145–3153 (2017). PMLR
  * Springenberg et al. [2014] Springenberg, J.T., Dosovitskiy, A., Brox, T., Riedmiller, M.: Striving for simplicity: The all convolutional net. arXiv preprint arXiv:1412.6806 (2014)
  * Zeiler and Fergus [2014] Zeiler, M., Fergus, R.: European conference on computer vision. In: European Conference on Computer Vision., Berlin/Heidelberg, Germany: Springer, pp. 818–833 (2014)
  * Selvaraju et al. [2017] Selvaraju, R.R., Cogswell, M., Das, A., Vedantam, R., Parikh, D., Batra, D.: Grad-cam: Visual explanations from deep networks via gradient-based localization. In: Proceedings of the IEEE International Conference on Computer Vision, pp. 618–626 (2017)
  * Friedman [2004] Friedman, E.J.: Paths and consistency in additive cost sharing. International Journal of Game Theory 32, 501–518 (2004)
  * Song et al. [2023] Song, J., Wang, G., Xing, L., Qian, J., Dai, L., Di, H.: Influencing factors of rock thermal conductivity and applicability evaluation of its mixing law predictive models. Geothermics 110, 102680 (2023)
  * Chen et al. [2021] Chen, C., Zhu, C., Zhang, B., Tang, B., Li, K., Li, W., Fu, X.: Effect of temperature on the thermal conductivity of rocks and its implication for in situ correction. Geofluids 2021, 1–12 (2021)
  * Brady et al. [2006] Brady, R.J., Ducea, M.N., Kidder, S.B., Saleeby, J.B.: The distribution of radiogenic heat production as a function of depth in the sierra nevada batholith, california. Lithos 86(3-4), 229–244 (2006)
  * Allis et al. [2018] Allis, R., Gwynn, M., Hardwick, C., Hurlbut, W., Moore, J.: Thermal characteristics of the forge site, milford, utah. Geothermal Resources Council Transactions 42, 15 (2018)
  * Cladouhos et al. [2015] Cladouhos, T.T., Petty, S., Swyer, M.W., Nordin, Y., Garrison, G., Uddenberg, M., Grasso, K., Stern, P., Sonnenthal, E., Foulger, G., et al.: Newberry egs demonstration: Phase 2.2 report. Technical report, AltaRock Energy, Seattle, WA (United States) (2015)
  * Purwamaska and Fulton [2023] Purwamaska, I., Fulton, P.M.: Preliminary constraints on thermal conditions within the cornell university borehole observatory (cubo) (2023)
  * Podgorney [1991] Podgorney, R.: Snake river plain forge: Well data for wo-2. Technical report, USDOE Geothermal Data Repository (United States); Idaho National Lab.(INL … (1991)
  * Blankenship [2017] Blankenship, D.A.: The proposed fallon forge site: Phase 2 update. Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States) (2017)
  * Shervais et al. [2013] Shervais, J.W., Schmitt, D.R., Nielson, D., Evans, J.P., Christiansen, E.H., Morgan, L., Pat Shanks, W., Prokopenko, A.A., Lachmar, T., Liberty, L.M., et al.: First results from hotspot: the snake river plain scientific drilling project, idaho, usa. Scientific Drilling 15, 36–45 (2013)
  * Harrison [1983] Harrison, W.E.: Geothermal resource assessment in oklahoma: Modification of a report prepared for us department of energy, division of geothermal energy, under contract no. de-as07-80id12172, october 1981. Technical report, Oklahoma Geological Survey (1983)