Introduction

The thermal impact assessment is a crucial step in a geothermal project. Particularly in regards of geothermal energy planning. In this regard, analytical solutions are straightforward tools for a preliminary impact assessment. In this article we explore three analytical solutions which can help to estimate the thermal impact caused by groundwater heat pump systems. The valididy of these (semi-)analytical solutions was discussed in details by Pophillat et al. (2020) [1]. Now, let’s dive into these equations.

Problem settings

Here we study the thermal impact of a hot water injection in an aquifer. The parameters of the problem read as follow:

• we study the thermal impact in a 2D plane x, y,
• the water injection rate is Q_{inj} and temperature difference between injection groundwater background temperature is \Delta T_{inj},
• the groundwater is flowing with an angle \alpha from positives x
• the injection is located in X_{0}, Y_{0} ,
• the background groundwater seepgae velocity is v_{a}, the effective porosity of the aquifer is n and its thickness is ( b ),
• we also need to define some hydraulic/thermal parameters such as : thermal conductivity of the aquifer ( \lambda _{m} ), volumetric heat capacity of the aquifer and of water ( C _{m,w} ), longitudinal and transverse dispersivities ( \alpha _{L,T} ),
• we calculate the temperature alteration after a time t of injection.

Definition of 2D(semi-)analytical models

To define our models of intereset, we first need to import some usefull librairies.

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from scipy.special import erf,erfc,erfcinv

Radial Heat Transport model

The radial heat transport model [2] is appropriate when the background groundwater flow velocity is almost zero. In that case, the injection of the heated water creates a radial thermal disturbance according to following equation:

with: r* = \left( 2 A_{T}t \right)^{1/2}, A_{T} = \frac{Q_{inj}}{2R \pi n b}, and R = \frac{C_{m}}{nC_{w}}.

def RHT(x,y,X0,Y0,t,Qinj,DTinj,b,va,n,C_m,C_w,alpha_L,alpha_T,lambda_m):
    R = C_m/(n*C_w)
    At = (1/R)*(Qinj/(2*np.pi*n*b))
    rp = np.sqrt(2*At*t)
    
    x = x - X0                                                          
    y = y - Y0 
    
    r = np.sqrt(x**2+y**2)
    
    return (DTinj/2)*erfc((r**2-rp**2)/(2*np.sqrt((4*alpha_L/3)*(rp**3)+(lambda_m/At/C_m)*(rp**4))))

Planar Heat Transport model

The planar advective heat transport model [3] accounts for the fact that injection may induce a high local hydraulic gradient around the injection well. In such a case, the geometry of this heat source cannot be considered only with a vertical line. Instead, the source is represented as an area in the yz-plane. Accordingly, in a 2D horizontal projection in the xy-plane, the heat source corresponds to a line perpendicular to the groundwater flow direction. Please note that this model is undefined upstream the injection (for x < 0).

with \Delta T_{0} = \frac{F_{0}}{v_{a}nC_{w}Y}, F_{0} = \frac{q_{h}}{b}, q_{h} = \Delta T_{inj} C_{w} Q_{inj}, D_{x,y} = \frac{\lambda _{m}}{n C_{w}} + \alpha _{L,T} v_{a}, and  Y = \frac{Q_{inj}}{2b v_{a} n}.

def PAHT(x, y, X0, Y0, alpha, t, Qinj, DTinj, 
         b, va, n, C_m, C_w, alpha_L, alpha_T, lambda_m):
    Y = Qinj/(2*b*va*n)
    qh = DTinj*C_w*Qinj
    F0 = qh/b
    DT0 = F0/(va*n*C_w*Y)
    Dx = lambda_m/n/C_w + alpha_L*va
    Dy = lambda_m/n/C_w + alpha_T*va
    R = C_m/(n*C_w)
    
    alpha_rad = -alpha*np.pi/180
    x1 = x - X0                                                          
    y1 = y - Y0 
    
    x2 = np.cos(alpha_rad)*x1 - np.sin(alpha_rad)*y1                  
    y2 = np.sin(alpha_rad)*x1 + np.cos(alpha_rad)*y1
    
    res = DT0/4*erfc((R*x2 - va*t)/(2*np.sqrt(Dx*R*t)))*\
    (erf((y2 + Y/2)/(2*np.sqrt(Dy*x2/va))) - erf((y2 - Y/2)/(2*np.sqrt(Dy*x2/va))))

    return res

Linear Heat Transport model

The linear model [4] is appropriate for higher groundwater flow velocities. It describes heat propagation from an injection well with transient
conditions, simulated as continuous line-source, considering background flow:

with r' = \left( x^{2} + y^{2} \frac{\alpha _{L}}{\alpha _{T}}\right) ^{1/2} .

def LAHT(x, y, X0, Y0, alpha, t, Qinj, DTinj,
         b, va, n, C_m, C_w, alpha_L, alpha_T, lambda_m):
    R = C_m/n/C_w
    
    alpha_rad = -alpha*np.pi/180
    x = x - X0                                                          
    y = y - Y0 
    
    x = np.cos(alpha_rad)*x - np.sin(alpha_rad)*y                  
    y = np.sin(alpha_rad)*x + np.cos(alpha_rad)*y 
    
    rp = np.sqrt(x**2+y**2*alpha_L/alpha_T)
    res = Qinj*DTinj/(4*n*b*va*np.sqrt(np.pi*alpha_T))*\
    np.exp((x-rp)/2/alpha_L)*(1/np.sqrt(rp))*\
    erfc((rp-va*t/R)/(2*np.sqrt(va*alpha_L*t/R)))
    return res

Creation of a grid over a region of interest using numpy

To apply one of these models, we need to define a grid over a region of interest. So we need to define the extent of interest (X_{min}, X_{max} and Y_{min},Y_{max} ) and to discretize this area. We can do that using linespace and meshgrid functions of the numpy librairy :

#definition of your grid from Xmin to Xmax, and from Ymin to Ymax
Xmin = -100
Xmax = 100
xgrid_len = 200

Ymin = -100
Ymax = 100
ygrid_len = 200

#You create your grif of interest
xi = np.linspace(Xmin, Xmax, xgrid_len)                                    
yi = np.linspace(Ymin, Ymax, ygrid_len)                                    
xi, yi = np.meshgrid(xi, yi)

Plotting the results using matplotlib

Before using our models over our grid, we need to define the parameters of the problem. So let’s take the values below as example, and let’s assume that we calculate the thermal impact after one year of hot water injection:

C_m = float(2888000) #volumetric heat capacity porous media (J/ kg / K)
C_w = float(4185000) #volumetric heat capacity water (J/ kg / K)
alpha_L = float(5) #longitudinal dispersivity (m)
alpha_T = float(0.5) #transverse dispersivity (m)
lambda_m = float(2.24) #thermal conductivity (W/m/K)

alpha = float(20) # groundwater flow angle
K = float(0.001) #permeability (m/s)
b = float(10) #aquifer thickness [m]
grad_h = float(0.001) #hydraulic gradient
n = float(0.2) #effective porosity

v0 = K*grad_h #darcy velocity
va = v0/n #seepage velocity
R = C_m/(n*C_w) #retardation factor

#We also define the location of the hot water injection
X0, Y0 = 20, 20
DTinj = 10. #temperature difference between pumping and reinjection
Qinj = 0.0001 #injection rate m3/s
time = 365*24*3600 # operation time in seconds (365 days)

Then, we can calculate the thermal impact over our grid using the three models previously defined. Please note that in the case of the radial model, there is not background groundwater velocity:

deltaT_rht =  RHT(xi, yi, X0, Y0, time, Qinj, DTinj, b, va, n, C_m, C_w, alpha_L, alpha_T, lambda_m)

deltaT_paht = PAHT(xi, yi, X0, Y0, alpha, time, Qinj, DTinj, b, va, n, C_m, C_w, alpha_L, alpha_T, lambda_m)

deltaT_laht = LAHT(xi, yi, X0, Y0, alpha, time, Qinj, DTinj, b, va, n, C_m, C_w, alpha_L, alpha_T, lambda_m)

Now we can create our figure using matplotlib. The three subplots represent the thermal impact calculation using the Radial, Planar and Linear model:

#We define the titles of our subplots
titles =["Radial model (no groundwater velocity)",
         "Planar model",
         "Linear model"]

results = [deltaT_rht, deltaT_paht, deltaT_laht]
 
fig, axs = plt.subplots(1,3, figsize=(23,6))
for i in range(3):
    ax = axs[i]
    ax.set_title(titles[i])
    ax.axis('equal')
    cf= ax.contourf(xi, yi, results[i], [1,2,3,4,5,6,7,8,9,10], cmap='viridis') # analytical contour (1 K disturbance)
    ax.grid(color='black', linestyle='-', linewidth=0.1)
    ax.scatter(X0, Y0, color='red', label = "Injection point")   
    ax.set_xlim(0, 100)
    ax.set_ylim(0, 100)
    ax.legend()
 
fig.colorbar(cf, ax=axs.ravel().tolist(), label = "Temperature disturbance [K]")
 
plt.show()

Illustation in case of multiple installations

In case you want to calculate the thermal impact caused by several installations, you can iterate over your installations and calculate the total thermal impact (of course you have to assume that thermal impacts are additive which is not rigourously the case. However, this assumption can help to have an idea of the thermal stress of an higly sollicited area).

Let’s choose the linear model and calculate the thermal impact of three installations:

#inistialize an array of 0 with the good size
deltaT_m = xi - xi 

#add the impact of the first instalation
X0_1, Y0_1 = 0, 0 #location of the first installation
deltaT_m += LAHT(xi, yi, X0_1, Y0_1, 20, time, 0.0002, 10, b, 5*va, n, C_m, C_w, alpha_L, alpha_T, lambda_m)

#add the impact of the second instalation
X0_2, Y0_2 = 0, 50 #location of the second installation
deltaT_m += LAHT(xi, yi, X0_2, Y0_2, 25, time, 0.0001, 10, b , 5*va, n, C_m, C_w, alpha_L, alpha_T, lambda_m)

#add the impact of the third instalation
X0_3, Y0_3 = -50, -50 #location of the third installation
deltaT_m += LAHT(xi, yi, X0_3, Y0_3, 0, time, 0.0001, 10, b, 5*va, n, C_m, C_w, alpha_L, alpha_T, lambda_m)

#Creation of the figure
fig, ax = plt.subplots(figsize=(7,7))

ax.set_title("Cumulative impact of multiple injections")
ax.axis('equal')
cf= ax.contourf(xi, yi, deltaT_m, [1,2,3,4,5,6,7,8,9,10], cmap='viridis')
ax.grid(color='black', linestyle='-', linewidth=0.1)
ax.scatter(X0_1, Y0_1, color='red', label = "Injection point 1")
ax.scatter(X0_2, Y0_2, color='yellow', label = "Injection point 2")
ax.scatter(X0_3, Y0_3, color='orange', label = "Injection point 3")   

ax.legend()

fig.colorbar(cf, label = "Temperature disturbance [K]")

plt.show()

Exporting the results in a shapefile using shapely, pandas and geopandas

It can be convenient to export our result in a shapefile. With a shapefile, we can project the results over a geographic map in our favourite GIS (The procedure to plot these kind of result on basemaps or on interactive folium map is not covered in this article but will be presented in a future post) and cross the result with other geographical data.

In the following a procedure to transform your numpy array associated to the thermal disturbance into a geoDataframe is presented. We’ll finally export this geoDataframe as a shapefile.

Let’s first import the librairies we need to do that job:

import pandas as pd
from shapely.geometry import Point
from geopandas import GeoDataFrame

Then, we can transform our numpy array calculated in the previous section into a pandas dataFrame:

dat = np.array([xi, yi, deltaT_m]).reshape(3, -1).T
df = pd.DataFrame(dat)
df.columns = ['X', 'Y', 'DeltaT']
df = df.dropna()
df.head()

Now, we can use shapely and geopandas to transforme the dataframe into a geodataframe. Fisrt a geometry field is created in our dataframe. Secondly, the geodataframe function is used sepecifying the geometry field of the dataframe.

df['geometry'] = df.apply(lambda x: Point((float(x.X), float(x.Y))), axis=1)
geopd = GeoDataFrame(df, geometry='geometry') #Here is you geodataframe
geopd.head()

Finally, our file can be exported using the to_file() method:

geopd.to_file('my_result.shp', driver='ESRI Shapefile')

Play with the script…

This post is available as a notebook and can be run on Goolge Colab. What you need is a google account.

Bibliography

[1] Pophillat, W., Attard, G., Bayer, P., Hecht-Méndez, J., & Blum, P. (2020). Analytical solutions for predicting thermal plumes of groundwater heat pump systems. Renewable Energy, 147, 2696-2707.

[2] Guimerà, J., Ortuño, F., Ruiz, E., Delos, A., & Pérez-Paricio, A. (2007, May). Influence of ground-source heat pumps on groundwater. In Conference Proceedings: European Geothermal Congress.

[3] Hähnlein, S., Molina-Giraldo, N., Blum, P., Bayer, P., & Grathwohl, P. (2010). Ausbreitung von Kältefahnen im Grundwasser bei Erdwärmesonden. Grundwasser, 15(2), 123-133.

[4] Kinzelbach, W. (1987). Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser. Oldenbourg.