# Deep learning based surrogate modeling for thermal plume prediction of
groundwater heat pumps

Kyle Davis1,∗, Raphael Leiteritz2,∗, Dirk Pflüger2, Miriam Schulte1 1\.
Simulation of Large Systems, IPVS, University of Stuttgart, 70569, Stuttgart,
Germany  2\. Scientific Computing, IPVS, University of Stuttgart, 70569,
Stuttgart, Germany  * Authors contributed equally<EMAIL_ADDRESS>stuttgart.de<EMAIL_ADDRESS>
###### Abstract

The ability for groundwater heat pumps to meet space heating and cooling
demands without relying on fossil fuels, has prompted their mass roll-out in
dense urban environments. In regions with high subsurface groundwater flow
rates, the thermal plume generated from a heat pump’s injection well can
propagate downstream, affecting surrounding users and reducing their heat pump
efficiency. To reduce the probability of interference, regulators often rely
on simple analytical models or high-fidelity groundwater simulations to
determine the impact that a heat pump has on the subsurface aquifer and
surrounding heat pumps. These are either too inaccurate or too computationally
expensive for everyday use. In this work, a surrogate model was developed to
provide a quick, high accuracy prediction tool of the thermal plume generated
by a heat pump within heterogeneous subsurface aquifers. Three variations of a
convolutional neural network were developed that accepts the known groundwater
Darcy velocities as discrete 2D inputs and predicts the temperature within the
subsurface aquifer around the heat pump. A data set consisting of 800
numerical simulation samples, generated from random permeability fields and
pressure boundary conditions, was used to provide pseudo-randomized Darcy
velocity fields as input fields and the temperature field solution for
training the network. The subsurface temperature field output from the network
provides a more realistic temperature field that follows the Darcy velocity
streamlines, while being orders of magnitude faster than conventional high-
fidelity solvers.

  * February 2023

##  1 Introduction

Addressing the challenges of climate change and growing energy costs are
forcing cities to increase their usage of renewable energies. The European
Union 2030 climate and energy framework requires a minimum reduction of 40% in
greenhouse gas emissions, with a minimum of 32% renewable energy usage [1].
Energy usage can be significantly reduced by reducing the space heating and
cooling demands of buildings, which have been met predominantly with fossil
fuel based resources. An increasingly popular alternative is to use shallow
groundwater heat pumps, which have proved to be an effective alternative
within cities.

Rolling out groundwater heat pumps on a city-wide scale is not without its
challenges. City regulators must decide how many heat pumps can be effectively
installed, how much can they be used and where can they be placed. The mass
roll-out of heat pumps within a confined urban environment often leads to
negative interference, where one heat pump changes the temperature of the
water within a shared subsurface aquifer, reducing the efficiency of
downstream heat pump systems. Additionally, the synergistic use of heat pump
systems are often over-looked due to the added complexity of determining their
mutual interactions. Therefore, monitoring any negative interactions and
improving the overall efficiency by optimizing the usage and location of the
heat pumps is imperative.

Numerous approaches to optimize heat pump usage and layout exist. Most methods
involve either performing expensive numerical groundwater simulations using
classical physics-based solver [2, 3, 4], or by using a cheaper approximation
of the heat transfer within the subsurface [5, 6]. The classical numerical
simulations offer highly accurate solutions for the subsurface thermal field,
especially in the case of heterogeneous groundwater properties but at the
expense of computational cost and increased run-time. The simulations often
require extensive model preparation, calibration and validation, before
running on large computing clusters. Only a few variations of the layout of
heat pumps can be simulated within a reasonable time frame, limiting the
practicality of this method for large cities with potentially thousands of
heat pumps. Alternatively, using simplified approximations of the subsurface
thermal profile allows for almost real-time solutions of the subsurface
thermal profile induced by the heat pump, but at the expense of accuracy and
only is valid for certain groundwater properties.

A surrogate model that is able to reproduce the classical numerical simulation
results, while being orders of magnitude cheaper, provides an attractive
solution to the problem. This can be used to create a virtual model of the
subsurface domain, allowing city and energy planners to perform optimization
studies that are otherwise computationally infeasible, or provide an online
monitoring tool of the subsurface of the city while having fast and accurate
results. Many virtual models, or digital twins, of cities typically focus only
on the above surface world and rarely model what happens underneath.

In order to provide the level of accuracy as the high-fidelity solver, the
surrogate model must be able to account for the subsurface heterogeneity.
Popular methods for creating a reduced order model (ROM) are projection-based
methods, such as proper orthogonal decomposition [7, 8, 9]. Projection methods
use snapshots of input-output pair information, where previous high-fidelity
simulation runs are used to obtain these pairs. Once a ROM is generated,
access to the underlying numerical model is typically required to build the
surrogate model. This is sometimes problematic if the underlying numerical
model is inaccessible, or if a truly black-box surrogate model is required.

Recently, deep learning for black-box surrogate modeling of physical problems
have been created using artificial neural networks (ANN) [10, 11], requiring
access to the input-output data pairs only and not the underlying numerical
model or internal simulation solver. ANNs have already been applied to
modeling reacting flows [12, 13], predicting airflow over airfoils [14] and
modeling fluid flow within porous media [15, 16]. In addition to surrogate
modeling, ANNs can provide coarse to fine mesh mapping, where the network
learns to approximate the fine scale information with access to the coarse
scale features only [17], providing a middle-ground feature between a ROM and
high-fidelity solution. Critical to the success of ANNs are the data, where
generating sufficient simulation data to train an ANN is often a stumbling
block for physics-based surrogate models. Physics-informed neural networks
have developed due to the underlying difficulty of obtaining expensive
numerical simulation data [18]. Instead of using the input-output data pairs,
the network is instead trained knowing that the output must satisfy a partial
differential equation (PDE) that defines the physical problem. Therefore, only
input data is required and the network loss function consists of the mismatch
between the output result and the PDE that it must satisfy. This has been
applied to numerous physics-based surrogate models where the governing
equations are known [19, 20, 21].

The objective of our study was to utilize artificial neural networks to create
a fast surrogate model that solves for the local subsurface thermal field due
to the presence of a groundwater heat pump, while accounting for the
subsurface heterogeneity of the permeability and Darcy velocity fields. The
surrogate model can be used to perform heat pump layout and usage
optimization, which may consist of thousands of layout configurations, or
serve as a fast approximation to build into an online evaluation tool for
groundwater heat pump management.

##  2 Method

###  2.1 Open-loop groundwater heat pumps

Shallow groundwater heat pumps are devices that transfer heat to and from a
shallow aquifer beneath the surface in order to provide space heating or
cooling. Open-loop groundwater heat pumps, depicted in figure 1, function by
extracting water from the subsurface aquifer, passes the water through a heat
exchanger, before re-injecting the water back into the subsurface. When
operating in heating mode, i.e., heating a building, the energy in the water
is passed into the building through the heat exchanger, thereby cooling the
water and re-injecting the fluid at a lower temperature back into the aquifer.
Conversely, in cooling mode the energy is passed from the building and into
the fluid, warming up the water and re-injecting the warmer fluid back into
the aquifer. The natural flow of the groundwater within the subsurface aquifer
causes the re-injected water, now at an elevated or reduced temperature to the
rest of the aquifer, to be pulled along with the movement of the groundwater,
causing a thermal plume to form downstream of the injection well. If this
plume reaches the extraction well of a nearby downstream heat pump, the
heightened or lowered temperature of the groundwater within the plume may
reduce the efficiency of the downstream system to the extent that it becomes
unusable. When a small urban area is densely populated with heat pumps, the
overall efficiency may be reduced more than if only a handful of well-placed
heat pumps were used.

Figure 1: Open loop groundwater heat pump, where the flow of groundwater
causes a plume to develop that stretches downstream of the injection well.
Image modified from [22]

###  2.2 Subsurface modelling

Before any heat pump layout optimization study or heat pump monitoring can be
performed, a suitable model for the groundwater temperature is required. There
are two popular methods that exist to model the thermal plume that develops
from an open-loop heat pump. The linear advective heat transport model (LAHM)
from [23] is a common analytical formulation that defines the change in
groundwater temperature around the heat pump over time compared to the far-
away background temperature

$\displaystyle\Delta T(x,y,t)=\frac{Q\cdot\Delta T_{inj}}{4\cdot n_{e}\cdot
M\cdot v_{a}\cdot\sqrt{\pi\cdot\alpha_{T}}}\cdot
exp\left(\frac{x-r}{2\cdot\alpha_{L}}\right)\cdot$
$\displaystyle\frac{1}{\sqrt{r}}\cdot erfc\left(\frac{r-v_{a}\cdot
t/R}{2\cdot\sqrt{v_{a}\cdot\alpha_{L}\cdot t/R}}\right),$ (1)

where $\Delta T(x,y,t)$ is the time-dependent temperature difference at
coordinates $x$ and $y$ (the heat pump is at coordinates $x$ = 0 and $y$ = 0)
at time $t$, $Q$ is the injection mass flow rate, $\Delta T_{inj}$ is the
difference between the background temperature and the injection temperature,
$v_{a}$ is the groundwater velocity at the heat pump injection well, $r$ is
the radial distance from the injection well, $M$ is the aquifer thickness, and
$\alpha_{L}$ and $\alpha_{T}$ are the longitudinal and tangential dispersivity
values. The disadvantage of the LAHM is that the plume is dependent only on
the velocity magnitude in subsurface aquifer at the heat pump location with
constant subsurface properties. Therefore, no heterogeneity of the subsurface
is accommodated. The uni-directional thermal plumes determined using the LAHM
are shown in figure 2 as examples of the solution.

The LAHM has already been integrated into an online evaluation tool shown in
figure 2 (a), that determines the iso-lines of the change in subsurface
temperature $\Delta T(x,y,t)$ due to the presence of a heat pump. This
provides a fast surrogate model, but does not model the heterogeneity as shown
in figure 2 (b), where the LAHM solution iso-lines are overlaid onto the
classical physics numerical simulation results using a 3D groundwater
mechanics solver. The classical solver accounts for subsurface heterogeneity,
causing non-uniform Darcy velocity streamlines and varying direction of the
thermal plume.

(a) Online evaluation tool

(b) Analytical vs. numerical solution

Figure 2: Temperature iso-lines of the LAHM analytical solution applied to an
online groundwater heat pump planning and optimization tool (left) and a
comparison between the LAHM analytical solution temperature iso-lines overlaid
onto the 2D numerical solution (right).

More complex 3D numerical groundwater models can provide a more accurate
solution for the groundwater temperature field. The governing equations for
the subsurface groundwater model is derived from [24]. The subsurface model
applies to a single phase and variably saturated fluid, defined on a single
space-time domain. The subsurface fluid flow is governed first by the
conservation of mass,

$\frac{\delta}{\delta t}\left(\varphi
s\eta\right)+\nabla\cdot\left(\eta\mathbf{q}\right)=Q_{w},$ (2)

with the porosity $\varphi$ [-], saturation ratio $s$ $[m^{3}m^{-3}]$, molar
density $\eta$ $[kmolm^{-3}]$, the Darcy velocity $\mathbf{q}$ and mass
source/sink term $Q_{w}$ $[kmol\cdot m^{-3}\cdot s^{-1}]$. The groundwater
temperature is modeled by the inclusion of the conservation of energy,

$\displaystyle\frac{\delta}{\delta t}\left(\varphi s\eta
U+\left(1-\varphi\right)\rho_{p}c_{p}T\right)$
$\displaystyle+\nabla\cdot\left(\eta\mathbf{q}H-\kappa\nabla T\right)=Q_{e},$
(3)

with the rock density $\rho_{r}$, heat capacity $c_{p}$ and thermal
conductivity $\kappa$ of the porous medium - fluid mixture and the energy
source/sink term $Q_{e}$. The fluid enthalpy $H$ is related to the internal
energy of the water through the expression

$U=H-\frac{P}{\eta}.$ (4)

The Darcy velocity $\mathbf{q}=(q_{x},q_{y},q_{z})^{T}$ in $[m\cdot s^{-1}]$
is defined as

$\mathbf{q}=-\frac{\mathbf{K}(s)}{\mu}\nabla\left(P-\rho gz\right),$ (5)

with the relative permeability field $\mathbf{K}(s)$ $[m^{2}]$, the viscosity
$\mu$ $[Pa\cdot s]$, subsurface water pressure $P$ $[Pa]$, gravitational
constant $g$ $[m\cdot s^{-2}]$ and the relative reference height $z$ $[m]$.
The constant pressure initial condition for the subsurface model is defined on
the boundary of the domain $\delta\Omega_{P}$, inducing the Darcy velocity
$\mathbf{q}=\mathbf{K}\nabla P$ throughout $\Omega_{P}$. The boundary
condition of a heat pump can be specified by either injecting energy into the
domain via the heat flux $Q_{e}$, or by injecting a mass of water $Q_{w}$ at a
predetermined temperature. The latter method more resembles the operation of
an open-loop heat pump, which pumps water back into the subsurface at a new
temperature that is dependent on the extraction well temperature.

The operation of commercial heat pumps utilize a constant temperature
difference across the heat exchanger, and alters the amount of energy passed
to/from the groundwater by increasing or decreasing the mass flow rate
extracted from the subsurface aquifer. The difference between the injection
well $T_{inj}$ and extraction well $T_{ext}$ temperature is defined as $\Delta
T=T_{inj}-T_{ext}=5$ when operation in cooling mode and $\Delta T=-5$ when
operating in heating mode. The extraction and injection mass flow rate is
defined from the amount of energy passed into-or-from the fluid and $\Delta T$
by

$\dot{m}=\frac{\dot{Q}}{c_{p}\Delta T}\,,$ (6)

with the energy transferred to/from the fluid $\dot{Q}$ $[J]$, the specific
heat of water (assumed constant) $c_{p}$ = 4184 $[J\cdot kg^{-1}\cdot K^{-1}]$
and the temperature difference $\Delta T$.

###  2.3 Neural network design

####  2.3.1 Input data:

The neural network operation is divided into two parts: offline and online
stages. The offline stage involves generating and processing the training
data, training the network and finally testing the network accuracy. The
online stage involves feeding the input to the network and obtaining the
output during normal usage of the network. The offline step is typically
computationally expensive, but is tolerable as this is only performed once in
order to generate a network. The online step is intended to be performed
thousands if not millions of times, thereby requiring a fast prediction
performance in the online step. To minimize the computational cost of the
network prediction step, the surrogate model must not rely on performing any
numerical simulation during the online step and therefore, all input data must
be readily available. A baseline simulation of the entire subsurface of Munich
already exists and contains the subsurface permeability, pressure and Darcy
velocity throughout the city’s subsurface. These data can simply be extracted
from the baseline simulation result and input into the network during the
online phase without having to perform any expensive numerical simulation.

For this study of building a surrogate model, we ignore the time-dependent
behavior of the thermal plume propagation and only consider the steady state
solution. Examining equation 2.2 and ignoring the time-derivative, the thermal
plume is dependent on an advection term and a diffusion term. We assume that
the plume behavior is dominated by the advection term of this equation, where
figure 2 strengthens this assumption as the thermal plume tends to follow the
velocity streamlines. As the subsurface Darcy velocity is readily available
and likely to provide meaningful information to the network in order to
predict the plume behavior, we use the Darcy velocity components $q_{x}$ and
$q_{y}$ as inputs into the network. Previous work on building a surrogate
model of the subsurface flow accepted the groundwater permeability as input
and predicted the groundwater pressure and Darcy velocities [19]. However,
this is only applicable for scenarios with a constant pressure boundary
condition, and therefore not applicable for general use such as with our
surrogate model.

####  2.3.2 Output data:

Obtaining the output results from the surrogate model must be orders of
magnitude faster to attain than the high-fidelity model in order to be useful,
while remaining accurate. For our specific application of optimization and
condition monitoring, only the temperature output is required for build a
practical surrogate model. The network input data (obtained from the baseline
simulation data), does not include the influence that the heat pump has on the
Darcy velocities or thermal field. The baseline simulations do not have the
required output data necessary to train the network using supervised learning,
which must be sourced elsewhere. Therefore, smaller 2D numerical simulations
are performed in order to generate input-output data pairs for the network
training and testing. However, these smaller numerical simulations are much
smaller than the baseline simulation case, such that the data generation time
is feasible for our model.

####  2.3.3 Network architecture:

The available input and output data structure must be considered when
designing the network architecture. The output of the 2D numerical groundwater
simulation, which will be used as input and output of the network, can be
treated as an image where the value at each pixel defines the Darcy velocity
(input) or temperature (output) at the cell centers of the finite volume mesh.
The input data are provided as a two-channel image
$\varphi_{in}\in\mathbb{R}^{2\times 65\times 65}$ of the Darcy velocity
$q_{x}$ and $q_{y}$ of size 65$\times$65 pixels each ($\varphi_{in,1}$ and
$\varphi_{in,2}$). The network output $\varphi_{out}\in\mathbb{R}^{65\times
65}$ is a single channel image of 65$\times$65 of the temperature field. The
location of the heat pump is always at the center pixel. An image size of
65$\times$65 was selected so that an equal number of pixels surround the
center pixel in each direction. The image-like nature of the input and output
data naturally allows the use of a convolutional neural network (CNN), which
is favorable for image like data. Previous work on predicting the behavior of
physical systems have shown that CNNs are favorable using image-like data [14,
19, 17, 25].

Within this work, we built and evaluated three variations of the TurbNet
architecture by Thuerney et. al. [14], which is also a variant of the U-Net
architecture by Ronneberger et. al. [26]. The TurbNet architecture features
skip connections between the encoding and decoding steps.

Each layer of the CNN applies a 2D convolution operation to the input image by
sweeping a fixed size kernel (e.g. 4x4) over the input and multiplying it with
the underlying image data. The kernel weights $\mathbf{W}_{c}$, are parameters
which are optimized during the network’s training procedure. By having these
weights shared for the whole image, the convolutional layer essentially learns
to pick up translation invariant features. A user-defined number of features,
each having its own kernel, can then be learned to create multiple outputs per
layer. Usually, after each convolutional layer an aggregation operation is
applied such as a pooling operation where multiple pixels are combined to
create a coarser output. Recursively continuing this process allows the
network to pick up different features at multiple scales.

For the U-Net architecture this process is reversed after reaching a pre-
defined bottleneck state where the coarsest features are represented. This
bottleneck is then followed by inverse convolutions to reconstruct higher-
resolution images for the output path. This is comparable to simple auto-
encoder architectures which also feature an encoder-bottleneck-decoder
structure.

The entire encoding-decoding procedure can be characterized as a mapping $S$,
from the input images and kernel weights to the prediction output
$\mathbf{\varphi}_{pred}$,

$S\left(\mathbf{W}_{c},\mathbf{\varphi}_{in,1},\mathbf{\varphi}_{in,2}\right)\mapsto\mathbf{\varphi}_{pred}.$
(7)

The network is trained to minimize the difference between the known output
$\mathbf{\varphi}_{out}$ and the network prediction output
$\mathbf{\varphi}_{pred}$ by varying the kernel weights $\mathbf{W}_{c}$ and
other network parameters $\theta$

$\mathop{\mathrm{arg\,min}}_{W_{c},\theta}\sum_{i}^{N_{data}}\|\mathbf{\varphi}_{pred}-\mathbf{\varphi}_{out}\|.$
(8)

Both the encoding and decoding steps are divided into multiple layers, which
each contain a rectified linear unit (ReLU) activation function with a max
pooling layer of 2$\times$2, a stride of $2$ for down-sampling the image size
and batch normalization. The image down-sampling halves the image pixel size
per direction (the first convolution layer only reduces the image size from
65$\times$65 to 64$\times$64). The total number of trainable parameters of the
network can be increased by increasing the number of layers or increasing the
number of initial features. This also affects the accuracy of the network.

The three architectures tested are:

  1. 1.

TurbNet-Geo: 6 layers with skip connections

  2. 2.

TurbNet-Geo-Light: 4 layers with skip connections.

  3. 3.

TurbNet-Geo-NoSkip-Light: 4 layers without skip connections (TurbNet-Geo-Light
without the connections).

With the latter two architectures, we test the model performance with regard
to fewer network parameters and the elimination of skip connections. The skip
connections add a concatenation step between the encoding and decoding layers
in order to retain features in the encoding steps. The TurbNet-Geo-Light
network architecture is illustrated in figure 3. The 2 channel input is
provided on the left, with the network output on the right.

Figure 3:  TurbNet-Geo architecture containing skip connections for the
temperature plume prediction. The network accepts a 2-channel input of the $x$
and $y$-direction Darcy velocity magnitudes and outputs a single channel of
the temperature field. The bottleneck contains 32$\times$ the amount of
initial features.

###  2.4 Data generation

In order to study whether the CNN is able to predict the thermal plume with
sufficient accuracy, a simple 2D groundwater model is utilized. This is
sufficient for two reasons: firstly, the network only accepts 2D images as
inputs and secondly, the thermal plume only influences the local area around a
heat pump. This allows for relatively cheap numerical simulations to be
performed to obtain the Darcy velocities without the influence of the heat
pump, followed by re-running the same model but with the heat pump activated
to obtain the thermal plume profile. As the thermal plume profiles are not
available from the current baseline simulation, it would be too costly to run
multiple large, high-fidelity simulations of a large area of the city, only to
obtain local temperature profiles around the heat pumps. Therefore, the data
generation procedure is divided into 4 steps:

  1. 1.

Input field generation - generates the numerical simulation input files with
randomized permeability fields and pressure boundary conditions,

  2. 2.

Without-heat-pump evaluation - generates the Darcy velocity fields without the
influence of a heat pump,

  3. 3.

With-heat-pump evaluation - the ”Without-heat-pump” simulation is repeated but
with the heat pump activated,

  4. 4.

Data manipulation - pre-processing of the input and output data for the
network.

All input and output training and testing data of the numerical groundwater
simulations were generated using the subsurface simulation software PFLOTRAN
v3 [27]. The data generation simulations were performed on a workstation with
an AMD Ryzen Threadripper 3960X 24-core Processor at the University of
Stuttgart, running Ubuntu 20.04.

####  2.4.1 Input field generation:

The prediction quality of the network is dependent on the quality of the
training data. Therefore, the purpose of the input field generation procedure
is to generate a variety of Darcy velocities and flow directions. To generate
the variety of velocity data, randomized permeability fields are combined with
randomized pressure gradient boundary conditions to generate seemingly
arbitrary Darcy velocities that are within a suitable velocity range. For each
data sample (one permeability field and one pressure gradient boundary
condition), two simulations are performed: one with no heat pump active and,
therefore, generating a velocity field without a heat pump. The second
simulation has an active heat pump, injecting water at a different temperature
than the background groundwater temperature to generate the thermal plume.

The simulation domain covers an area of 130m$\times$130m$\times$2m, divided
into 65$\times$65$\times$1 finite volume grid cells (each grid cell is
2m$\times$2m$\times$2m) for a total of 4225 cells. The heat pump is placed at
location $(33,33)$ to be in the center of the domain. A permeability field is
generated by placing pilot points on a uniformly space grid of either
4$\times$4 (16 pilot points) or 6$\times$6 (36 pilot points) within the
simulation domain. A random value is generated for each pilot point between
1.13 $\times 10^{-7}$ and 3.77 $\times 10^{-11}$. The values at the pilot
points are interpolated onto the PFLOTRAN mesh using radial basis function
interpolation with thin-plate-splines basis function, to create the randomized
permeability field. The interpolation step from a coarse mesh of pilot points
to a fine PFLOTRAN mesh avoids sudden changes in the permeability, which could
occur if a random value was assigned to each mesh cell. The large difference
in orders of magnitudes creates both high and low permeability regions, that
can allow for the non-uniform velocity streamlines. The pressure gradient
boundary condition is generated by randomly selecting two values between
$[-0.0006,0.0006]$ and applying these in the $x$-direction and $y$-direction.
The magnitude of the values generates realistic groundwater Darcy velocities
and allow the direction of flow at the heat pump location to be in any
360$\degree$ direction.

Finally, the combination of each permeability field and pressure boundary
condition resulted in 800 unique input data samples. The small PFLOTRAN
simulation domain allowed for the fast generation of data, with each set of 25
samples taking approximately $180s$.

####  2.4.2 Without-heat-pump evaluation:

The input data for the CNN requires the Darcy velocities without the influence
of a heat pump in the region. This recreates the conditions of extracting the
velocities from a region of the baseline simulation domain without a heat
pump. Therefore, the heat pump mass flow rate was set to zero and the PFLOTRAN
simulation was run for a period of 365 days, where the Darcy velocities were
extracted at day 365 to obtain a pseudo steady-state solution.

####  2.4.3 With-heat-pump evaluation:

Supervised training of the network requires the correct temperature field to
evaluate the prediction accuracy. Therefore, each simulation performed in the
”Without-heat-pump evaluation” was rerun with the heat pump activated. The
heat pump mass flow rate was set to 0.05 $l/s$ and an injection temperature of
15$\degree$C, against the background temperature of 10$\degree$C.

####  2.4.4 Data pre-processing:

The accuracy of the CNN model can be greatly improved by suitable pre-
processing of the data. Firstly, the background of 10$\degree$C is subtracted
from the temperature output $T_{Target}$ such that the far away temperature is
around 0$\degree$C. Next, each data sample field (Darcy velocities $q_{x}$ and
$q_{y}$ and temperature field $T_{Target}$) is normalized to the range
$[-1,1]$ over the whole data set. To obtain the correct output, the inverse
operation is applied to the temperature prediction to obtain the final output.

####  2.4.5 Network limitations:

The current network design still suffers from some practical limitations.
Firstly, the network can only predict the groundwater temperature under the
assumption that the background temperature is completely uniform. Secondly, it
is assumed that the injection temperature is exactly 5$\degree$C higher than
the background temperature and with a constant flow rate of 0.05 l/s. However,
the purpose of this study was to determine whether CNNs are capable of being a
reasonably good surrogate model. Tackling further practical usage aspects are
planned for future studies.

###  2.5 Experimental setup

The three networks were trained using the Adam optimizer [28] for 50,000
epochs with a fixed learning rate of $0.0005$ and a batch size of $64$. Each
network was trained using $4$, $8$, $16$ and $32$ initial features, with the
total number of trainable parameters per network shown in table 1.

Table 1: Number of trainable parameters for the convolutional neural network. The number of parameters increases when increasing the number of initial features extracted in the first convolutional layer. Init. Feat. |  | TurbNetGeo | TurbNetGeo-Light | TurbNetGeo-NoSkip-Light
---|---|---|---|---
4 |  | 77,033 | 18,697 | 17,221
8 |  | 306,257 | 73,873 | 68,041
16 |  | 1,221,281 | 293,665 | 270,481
32 |  | 4,877,633 | 1,171,009 | 1,078,561

A total of 800 samples were generated, of which 650 were used for training and
150 for testing (validation). A data-driven loss function was used to compare
the network output with the actual results from the numerical simulation
training data. The mean squared error is defined as

$MSE=\frac{1}{N_{data}}\sum_{i}^{N_{data}}\left(T^{i}_{Pred}-T^{i}_{Target}\right)^{2},$
(9)

where $T^{i}_{Pred}$ is the temperature prediction for sample $i$,
$T^{i}_{Target}$ is the known solution, and $N_{data}$ is the total number of
data samples used for either training or testing.

The neural networks were all built using PyTorch and trained using an NVIDIA
GeForce RTX 3090 GPU with 24Gb RAM at the University of Stuttgart. The
training and testing data, as well as the Python code for the network is
provided in the DaRUS dataset ”Replication Data for: Geothermal-ML -
predicting thermal plume from groundwater heat pumps” [29].

##  3 Results and Discussion

The following section provides the results and discussion for all three
networks: TurbNet-Geo (TNG), TurbNet-Geo-Light (TNG-L) and TurbNet-Geo-NoSkip-
Light (TNG-NS-L). The training and testing loss for the three network
architectures and varying initial features are shown in table 2 at the
50,000th epoch.

Table 2: Training and testing loss of all networks with varying number of initial features. The three networks were trained with 4, 8, 16 and 32 initial features (’Init. Feat.’), which varies the total number of training parameters in the network. The training loss (’Loss’) and the testing loss (’Test Loss’) are evaluated at the 50,000th epoch.  | TurbNetGeo |  | TurbNetGeo-Light |  | TurbNetGeo-NoSkip-Light
---|---|---|---|---|---
Init. Feat. | Loss | Test Loss |  | Loss | Test Loss |  | Loss | Test Loss
4 |  1.19 $\cdot 10^{-5}$ | 0.0659 |  |  1.38 $\cdot 10^{-5}$ | 0.0556 |  |  2.40 $\cdot 10^{-5}$ | 0.0224
8 |  6.22 $\cdot 10^{-6}$ | 0.0455 |  |  3.93 $\cdot 10^{-6}$ | 0.0423 |  |  1.26 $\cdot 10^{-6}$ | 0.0173
16 |  2.22 $\cdot 10^{-7}$ | 0.0812 |  |  3.31 $\cdot 10^{-7}$ | 0.0318 |  |  3.96 $\cdot 10^{-7}$ | 0.0171
32 |  2.16 $\cdot 10^{-7}$ | 0.0433 |  |  2.17 $\cdot 10^{-7}$ | 0.0270 |  |  2.66 $\cdot 10^{-7}$ | 0.0159

The training loss for all three networks with 4 and 8 initial features are at
least an order of magnitude larger than with 16 and 32 initial features. The
TNG-NS-L has the largest training loss at 50,000 epochs for both 16 and 32
initial features. The TNG network has the lowest training loss for both 16 and
32 initial features, but only marginally lower with 4 times the number of
trainable parameters than the other two networks, making each network
evaluation more expensive. Despite having fewer trainable parameters than the
other networks, the TNG-NS-L has the lowest testing loss. Overall, the three
network architectures show similar training losses for the same number of
initial features. This indicates that even the smallest network architecture
is already expressive enough to capture the prediction task at hand.

The training loss across all 50,000 epochs for all three network architectures
with 16 and 32 initial features is shown in Figure 4. For each case, the loss
value is only plotted for every 500th epoch, with no smoothing performed on
the data. Examining the training loss over time in Figure 4, there is little
difference between the networks with an equivalent number of 16 or 32 initial
features. It is clear that the number of initial features influences the
training loss more than the number of network layers.

The testing loss, measured every 10,000th epoch, is shown for all three
network architectures with 16 and 32 initial features in Figure 5. Minor, if
any, improvement in the testing loss is observed after the first test at
10,000 epochs, indicating that additional training does not aid in improving
the real-world capability of the network on predicting the thermal plume on
unseen Darcy velocity data. This would perhaps require more training and
testing data to improve, or a fundamental change to the network is required to
reduce the testing loss.

Figure 4:  Training loss for all three network architectures TurbNet-Geo
(’TNG’), TurbNet-Geo-Light (’TNG-Light’) and TurbNet-NoSkip-Light (’TNG-
NoSkip’), were evaluated for 16 and 32 initial features, denoted as ’-16’ and
’-32’, respectively. The training loss is plotted at every 200th epoch.
Figure 5: Testing loss for all three network architectures TurbNet-Geo
(’TNG’), TurbNet-Geo-Light (’TNG-Light’) and TurbNet-NoSkip-Light (’TNG-
NoSkip’), were evaluated for 16 and 32 initial features, denoted as ’-16’ and
’-32’, respectively. The testing loss is plotted at every 5000th epoch.

For each network, a total of 150 test samples were evaluated and categorized
according to the maximum error $\epsilon_{max}$ across all pixels:

  1. (i)

good: $\epsilon_{max}<$ 0.5

  2. (ii)

medium: 0.5 $<\epsilon_{max}<$ 1.0

  3. (iii)

bad: 1.0 $<\epsilon_{max}$

The total number of predictions categorized as ”good”, ”medium” and ”bad” for
each network and varying number of initial features, is shown in table 3. The
TNG network performs only moderately better than TNG-L, with 24 and 22 ”good”
prediction samples, respectively. The TNG-NS-L has 15 ”good” predictions,
while the TNG is a might larger network with more trainable parameters. To
reduce the network size as much as possible while maintaining reasonably good
predictions, we focus on the TNG-L and TNG-NS-L networks as the preferred
networks for further analysis. Comparing these two, the TNG-L has more test
samples categorized as good for 16 and 32 initial features and more samples
categorized as medium for 16 initial features. Even though the testing loss is
lower for the TNG-NS-L network, there may be few pixels where the error is
large, forcing test samples to be categorized as ”medium” or ”bad”. However,
we select the TNG-L with 32 initial features as the current default network
for the rest of the analysis and compare the analysis to this network.

Table 3: Number of samples classified as ’good’, ’medium’ and ’bad’ predictions for all three network architectures. Any prediction with the maximum absolute error $|\epsilon_{max}|$¡ 0.5$\degree$C is defined as ’good’, any prediction with 0.5$\degree$C ¡ $|\epsilon_{max}|$ ¡ 1$\degree$C is defined as ’medium’ and all others are defined as bad. | TurbNetGeo | TurbNetGeo-Light | TurbNetGeo-NoSkip-Light
---|---|---|---
Init. Feat. | Good | Medium | Bad | Good | Medium | Bad | Good | Medium | Bad
4 | 0 | 5 | 145 | 0 | 3 | 147 | 0 | 1 | 149
8 | 1 | 50 | 99 | 1 | 28 | 121 | 1 | 43 | 106
16 | 7 | 81 | 62 | 14 | 56 | 80 | 8 | 64 | 78
32 | 24 | 66 | 60 | 22 | 63 | 65 | 15 | 69 | 66

Table 3 categorizes the predictions for the test set into ”good”, ”medium” and
”bad” predictions based on the maximum point-wise error. For the TurbNetGeo-
Light (TNG-L) and TurbNetGeo-NoSkip-Light (TNG-NS-L) networks with 16 and 32
initial features, box plots (box and whisker plot) of the per-pixel error
magnitude across all 150 test samples are plotted in Figure 6(A) - (D). Each
box provides the median $Q_{2}$ (middle black line), lower quartile $Q_{1}$
and upper quartile $Q_{3}$ (box edges). The top and bottom edges (whiskers)
define the 1.5$\times$ interquartile range (IQR), i.e., 1.5$\cdot(Q_{3}$ \-
$Q_{1})$ values of the error magnitude. For each test (A) to (D), the box
plots were generated using values where the error was more than 0.2$\degree$C
to ignore minor fluctuations of the prediction output error in regions where
the background temperature was 10$\degree$C, i.e., to ignore errors far away
from the thermal plume. Furthermore, all points within an inner region of 7
pixels from the heat pump injection site in each direction are shown in the
’Inner’ plots and all other points are shown in the ’Outer’ plots. This
highlights if the error is, in general, larger near the heat pump injection
site or further downstream the thermal plume. For all box plots, the median
error is below 0.4$\degree$C, indicating that at least half of the pixels
where the error is above 0.2$\degree$C is also below 0.5$\degree$C (the ”good”
criterion). The largest upper interquartile point occurs for the TNG-L-16
’Inner’ with a magnitude of 0.73$\degree$C, meaning 75% of error values are
beneath this value. As these plots also include the ”bad” predictions, this
shows that the networks are able to provide reasonably good predictions, and
only few pixels have a large error that cause the ”bad” categorization. The
’Inner’ box plots have a larger interquartile range and top whisker than the
’All’ plot as the error tends to be slightly larger in the middle of the
domain, leaving fewer small error values to pull the box plot down.

(a) TNG-L-16

(b) TNG-L-32

(c) TNG-NS-L-16

(d) TNG-NS-L-16

Figure 6: Per-pixel error magnitude box plots for the TurbNet-Geo-Light and
TurbNet-Geo-NoSkip-Light network architectures with 16 and 32 initial features
for all 150 test samples. Each box provides the median $Q_{2}$, lower quartile
$Q_{1}$ and upper quartile $Q_{3}$. The top and bottom edges (whiskers) define
the 1.5 interquartile range (1.5$\cdot(Q_{3}$ \- $Q_{1})$) values of the error
magnitude. Only error values larger than 0.2$\degree$C were added to the
dataset for the box plots.

Examining the training and testing loss conveys how well the network performs
overall. Additionally, the box plots tells use how many pixels are considered
good, medium and bad. However, it does not indicate where or how the greatest
errors occur, i.e., whether it occurs at the GWHP location itself or if the
thermal plume is unable to follow the streamlines. Therefore, a qualitative
assessment of the network’s performance is required to gain further insight.
Four different network predictions for the TNG-L network with 32 initial
features from the ”good”, ”medium” and ”bad” categories are shown in Figure 7,
Figure 8 and Figure 9, respectively. Finally, the difference between the
TurbNetGeo-Light (TNG-L) and the TurbNetGeo-NoSkip-Light (TNG-NS-L) with 32
initial features from the ”bad” category are exemplified on four test samples
in Figure 10 and Figure 11.

Figure 7: TurbNet-Geo-Light: ”good” network prediction for four test samples
with 32 initial features. The test sample numbers listed from top to bottom:
35, 96, 100 and 136.

In Figure 7 the maximum error occurs at the heat pump location for all four
samples. However, the second sample (second row) has very small errors in the
thermal plume in comparison to the others. The top and bottom samples indicate
a good ability to follow the streamline as its path changes. The similar
looking middle samples show how the thermal plume is able to spread out when
the streamlines diverge (second row), whereas the thermal plume remains narrow
when the streamlines remain straight (third row). As expected, the ”good”
category performs well in general.

Figure 8: TurbNet-Geo-Light: ”medium” network prediction for four test samples
with 32 initial features. The test sample numbers listed from top to bottom:
16, 76, 99 and 12.

Four test samples from the ”medium” category for the TNG-L with 32 initial
features are shown in Figure 8. The top two samples have a very low error
magnitude in the thermal plume, and was categorized into ”medium” due to a
small error spike at the heat pump location. The large error occurs when the
GWHP location in the plume prediction is shifted one or two pixels compared to
the target solution, shifting the maximum temperature pixel. If shifted in the
wrong direction, i.e., if the plume extends downwards but the GWHP pixel
location is shifted upwards, a large difference exists at this shifted
location only and artificially inflates the maximum error value. However, the
error values within the thermal plumes themselves are well within acceptable
limits. The bottom two samples show that the plume can morph with the
streamlines, but not well enough to keep the error below 0.5$\degree$C.

Figure 9: TurbNet-Geo-Light: ”bad” network prediction for four test samples
with 32 initial features. The test sample numbers listed from top to bottom:
20, 52, 66 and 139.

Four test samples from the ”bad” category for the TNG-L with 32 initial
features are shown in Figure 9. Each sample indicates various problems with
some predictions. The top sample cannot accurately capture the temperature far
downstream, whereas the second sample does not follow the streamline at all,
but cuts across it instead. This is in contrast to our assumption that the
thermal plume follows the streamline. The third sample cannot capture the high
temperature at the heat pump location, where the plume is also relatively wide
due to having a low Darcy velocity at the heat pump. This causes the higher
temperature to diffuse outwards instead of being dragged downstream. Finally,
the bottom sample is unable to capture the thermal profile near the heat pump.

Another four test samples from the ”bad” category for the TNG-NS-L (first and
third image) and TNG-L (second and fourth image) with 32 initial features,
comparing the ability of the two networks, are shown in Figure 10 and Figure
11, each with two samples. For sample 20 in Figure 10, the TNG-NS-L network
provides a better prediction than the TNG-L network, which cannot predict the
entire plume correctly. The largest error occurs further downstream of the
GWHP. However, the TNG-L network provides a better prediction of the
downstream plume compared to the TNG-NS-L for sample 51.

In Figure 11, the prediction for the first sample is similar for both
networks, where the largest error is close to the GWHP. However, the TNG-NS-L
network has a wider plume and does not quite follow the streamline as well as
the TNG-L network. The final prediction is difficult for both networks. Both
fail to predict the maximum temperature directly downstream of the GWHP and
also over-predict the width of the plume.

By comparing the two networks, both have their own positives and negatives.
The TNG-NS-L sometimes outperforms the TNG-L, but it is not clear that it is
necessarily better for all cases. The only clear benefit of the TNG-NS-L would
be the simpler design and fewer trainable parameters for training. The
provided examples are only a small set of samples that are available for
comparison and were chosen from the set of 150 samples to highlight problems
that may occur for the network prediction. Even for the ”bad” category, the
predictions are not completely unusable, but there are edge cases where the
prediction would not be suitable to use in place of the high-fidelity solver.
However, the CNN is suitable as an initial pre-processing step for the high-
fidelity optimization problem, where highly interacting heat pumps can be
identified. Therefore, the recommended network is the TNG-L with 32 initial
features, followed by TNG-NS-L with 32 initial features. This network is
suitable to implement inside the online evaluation tool of Fig to complement
the LAHM.

(a) TNG-NS-L: sample 20

(b) TNG-L: sample 20

(c) TNG-NS-L: sample 51

(d) TNG-L: sample 51

Figure 10: TurbNet-Geo-Light versus TurbNet-Geo-NoSkip-Light: ”bad” network
prediction for two test samples with 32 initial features.

(a) TNG-NS-L: sample 137

(b) TNG-L: sample 137

(c) TNG-NS-L: sample 50

(d) TNG-L: sample 50

Figure 11: TurbNet-Geo-Light versus TurbNet-Geo-NoSkip-Light: ”bad” network
prediction for two samples with 32 initial features.

##  4 Conclusion

In this paper, we introduced a novel method for predicting the thermal plume
downstream propagation from open-loop groundwater heat pump injection wells.
Conventional methods for modeling heat pumps are either too expensive or lack
a suitable level of accuracy. By understanding the governing equations for the
subsurface water temperature and studying readily available simulation data, a
seemingly simplistic convolutional neural network was built that accepts the
subsurface Darcy velocities as input and outputs the thermal profile due to
the presence of a heat pump. The network was based on a modified U-Net
architecture to create three network variants, one comprised of 6 layers with
skip connections (TurbNet-Geo), one with 4 layers with skip connections
(TurbNet-Geo-Light) and one with 4 layers without skip connections (TurbNet-
NoSkip-Light).

Training and testing data was generated for small 2D domains, by creating
random permeability and pressure gradient boundary conditions, and performing
numerical groundwater simulations with PFLOTRAN to obtain the Darcy velocities
and temperature field. A total of 800 input-output samples were generated,
with 650 samples used to train the network and 150 samples used to test the
network prediction.

The large TurbNet-Geo architecture was found to be only marginally more
accurate than the smaller network architectures, but requiring nearly 4 times
the number of trainable parameters. The prediction results were classified
into good, medium and bad according to a specified criterion of
$|\epsilon_{max}|<$ 0.5, 0.5 $<|\epsilon_{max}|<$ 1.0 and 1.0
$<|\epsilon_{max}|$, respectively. It was found that only a few pixels in the
prediction caused the predictions to be classified as medium or bad. Notably,
75% of all pixels where the error was larger than 0.2$\degree$C were also
below 0.73$\degree$C, within the ”medium” classification. Examining the
thermal plume profiles of each network, each architecture showed good
agreement between the prediction and the ground truth for most of the test
samples. There were a few samples where the network struggled to reasonably
determine the shape of the plume following the streamlines. However, this was
limited to a few samples only and were presented in the discussion. Overall,
the networks performed well as a fast surrogate model alternative to the high-
fidelity solver. Future work will include training the network on larger
datasets, larger domains, and more complex 3D domains where generating the
training data is significantly more expensive.

##  5 Funding

Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
under Germany’s Excellence Strategy—EXC 2075—390740016. We acknowledge the
support by the Stuttgart Center for Simulation Science (SimTech).

##  6 Conflicts of Interest

The funders had no role in the design of the study, in the collection,
analyses, or interpretation of data; in the writing of the manuscript, or in
the decision to publish the results.

##  7 References

## References

  * [1] Commission E, Centre JR, Economidou M, Ringel M, Valentova M, Zancanella P, et al.  National energy and climate plans for 2021-2030 under the EU energy union : assessment of the energy efficiency dimension.  Publications Office; 2020.  Available from: https://op.europa.eu/en/publication-detail/-/publication/2c751d22-4f07-11eb-b59f-01aa75ed71a1/language-en/format-PDF/source-271999318.
  * [2] Attard G, Bayer P, Rossier Y, Blum P, Eisenlohr L.  A novel concept for managing thermal interference between geothermal systems in cities.  Renewable Energy. 2020;145:914-24.  Available from: https://doi.org/10.1016/j.renene.2019.06.095.
  * [3] García-Gil A, Moreno MM, Schneider EG, Marazuela MÁ, Abesser C, Lázaro JM, et al.  Nested shallow geothermal systems.  Sustainability (Switzerland). 2020;12(12).
  * [4] Meng B, Vienken T, Kolditz O, Shao H.  Evaluating the thermal impacts and sustainability of intensive shallow geothermal utilization on a neighborhood scale : Lessons learned from a case study.  Energy Conversion and Management. 2019;199(July):111913.  Available from: https://doi.org/10.1016/j.enconman.2019.111913.
  * [5] Beck M, Bayer P, de Paly M, Hecht-Méndez J, Zell A.  Geometric arrangement and operation mode adjustment in low-enthalpy geothermal borehole fields for heating.  Energy. 2013;49(1):434-43.  Available from: http://dx.doi.org/10.1016/j.energy.2012.10.060.
  * [6] De Paly M, Hecht-Méndez J, Beck M, Blum P, Zell A, Bayer P.  Optimization of energy extraction for closed shallow geothermal systems using linear programming.  Geothermics. 2012;43:57-65.  Available from: http://dx.doi.org/10.1016/j.geothermics.2012.03.001.
  * [7] Cardoso MA, Durlofsky LJ, Sarma P.  Development and application of reduced-order modeling procedures for subsurface flow simulation.  International Journal for Numerical Methods in Engineering. 2009;77(9):1322-50.  Available from: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2453.
  * [8] Pasetto D, Guadagnini A, Putti M.  POD-based Monte Carlo approach for the solution of regional scale groundwater flow driven by randomly distributed recharge.  Advances in Water Resources. 2011;34(11):1450-63.  Available from: https://www.sciencedirect.com/science/article/pii/S0309170811001321.
  * [9] Li X, Chen X, Hu BX, Navon IM.  Model reduction of a coupled numerical model using proper orthogonal decomposition.  Journal of Hydrology. 2013;507:227-40.  Available from: https://www.sciencedirect.com/science/article/pii/S0022169413006562.
  * [10] Vinuesa R, Brunton SL.  Enhancing computational fluid dynamics with machine learning.  Nat Comput Sci. 2022;2:358-66.  Available from: https://doi.org/10.1038/s43588-022-00264-7.
  * [11] Chen W, Wang Q, Hesthaven JS, Zhang C.  Physics-informed machine learning for reduced-order modeling of nonlinear problems.  Journal of Computational Physics. 2021;446:110666.  Available from: https://www.sciencedirect.com/science/article/pii/S0021999121005611.
  * [12] Laubscher R. Utilization of artificial neural networks to resolve chemical kinetics in turbulent fine structures of an advanced CFD combustion model; 2017\.
  * [13] du Toit P. An artificial intelligence approach for biomass devolatilisation in an industrial CFD model with advanced turbulence-chemistry interaction; 2018\.
  * [14] Thuerey N, Weißenow K, Prantl L, Hu X.  Deep learning methods for Reynolds-averaged Navier-Stokes simulations of airfoil flows.  American Institute of Aeronautics and Astronautics. 2019;58:686 707.
  * [15] Wang Y, Lin G.  Efficient deep learning techniques for multiphase flow simulation in heterogeneous porousc media.  Journal of Computational Physics. 2020;401:108968.  Available from: https://www.sciencedirect.com/science/article/pii/S0021999119306734.
  * [16] Tang M, Liu Y, Durlofsky LJ.  Deep-learning-based surrogate flow modeling and geological parameterization for data assimilation in 3D subsurface flow.  Computer Methods in Applied Mechanics and Engineering. 2021;376:113636.  Available from: https://www.sciencedirect.com/science/article/pii/S0045782520308215.
  * [17] Gao H, Sun L, Wang JX.  Super-resolution and denoising of fluid flow using physics-informed convolutional neural networks without high-resolution labels.  Physics of Fluids. 2021;33(7).
  * [18] Raissi M, Perdikaris P, Karniadakis GE.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics. 2019;378:686-707.  Available from: https://www.sciencedirect.com/science/article/pii/S0021999118307125.
  * [19] Zhu Y, Zabaras N, Koutsourelakis PS, Perdikaris P.  Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data.  Journal of Computational Physics. 2019;394:56-81.  Available from: https://www.sciencedirect.com/science/article/pii/S0021999119303559.
  * [20] Kashefi A, Rempe D, Guibas LJ.  A point-cloud deep learning framework for prediction of fluid flow fields on irregular geometries.  Physics of Fluids. 2021;33(2).
  * [21] Laubscher R.  Simulation of multi-species flow and heat transfer using physics-informed neural networks Simulation of multi-species flow and heat transfer using physics-informed neural networks.  Phys Fluids. 2021;33(July).
  * [22] Pavlova A, Hansen J, Obermeyer H, Pavlova I.  Geothermal heat in a heat pump use.  IOP Conference Series: Earth and Environmental Science. 2016 9;43:012025.  Available from: https://doi.org/10.1088/1755-1315/43/1/012025.
  * [23] Kinzelbach W.  Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser.  2nd ed. Munich: Oldenbourg; 1987.
  * [24] PFLOTRAN. PFLOTRAN Documentation: Theory Guide - TH Mode (Thermal-Hydrologic); 2022\.  Available from: https://documentation.pflotran.org/theory_guide/mode_th.html.
  * [25] Tang M, Liu Y, Durlofsky LJ.  A deep-learning-based surrogate model for data assimilation in dynamic subsurface flow problems.  Journal of Computational Physics. 2020;413:109456.
  * [26] Ronneberger O, Philipp F, Thomas B.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  In: Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015. Cham: Springer International Publishing; 2015. p. 234-41.
  * [27] Hammond GE, Lichtner PC, Mills RT.  Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN.  Water Resources Research. 2014;50:208-28.
  * [28] Kingma DP, Ba J.  Adam: A Method for Stochastic Optimization. 2014.  Available from: https://arxiv.org/abs/1412.6980.
  * [29] Davis K, Schulte M.  Replication Data for: Geothermal-ML - predicting thermal plume from groundwater heat pumps. 2022.  Available from: https://doi.org/10.18419/darus-3184.


[a]L. P. Taylor

# Design and performance of the prototype Schwarzschild-Couder Telescope
camera

###### Abstract

The Cherenkov Telescope Array (CTA) is the next-generation ground-based
observatory for very-high-energy gamma-ray astronomy. An innovative 9.7 m
aperture, dual-mirror Schwarzschild-Couder Telescope (SCT) design is a
candidate design for CTA Medium-Sized Telescopes. A prototype SCT (pSCT) has
been constructed at the Fred Lawrence Whipple Observatory in Arizona, USA. Its
camera is currently partially instrumented with 1600 pixels covering a field
of view of 2.7 degrees square. The small plate scale of the optical system
allows densely packed silicon photomultipliers to be used, which combined with
high-density trigger and waveform readout electronics enable the high-
resolution camera. The camera’s electronics are capable of imaging air shower
development at a rate of one billion samples per second. We describe the
commissioning and performance of the pSCT camera, including trigger and
waveform readout performance, calibration, and absolute GPS time stamping. We
also present the upgrade to the camera, which is currently underway. The
upgrade will fully populate the focal plane, increasing the field of view to 8
degree diameter, and lower the front-end electronics noise, enabling a lower
trigger threshold and improved reconstruction and background rejection.

##  1 Introduction

Figure 1: The full pSCT camera is comprised of 177 modules divided into 9
sectors. Each sector can hold up to 25 modules with the sectors in each corner
containing fewer modules such that the final focal plane will be roughly
circular. Each module contains 64 image pixels with the full camera containing
11,328 pixels. Currently, only the center sector is instrumented. 24 modules
are installed, with the central module not yet installed to accommodate
optical alignment hardware.

Earth’s atmosphere is opaque to very-high-energy (VHE) photons, however these
photons initiate air showers upon hitting the atmosphere. The shower
constituents move faster than the speed of light in the atmosphere, resulting
in Cherenkov radiation which can be detected at the ground by Imaging
atmospheric Cherenkov Telescopes (IACTs). The Cherenkov Telescope Array (CTA)
is an IACT array with one site in each hemisphere. CTA will feature three
telescope sizes and cover energies between 20 GeV and 300 TeV [1]. Most IACTs
use Davies-Cotton (single-mirror) optics. However, the prototype
Schwarzschild-Couder Telescope (pSCT), a candidate for the CTA Medium Sized
Telescope, will use dual-mirror optics. The Schwarzschild-Couder optics
produce a wide field of view, an excellent optical point spread function, and
a smaller plate scale when compared to single-mirror designs [2, 3].

This small plate scale allows the use of silicon photomultipliers (SiPMs) in
place of traditional photomultiplier tubes. SiPMs are smaller than traditional
PMTs enabling improved image resolution. The current pSCT is instrumented with
25 modules (1600 pixels) and has a 2.7∘ field of view. An upgrade to the pSCT
is underway which is expected to be complete by the end of 2022. The upgrade
will include improvements to the camera consisting of an update of the camera
module front-end electronics (FEE), the camera backend electronics, and the
module SiPMs as well as fully populating the focal plane (increasing the
number of pixels to 11,328 and the field of view to 8∘ diameter.

##  2 Commissioning and Performance

Commissioning of the pSCT camera is an ongoing process and in the last two
years, significant improvements have been made in the quality of pSCT data
taking. Among these are validation of the trigger path mapping, increasing the
trigger path dynamic range, and the inclusion of all modules in the trigger.
Concurrent with commissioning work, initial observations of the Crab Nebula
were made. The observation and detection of the Crab Nebula by the pSCT are
described in detail in [4, 5].

###  2.1 Validation of the Trigger Path

Figure 2: Circuit diagram for one trigger pixel. [6] In TARGET7, four (2x2)
image pixels are amplified individually via two inverting amplifiers (with
Vofs1 and Vofs2 as reference) and then summed with an inverting summing
amplifier (with PMTref4 as reference) to produce the signal for one trigger
pixel. The output of the summing inverting amplifier (Node C) originally had a
desired set-point of 1.25 V at baseline. This desired set-point is achieved by
tuning the reference values Vofs1, Vofs2, and PMTref4 such that Node C reaches
this value for a baseline signal. This desired set-point was changed to 2.1 V
resulting in a greater dynamic range for our expected signal (downgoing at
Node C). Note that due to inversion of the signal along the amplifier chain, a
higher Thresh value corresponds to a lower physical threshold.

(a) Design of the physical mask. Each module is represented by a group of 64
(8x8) image pixels (dots). Red empty circles highlight the image pixels used
to align the mask using the data path image. Yellow filled circles are the
trigger pixels used for checking the trigger pattern. Only one group of
trigger pixels is capable of satisfying the three-fold coincidence required to
produce a backplane trigger.

(b) Realization of the physical mask used to cover the focal plane.
Constructed from Cardboard and illuminated from behind. Each dot is the center
of a 6x6 $\mathrm{mm}^{2}$ image pixel; 4 image pixels make up a trigger
pixel. The corner slits are used for visual alignment, while the single-pixel
openings in each corner aid with alignment using the data path. Ten trigger
pixels were chosen to remain uncovered.

(c) Data-path image taken with camera mask installed on the camera. The color
scale is proportional to the signal amplitude in each pixel. Unmasked pixels
are clearly visible in yellow/green and the reference pixels in each corner
are well aligned. All pixels visible in the data path are also visible in the
trigger hit map in Figure 3(d).

(d) A binary image of module triggers (hit map). A yellow square indicates the
pixel has triggered while a white square indicates that it has not. The hit
map was read out concurrently with the event image in Figure 3(c). Only the
cluster of three trigger pixels on the left side passes the 3-fold coincidence
logic and produces a backplane trigger.

Figure 3: An asymmetric cardboard mask was constructed and attached to the
camera focal plane. The mask was designed such that there was only one
grouping of pixels capable of producing a camera trigger (which uses three-
fold-coincidence). A moderate intensity flasher was used to illuminate the
camera through the mask. Event images from the data path and hit maps from the
trigger path both match the mask pattern.

The pSCT trigger system contains two trigger levels: module triggers and
backplane triggers (when integrated into an array a third telescope trigger
will also be included). The module level trigger system (Figure 2) monitors
groups of four (2x2) image pixels (called a trigger pixel) for a threshold
(Thresh) crossing of their analog sum. If this value is crossed for any
trigger pixel a module trigger is sent to the backplane. The backplane then
collects these module trigger signals from each of the camera modules. The FEE
programmable pulse width (currently set to $\sim$10 ns) determines the
coincidence resolving time of the pattern trigger. A module trigger overlap of
approximately 1.5 ns is sufficient to generate a coincidence.

The backplane then forms backplane level triggers via a single trigger FPGA to
form coincidences from the 400 module trigger inputs (16 module trigger inputs
for each of the 25 modules). When three adjacent trigger pixels produce a
module trigger at the same time, a backplane trigger is produced, prompting
the camera to read out the event. A pixel is considered adjacent if it is
orthogonal-adjacent or diagonal-adjacent (i.e. a pixel can have up to eight
adjacent pixels). A timestamp along with a hit map (a binary image of module
triggers) is attached to the backplane trigger.

In order to validate this trigger system, a physical mask was constructed out
of cardboard and affixed to the camera focal plane. The mask was aligned
visually with the focal plane using three slits at the corners of the focal
plane. Four small holes were placed at each corner for data-path alignment.
The mask was designed such that ten trigger pixels would be illuminated in an
asymmetric pattern so that any mapping irregularities would be easy to spot.
It was also designed such that only one grouping of trigger pixels would
produce a backplane trigger. This grouping contained trigger pixels from three
different modules in order to ensure that cross-module triggering was
effective.

The camera was then illuminated by an LED flasher mounted on the secondary
mirror (running at 10 Hz) and allowed to trigger using the usual three-fold
coincidence. hit maps were read out concurrently. Data path images and hit
maps were then compared to the mask pattern. Figure 3 shows a single example.
The comparison showes that the data path mapping and trigger path mapping are
consistent both with each other and the physical layout of the camera.

###  2.2 Improving Dynamic Range

Figure 4: Two rate scans, comparing Thresh value to the trigger rate, each
with all modules included in the trigger and a flasher running at 10 Hz at
constant intensity. The baseline value of the output of the summing inverting
amplifier (Node C) shown in Figure 2 has been changed from 1.25 V to 2.1 V in
order to improve the dynamic range of our expected signal. The upper (blue)
rate scan was taken with the previous value of 1.25 V. The lower (red) rate
scan was taken with a new value of 2.1 V. The new configuration shows a
significant improvement in the dynamic range.

Modules contain 64 image pixels. In the trigger path these image pixels are
grouped into trigger pixels before trigger logic is applied. A trigger pixel
is comprised of the analog sum of four image pixels. The individual signals
are combined using a summing amplifier with the value PMTref4 used as
reference (see Figure 2). Nominally, the output of the summing amplifier (Node
C) can vary between 0–2.5 V. However, the range of linearity is 0.4–2.1 V. The
output of this amplifier is fed into a comparator where the signal is compared
with the Thresh value to determine if there is a trigger. Due to inversions
along the amplifier chain a high Thresh value corresponds to a low physical
threshold.

Originally, the Node C was selected to be 1.25 V at baseline, the midpoint of
the possible range. This desired output was achieved by tuning the reference
values Vofs1, Vofs2, and PMTref4 for a baseline input. Using the midpoint
value meant that both positive and negative signals had the same dynamic
range. However, pulses coming into the inverting summing amplifier are
unidirectional, meaning half of this dynamic range was not being used. For
this reason we decided to raise the baseline value of Node C to 2.1 V as that
was the highest value at which the summing amplifier maintained linearity.
This change effectively extended the trigger dynamic range for our expected
signal. Two rate scans (comparing trigger rate vs Thresh value) were taken in
order to illustrate this effect. Both rate scans were taken with the flasher
running at 100 Hz. One was taken with the old settings but with all modules in
the trigger (when typically only US modules would be included - See section
2.3) and one with the new settings. Figure 4 shows these rate scans. Prior to
the adjustment the camera exhibited a very narrow trigger dynamic range,
quickly entering the baseline noise. The new setting improved the dynamic
range significantly, clearly showing the flasher plateau.

###  2.3 Including All Modules in the Trigger

Currently the camera includes two types of SiPMs. Fifteen of the modules (US)
use Hamamatsu SiPMs (model S12642-0404PA-50(X)) while the remaining nine
(INFN) use FBK NUV-HD SiPMs. [7] INFN modules have approximately twice the
gain of US modules, due to a two-times larger FEE electronics gain combined
with different SiPM gains. Because of this the conversion between Thresh value
and photoelectrons for both modules is different. This difference in threshold
meant that for the same Thresh value INFN modules would dominate the trigger.
This fact combined with the geometry of the camera, which houses INFN modules
on the edge of the focal plane, meant that it was beneficial to exclude INFN
modules from the trigger during Crab observations. In an effort to update the
trigger the Thresh to photoelectron conversion was found for both US and INFN
modules with the intention of selecting two Thresh values (one for US and one
for INFN modules) which corresponded to the same photoelectron threshold.
However, there is no threshold at which both US and INFN modules can be
included in the trigger. Rather than continue to leave the INFN modules out of
the trigger completely, it was decided that they would be included at a thresh
value which keeps them subdominant to the US modules (approximately one tenth
the trigger rate). The INFN modules were included in the trigger starting
March 17, 2021.

###  2.4 Temperature Dependent Trigger Tuning

At the module level, the trigger path consists of two inverting amplifiers and
an inverting summing amplifier (Figure 2). Node C is desired to be 2.1 V when
the baseline of the signal is at its nominal value (see Section 2.2). This
output depends on the values Vofs1, Vofs2, and PMTref4, which are the
reference values for the two inverting amplifiers and the inverting summing
amplifier, respectively. Due to slight differences in component values between
circuits, these values must be tuned to reach the desired output. Originally,
this tuning was completed once for each module in the camera in the lab and
the resulting reference values were used for all data taking.

It was discovered that temperature had a significant impact on these reference
values resulting in an apparent increase in noise as the temperature of the
modules decreased. This effect has been mitigated by implementing temperature-
dependent trigger tuning (with a precision of 0.25∘ C). Each module FEE has
four temperature sensors and the temperature of the modules is determined by
averaging these four temperature readings. Each module is tuned at a variety
of temperatures and the resulting reference values are stored in a database.

During observations a five minute temperature equilibration time occurs prior
to each data taking run. At the end of this time the temperature of each
module is determined. From this temperature the correct reference values are
found using the temperature dependent database and used during subsequent data
taking.

##  3 Camera Upgrade

(a) Charge spectrum of existing modules using Hamamatsu SiPMs and TARGET 7
electronics. An LED flasher was used to illuminate the entire module at once.
A synchronous hardware trigger was used to trigger the module. [8]

(b) Charge spectrum produced from upgraded modules using FBK SiPMs and TARGETC
electronics. Using a synchronous hardware trigger, a laser source combined
with a moving stage was used to illuminate one pixel at a time.

Figure 5: The upgraded modules have lower electronic noise and therefore show
a significant improvement in the charge spectrum resolution.

The camera upgrade will include an upgrade to both the SiPMs and front-end
electronics. At present the camera uses a mixture of FBK NUV-HD and Hamamatsu
SiPMs. After the upgrade the entire camera will use newer FBK NUV-HD SiPMs.
These SiPMs will have improved photon detection efficiency, lower temperature
dependence, and lower breakdown voltage than those modules currently used in
the camera. The new FEE modules will incorporate next-generation TARGET C
ASICs for digitization and sampling as well as a separate TARGET 5 Trigger
Extension ASIC (T5TEA) for triggering. Figure 5 shows the charge distribution
of the current and upgraded camera modules, illustrating a significant
decrease in electronic noise. To connect the SiPMs to the FEE, a new custom
ASIC called the SMART chip has been developed. The SMART chip will allow for
improved SiPM bias voltage control as well as individual current measurements.
It will also provide shaping and amplification of the analog signals. In
addition to these upgrades the number of modules will also be increased from
25 to 177 in order to fully populate the focal plane. This will increase in
the field of view from 2.7∘ square to 8.0∘ diameter.

In order to accommodate a fully populated focal plane the backplane PCB will
also be redesigned. The current backplane was designed to accommodate 32
modules. The upgraded design will accommodate up to 25 modules and a total of
nine backplanes will be installed in the camera (one for each sector). See
Figure 1. DACQ boards will be mounted on top of their associated backplane.
The heat management system will also be expanded in order to accommodate the
increased number of modules in the camera. Additionally, the LED flashers will
be upgraded, allowing for synchronous triggering with the backplanes and
improved performance. [9]

The final SCT will include a Distributed Intelligent Array Trigger (DIAT)
board which will provide GPS time-tagging as well as a synchronized 62.5 MHz
clock to both backplanes and FEEs. [10] In addition to a DIAT board the pSCT
upgrade will include the installation of a stand-alone time-tagging system.
This system is divided into two modules: A laser-diode module located in the
camera and a photo-detector/GPS module located near the camera server. A fiber
optic communications cable will connect the two modules. Fast pulses (up to 5
kHz from camera backplane triggers) will be generated by a laser diode and
then fed to a discriminator where a logic pulse is generated. A GPS time-stamp
will then be attached to each pulse using custom-developed software/firmware.
The system will provide stable timestamps for camera triggers to a GPS
accuracy of 5 ns or better.

The pSCT camera upgrade has two main objectives: maximizing the camera’s field
of view, and upgrading the underlying electronics, enabling a lower trigger
threshold and improved reconstruction and background rejection. The camera is
scheduled to be upgraded by December of 2022.

## References

  * [1] A. Acharyya et al., _Monte carlo studies for the optimisation of the cherenkov telescope array layout_ , _Astroparticle Physics_ 111 (2019) 35.
  * [2] V. Vassiliev and S. Fegan, _Schwarzschild-couder two-mirror telescope for ground-based $\gamma$-ray astronomy_, in _30th International Cosmic Ray Conference_ , vol. 3, pp. 1445–1448, PoS (ICRC2008) 782, 2008.
  * [3] D. Ribeiro, _Prototype schwarzschild-couder telescope for the cherenkov telescope array: Commissioning the optical system_ , in _these proceedings_ , PoS (ICRC2021) 717, 2021.
  * [4] B.A.W. Mode, _Detection of the crab nebula by the prototype schwarzschild-couder telescope_ , in _these proceedings_ , PoS (ICRC2021) 830, 2021.
  * [5] C. Adams et al., _Detection of the Crab Nebula with the 9.7 m prototype Schwarzschild-Couder telescope_ , _Astroparticle Physics_ 128 (2021) 102562.
  * [6] M. Wachtendonk, _Characterization and Reduction of Noise in the Prototype Schwarzschild-Couder Telescope Camera_ , Master’s thesis, University of Wisconsin - Madison, 2018.
  * [7] A. Gola et al., _Nuv-sensitive silicon photomultiplier technologies developed at fondazione bruno kessler_ , _Sensors_ 19 (2019) 308\.
  * [8] M. Winter, _Characterization of the Charge Spectrum of the CTA pSCT_ , Master’s thesis, University of Wisconsin - Madison, 2017.
  * [9] C. Adams et al., _Upgrading the Prototype Schwarzschild-Couder Telescope Camera to a Wide-Field, High-Resolution Instrument_ , in _36th International Cosmic Ray Conference_ , PoS (ICRC2019) 742, 2019.
  * [10] H. Dickinson et al., _An image-based array trigger for imaging atmospheric Cherenkov telescope arrays_ , _Nuclear Instruments and Methods in Physics Research A_ 891 (2018) 6 [1802.05715].

## Full Authors List: CTA SCT Consortium

C. B. Adams1, G. Ambrosi2, M. Ambrosio3, C. Aramo3, P. I. Batista4, W.
Benbow5, B. Bertucci2,6, E. Bissaldi7,8, M. Bitossi9, A. Boiano3, C.
Bonavolontà3,10, R. Bose11, A. Brill1, A. M. Brown12, J. H. Buckley11, R. A.
Cameron13, M. Capasso14, M. Caprai2, C. E. Covault15, D. Depaoli16,17, L. Di
Venere7,8, M. Errando11, S. Fegan18, Q. Feng14, E. Fiandrini2,19, A.
Furniss20, A. Gent21, N. Giglietto7,8, F. Giordano7,8, R. Halliday22, O.
Hervet23, T. B. Humensky24, S. Incardona25,26, M. Ionica2, W. Jin27, D.
Kieda28, F. Licciulli8, S. Loporchio7,8, G. Marsella25,26, V. Masone3, K.
Meagher21,29, T. Meures29, B. A. W. Mode29, S. A. I. Mognet30, R. Mukherjee14,
A. Okumura31, N. Otte21, F. R. Pantaleo7,8, R. Paoletti32,9, G. Pareschi33, F.
Di Pierro17, E. Pueschel4, D. Ribeiro1, L. Riitano29, E. Roache5, D. Ross34,
J. Rousselle35, A. Rugliancich9, M. Santander27, R. Shang36, L. Stiaccini32,9,
H. Tajima31, L. P. Taylor29, L. Tosti2, G. Tovmassian37, G. Tripodo25,26, V.
Vagelli38,2, M. Valentino10,3, J. Vandenbroucke29, V. V. Vassiliev36, J. J.
Watson39, R. White40, D. A. Williams23, A. Zink41,

1Physics Department, Columbia University, New York, NY 10027, USA. 2INFN Sezione di Perugia, 06123 Perugia, Italy. 3INFN Sezione di Napoli, 80126 Napoli, Italy. 4Deutsches Elektronen-Synchrotron, Platanenallee 6, 15738 Zeuthen, Germany. 5Center for Astrophysics | Harvard & Smithsonian, Cambridge, MA 02138, USA. 6Dipartimento di Fisica e Geologia dell’Università degli Studi di Perugia, 06123 Perugia, Italy. 7Dipartimento Interateneo di Fisica dell’Università e del Politecnico di Bari, 70126 Bari, Italy. 8INFN Sezione di Bari, 70125 Bari, Italy. 9INFN Sezione di Pisa, 56127 Pisa, Italy. 10CNR-ISASI, 80078 Pozzuoli, Italy. 11Department of Physics, Washington University, St. Louis, MO 63130, USA. 12Dept. of Physics and Centre for Advanced Instrumentation, Durham University, Durham DH1 3LE, United Kingdom. 13Kavli Institute for Particle Astrophysics and Cosmology, SLAC National Accelerator Laboratory, Stanford University, Stanford, CA 94025, USA. 14Department of Physics and Astronomy, Barnard College, Columbia University, NY 10027, USA. 15Department of Physics, Case Western Reserve University, Cleveland, Ohio 44106, USA. 16Dipartimento di Fisica dell’Università degli Studi di Torino, 10125 Torino, Italy. 17INFN Sezione di Torino, 10125 Torino, Italy. 18LLR/Ecole Polytechnique, Route de Saclay, 91128 Palaiseau Cedex, France. 19Y. 20Department of Physics, California State University - East Bay, Hayward, CA 94542, USA. 21School of Physics & Center for Relativistic Astrophysics, Georgia Institute of Technology, Atlanta, GA 30332-0430, USA. 22Dept. of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA. 23Santa Cruz Institute for Particle Physics and Department of Physics, University of California, Santa Cruz, CA 95064, USA. 24Science Department, SUNY Maritime College, Throggs Neck, NY 10465. 25Dipartimento di Fisica e Chimica "E. Segrè", Università degli Studi di Palermo, via delle Scienze, 90128 Palermo, Italy. 26INFN Sezione di Catania, 95123 Catania, Italy. 27Department of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA. 28Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA. 29Department of Physics and Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin, Madison, WI 53706, USA. 30Pennsylvania State University, University Park, PA 16802, USA. 31Institute for Space–Earth Environmental Research and Kobayashi–Maskawa Institute for the Origin of Particles and the Universe, Nagoya University, Nagoya 464-8601, Japan. 32Dipartimento di Scienze Fisiche, della Terra e dell’Ambiente, Università degli Studi di Siena, 53100 Siena, Italy. 33INAF - Osservatorio Astronomico di Brera, 20121 Milano/Merate, Italy. 34Space Research Centre, University of Leicester, University Road, Leicester, LE1 7RH, United Kingdom. 35Subaru Telescope NAOJ, Hilo HI 96720, USA. 36Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA. 37Instituto de Astronomía, Universidad Nacional Autónoma de México, Ciudad de México, Mexico. 38Agenzia Spaziale Italiana, 00133 Roma, Italy. 39Deutsches Elektronen-Synchrotron, Platanenallee 6, D-15738 Zeuthen, Germany. 40Max-Planck-Institut für Kernphysik, P.O. Box 103980, 69029 Heidelberg, Germany. 41Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen Centre for Astroparticle Physics, D 91058 Erlangen, Germany.


IEEEexample:BSTcontrol

# DaPPA: A Data-Parallel Framework for Processing-in-Memory Architectures

Geraldo F. Oliveira ⋆ Alain Kohli ⋆ David Novo ‡ Juan Gómez-Luna ⋆ Onur Mutlu
⋆

⋆ETH Zürich ‡LIRMM Univ. Montpellier  CNRS

##  1 Motivation & Problem

The increasing prevalence and growing size of data in modern applications have
led to high costs for computation in traditional _processor-centric computing_
systems. To mitigate these costs, the _processing-in-memory_ (PIM) [1, 2, 3,
4, 5, 6] paradigm moves computation closer to where the data resides, reducing
the need to move data between memory and the processor. Even though the
concept of PIM has been first proposed in the 1960s [7, 8], real-world PIM
systems have only recently been manufactured [9, 10, 11, 12, 13]. The UPMEM
PIM system [9, 10, 14] is the first PIM architecture to become commercially
available. It consists of _UPMEM modules_ , which are standard DDR4-2400 DIMMs
with 16 PIM chips. A _PIM chip_ consists of eight small multithreaded general-
purpose in-order processors called DPUs. Each _DPU_ has exclusive access to a
64 MB DRAM bank (called _MRAM_), a 24 kB instruction memory (called _IRAM_),
and a 64 kB scratchpad memory (called _WRAM_). A common UPEM-capable system
has 20 DRAM modules with 128 DPUs and 8 GB of memory each, totaling 2,560 DPUs
with 160 GB of memory.

To program the DPUs in a UPMEM-capable system, UPMEM has developed a single-
instruction multiple-thread (SIMT) programming model. The programming model
uses a C-like interface and exposes to the programmer a series of APIs to
manage data allocation and data movement between the host CPU/DPUs and within
the memory hierarchy of the DPUs. A programmer needs to follow four main steps
to implement a given application targeting the UPMEM system. The programmer
needs to: (i) partition the computation (and input data) across the DPUs in
the system, _manually_ exposing thread-level parallelism (TLP) to the system;
(ii) distribute (copy) the appropriate input data from the CPU’s main memory
into the DPU’s memory space; (iii) launch the computation kernel that the DPUs
will execute; and (iv) gather (copy) output data from the DPUs to the CPU main
memory once the DPUs execute the kernel.

Even though UPMEM’s programming model resembles that of widely employed
architectures, such as GPUs, it requires the programmer to (i) have prior
knowledge of the underlying UPMEM hardware and (ii) manage data movement at a
fine-grained granularity _manually_. Such limitations can difficult the
adoption of PIM architectures in general-purpose systems. Therefore, our goal
in this work is to ease programmability for the UPMEM architecture, allowing a
programmer to write efficient PIM-friendly code _without_ the need to manage
hardware resources explicitly.

##  2 DaPPA: A Data-Parallel PIM Framework

To ease the programmability of PIM architectures, we propose DaPPA (data-
parallel processing-in-memory architecture), a framework that can, for a given
application, _automatically_ distribute input and gather output data, handle
memory management, and parallelize work across the DPUs. The _key idea_ behind
DaPPA is to remove the responsibility of managing hardware resources from the
programmer by providing an intuitive data-parallel pattern-based programming
interface [15, 16] that abstracts the hardware components of the UPMEM system.
Using this key idea, DaPPA transforms a data-parallel pattern-based
application code into the appropriate UPMEM-target code, including the
required APIs for data management and code partition, which can then be
compiled into a UPMEM-based binary _transparently_ from the programmer. While
generating UPMEM-target code, DaPPA implements several code optimizations to
improve end-to-end performance.

###  2.1 DaPPA Overview

Figure 1 shows an overview of our DaPPA framework. DaPPA takes as input C/C++
code, which describes the target computation using a collection of data-
parallel patterns and DaPPA’s programming interface, and generates as output
the requested computation. DaPPA consists of three main components: (i)
DaPPA’s data-parallel pattern APIs, (ii) DaPPA’s dataflow programming
interface, and (iii) DaPPA’s dynamic template-based compilation.

Figure 1: Overview of the DaPPA framework.

Data-Parallel Pattern APIs. DaPPA’s data-parallel pattern APIs ( 1 in Figure
1) are a collection of pre-defined functions that implement high-level data-
parallel pattern primitives. Each primitive allows the user to express how
data is transformed during computation. DaPPA supports five primary data-
parallel pattern primitives, including: (i) map, which applies a function $f$
to each individual input element $i$, producing unique output elements
$y_{i}=f(x_{i})$; (ii) filter, which selects input elements based on a
predicate; (iii) reduce, which reduces input elements to a scalar; (iv)
window, which _maps_ and output element as the _reduction_ of $W$
_overlapping_ input elements; (v) group, which _maps_ and output element as
the _reduction_ of $G$ _non-overlapping_ input elements. The user can combine
all five data-parallel primitives to describe complex data transformations in
an application. DaPPA is responsible for translating and parallelizing each
data-parallel primitive to efficient CPU and UPMEM code.

Dataflow Programming Interface. DaPPA exposes a dataflow-based programming
interface to the user ( 2 in Figure 1). In this programming interface, the
main component is the Pipeline class, which represents a sequence of data-
parallel patterns that will be executed on the DPUs. A given Pipeline has one
or more stages. Each stage utilizes a given data-parallel pattern primitive to
transform input operands following a user-defined computation. Stages are
executed in order, in a pipeline fashion.

Dynamic Template-Based Compilation. DaPPA uses a dynamic template-based
compilation ( 3 in Figure 1) to generate DPU code in two main steps. In the
first step, DaPPA creates a base DPU code based on a basic skeleton of a DPU
application. In the second step, DaPPA uses a series of transformations to (i)
extract the required information that will be fed to the DPU code template
from the user program; (ii) calculate the appropriate offsets used when
managing data across MRAMs and WRAMs; and (iii) divide computation between CPU
and DPUs.

Putting All Together. Using DaPPA’s data-parallel pattern APIs, data-flow
programming interface, and dynamic template-based compilation, the user can
quickly implement and deploy applications to the UPMEM system without any
knowledge of the underlying architecture. Figure 1 showcases an example of
implementing a simple vector dot product application using DaPPA. In this
example, the user defines a Pipeline with two stages: a map stage and a reduce
stage. DaPPA generates the appropriate binary for the UPMEM system, executes
the target computation in the DPUs, and copies the final output from the DPUs
to the CPU.

##  3 Key Results & Contributions

Methodology. To demonstrate DaPPA’s benefits, we implement a subset of the
workloads (i.e., vector addition, select, reduce, unique, imagine histogram
small, and gemv) presented in the UPMEM-based PrIM benchmark suite [17] using
our data-parallel pattern model. We conduct our evaluation on a UPMEM PIM
system that includes a 2-socket Intel Xeon Silver 4110 CPU at 2.10 GHz (host
CPU), standard main memory (DDR4-2400) of 128 GB, and 20 UPMEM PIM DIMMs with
160 GB PIM-capable memory and 2560 DPUs. We compare DaPPA’s performance and
programming complexity to that of the hand-tuned implementations present in
PrIM.

Key Results. First, compared to the hand-tuned PriM workloads, DaPPA improves
end-to-end performance by 2.1$\times$, on average across all six workloads
(min. 0.8$\times$, max. 10.6$\times$). DaPPA’s performance improvement is due
to code optimizations, such as parallel data transfer and workload partition
between CPU and DPUs. Second, DaPPA _significantly_ reduces programming
complexity (measured using line-of-code) on average by 94.4% (min. 92.3%, max.
96.1%). We conclude that DaPPA is an efficient framework that eases the
programmability of PIM architectures.

We make the following key contributions:

  * •

To our knowledge, this is the first work to propose a data-parallel pattern-
based framework to generate code for the UPMEM architecture _automatically_.

  * •

We propose DaPPA (data-parallel processing-in-memory architecture), a
framework that _automatically_ distributes input and gathers output data,
handles memory management, and parallelizes work across DPUs.

  * •

We equip DaPPA with a series of code optimizations that improve the
performance of workloads running on the UPMEM system.

  * •

We evaluate DaPPA using six workloads from the PrIM benchmark suite, and we
observe that DaPPA improves performance by 2.1$\times$ and reduces line-of-
code by 94.4%, on average, compared to the hand-tuned PrIM workloads.

## References

  * [1] S. Ghose _et al._ , “Processing-in-Memory: A Workload-Driven Perspective,” _IBM JRD_ , 2019.
  * [2] O. Mutlu _et al._ , “A Modern Primer on Processing in Memory,” _Emerging Computing: From Devices to Systems - Looking Beyond Moore and Von Neumann_ , 2021.
  * [3] G. F. Oliveira _et al._ , “DAMOV: A New Methodology and Benchmark Suite for Evaluating Data Movement Bottlenecks,” _IEEE Access_ , 2021.
  * [4] S. Ghose _et al._ , “The Processing-in-Memory Paradigm: Mechanisms to Enable Adoption,” in _Beyond-CMOS Technologies for Next Generation Computer Design_ , 2019.
  * [5] O. Mutlu _et al._ , “Enabling Practical Processing in and Near Memory for Data-Intensive Computing,” in _DAC_ , 2019.
  * [6] O. Mutlu, “Memory Scaling: A Systems Architecture Perspective,” in _IMW_ , 2013.
  * [7] W. H. Kautz, “Cellular Logic-in-Memory Arrays,” _IEEE TC_ , 1969.
  * [8] H. S. Stone, “A logic-in-memory computer.”
  * [9] UPMEM, “UPMEM Website,” https://www.upmem.com, 2023.
  * [10] UPMEM, “Introduction to UPMEM PIM. Processing-in-memory (PIM) on DRAM Accelerator (White Paper),” 2018.
  * [11] Y.-C. Kwon _et al._ , “25.4 A 20nm 6GB Function-In-Memory DRAM, Based on HBM2 with a 1.2 TFLOPS Programmable Computing Unit Using Bank-Level Parallelism, for Machine Learning Applications,” in _ISSCC_ , 2021.
  * [12] S. Lee _et al._ , “Hardware Architecture and Software Stack for PIM Based on Commercial DRAM Technology: Industrial Product,” in _ISCA_ , 2021.
  * [13] L. Ke _et al._ , “Near-Memory Processing in Action: Accelerating Personalized Recommendation with AxDIMM,” _IEEE Micro_ , 2021.
  * [14] J. Gómez-Luna _et al._ , “Benchmarking a New Paradigm: An Experimental Analysis of a Real Processing-in-Memory Architecture,” arXiv:2105.03814 [cs.AR], 2021.
  * [15] M. I. Cole, _Algorithmic Skeletons: Structured Management of Parallel Computation_. Pitman London, 1989.
  * [16] M. Cole, “Bringing Skeletons Out of the Closet: A Pragmatic Manifesto for Skeletal Parallel Programming,” _Parallel Computing_ , 2004.
  * [17] SAFARI Research Group, “PrIM Benchmark Suite,” https://github.com/CMU-SAFARI/prim-benchmarks.