Commit 4ff1a8f6 authored by Denis Anikiev's avatar Denis Anikiev
Browse files


parent b95d587d
%% Cell type:markdown id: tags:
# Py4HIP: Python tool for Heat In Place calculations
**Py4HIP** is a **Py**thon tool **for** **H**eat **I**n **P**lace calculations.
This tool was created to perform **Heat In Place (HIP)** calculations for your region of choice.
%% Cell type:raw id: tags:
%% Cell type:markdown id: tags:
Copyright 2022
- Judith Bott (
- Laureen Benoit (
- Nora Koltzer (
- Denis Anikiev (
- Helmholtz Centre Potsdam - GFZ German Research Centre for Geosciences (GFZ,
Licensed under the EUPL-1.2 or later
%% Cell type:markdown id: tags:
Calculating the Heat In Place (HIP) for certain geological units or confined reservoirs is a standard method performed in various regions to assess the spatial variability of geothermal potential ([Nathenson, 1975](#Nathenson1975); [Muffler and Cataldi, 1978](#MufflerCataldi1978); [Garg and Combs, 2015](#GargCombs2015)).
The respective implementation in **Py4HIP** is based on a volumetric ($V$) quantification of contained energy ($H$) after [Muffler and Cataldi (1978)](#MufflerCataldi1978), where the geological unit at hand is considered as spatially variable in terms of its thickness ($M$), porosity ($\phi$), as well as density ($\rho$) and specific heat capacity ($c_p$) (for the solid rock ($m$) and brine ($brine$)).
Finally, $H$ represents (i) the excess energy stored under mean temperature ($T_{mean}$) conditions with respect to a reference temperature ($T_{ref}$, here taken as the temperature at the Earth's surface ($Z_{topo}$)); (ii) the sum of stored heat in the solid and fluid parts of the rock:
<a id='eq1'></a>
$$ H = V ((1 - \phi) \rho_{m} c_{p,m} + \phi \rho_{brine} c_{p, brine}) (T_r - T_{ref}) $$
%% Cell type:markdown id: tags:
<a id='fig1'></a>
![Fig. 1](./Structure.png)
%% Cell type:markdown id: tags:
The mean (or reservoir) temperature ($T_r$) is defined as the temperature measured or modelled at the mean depth ($Z_{mean, XY}$) of the geological unit at position XY ([Fig. 1](#fig1)).
This depth is obtained ([van Wees et al., 2012](#WeesKronimusEtAl2012)) by subtracting half of the unit’s thickness ($\frac{M_{XY}}{2}$) from the top of the layer ($Z_{min, XY}$):
$$ Z_{mean, XY} = Z_{min, XY} - \frac{M_{XY}}{2}$$
Due to its strong temperature dependency, the specific heat capacity of the solid rock components ($c_{p,m}$) at reservoir conditions is calculated using an empirical formulation ([Bär, 2012](#Baer2012); modified after [Vosteen and Schellschmidt, 2003](#VosteenSchellschmidt2003)) and considering the respective reservoir temperature ($T_r$):
$$ c_{p,m}(T_r) = 6.295 \cdot 10^{-6} \cdot T_r^{-3} - 4.99 \cdot 10^{-3} \cdot T_r^2 + 1.71 \cdot T_r + c_{p,m,T20} $$
As the density of the pore fluid $\rho_{brine}$ is controlled by salinity $S$, temperature $T$, and pressure $p$, the following two empirical formulas of [Batzle and Wang (1992)](#BatzleWang1992) are used to correct the density values at standard conditions.
First one accounts for the effects of adjusted state parameters (obtaining the density of fresh water, $\rho_{fw}$, at reservoir conditions):
$$ \rho_{fw} = 1 + 10^{-6} \cdot ( -80 T - 3.3 T^2 + 0.00175 T^3 + 489p - 2Tp + 0.016 T^2 p - 1.3 \cdot 10^{-5} T^3 p - 0.333 p^2 - 0.002 T p^2)$$
where $p$ is calculated as the hydrostatic pressure.
Finally, **Py4HIP** accounts for the salinity effect ($\rho_{brine}$):
$$ \rho_{brine} = \rho_{fw} + S \cdot (0.668 + 0.44 S + 10^{-6} \cdot (300p - 2400pS + T \cdot (80 + 3T - 3300S - 13p + 47pS)))$$
The volume of rock ($V$) required to calculate the stored thermal energy ($H$; [eq. 1](#eq1)) is equal to the product of the total thickness of the reservoir unit ($M$) at the respective location ($XY$) and one square meter.
Hence, the energy values provided by **Py4HIP** as ASCII lists (and optionally as map representations, see below) correspond to the stored energy in $J/m^2$.
%% Cell type:markdown id: tags:
## Getting started
Before you can run the Heat-In-Place calculations, you need to undertake the following steps:
- Import relevant python packages
- Define your folder locations (paths)
- Ensure your files are in correct format
- Define steady parameters
%% Cell type:markdown id: tags:
### Import python packages
Before you can run the calculation you need to import the following python packages into your python environment:
%% Cell type:code id: tags:
``` python
from IPython.display import display, HTML
import numpy as np
import pandas as pd
from scipy.interpolate import griddata
from pathlib import Path
from ipywidgets import widgets
import matplotlib.pyplot as plt
import as cm
from mpl_toolkits.axes_grid1 import make_axes_locatable
%% Cell type:markdown id: tags:
### Define paths
Provided example is organized as follows:
- folder `input_path` contains test data files
- `input_file` represents the file name to be loaded
- `output_path` is the folder to store the results
- `output_path_fig` is the folder to store the resulting figures
Both `output_path` and `output_path_fig` are created after first execution of the script.
By default the script does not save the figures, in case you want to generate figures as well, enable it in the checkbox further down in this notebook.
User is free to adjust this script to its own needs.
%% Cell type:code id: tags:
``` python
input_path = './data/' # this folder contains test files
input_file = 'test_dataset.csv' # name of the file to load
output_path = './results/' # this folder will be created and will contain results
output_path_fig = output_path + 'figures/' # this folder will be created and will contain figures
# Create paths:
Path(output_path_fig).mkdir(parents=True, exist_ok=True) # this creates both output paths at once
%% Cell type:markdown id: tags:
### File setup
In order to run the caculation smooothly, your input files should follow a specific structure that goes well with the formulas used in this tool.
The desirable file structure of your input data table must contain the columns shown here:
| <div style="width:50px">x</div> | <div style="width:50px">y</div> | <div style="width:50px">$d_{mean}$ </div>| <div style="width:50px">$M_{total}$ </div> | <div style="width:50px">$T_{mean}$ </div> |
| 0 | 0 | 0 | 161.0944798 | 12.69794223 |
| 1004.016064 | 0 | 0 | 160.4305277 | 12.67989258 |
| 2008.032129 | 0 | 0 | 159.1173317 | 12.64417627 |
| ... | ... | ... | ... | ... |
%% Cell type:markdown id: tags:
### Define parameters for calculation
Define several steady parameters for your calculation:
| <div style="width:100px">Variable</div> | <div style="width:200px">Description</div> | <div style="width:100px">Unit</div>| <div style="width:100px">SI-unit</div> |
| $g$ | standard gravity | $\frac{m}{s^2}$ | $\frac{m}{s^2}$ |
| $atm$ | atmospheric pressure | $Pa$ | $\frac{kg}{ms^2}$ |
| $T_{ref}$ | temperature at topography,model reference temperature |$°C$ | $K$ |
| $\phi$ | porosity | % | - |
| $\rho_{brine\ 0}$ | density of brine (empirical) | $\frac{kg}{m^3}$ | $\frac{kg}{m^3}$ |
| $Cp_{brine}$ | speciifc heat capacity of brine | $\frac{J}{kgK}$ | $\frac{m^2}{Ks^2}$ |
| $S$ | Salinity of brine | $\frac{kg}{m^3}$ | $\frac{kg}{m^3}$ |
| $\rho_{solid}$ | density of solid | $\frac{kg}{m^3}$ | $\frac{kg}{m^3}$ |
| $Cp_{solid\ 0}$ | speciifc heat capacity of solid (empirical) | $\frac{J}{kgK}$ | $\frac{m^2}{Ks^2}$ |
These steady parameters need to be defined here:
%% Cell type:code id: tags:
``` python
g = 9.81 #[m/s^2]
atm = 101325 #[Pa; kg/(m*s^2)]
Tref = 8 #[°C; K]
Phi = 0.14 #[%]
Rho_brine_0 = 1040 #[kg/m^3]
Cpm_brine = 3.925*1000 #[J/(kg*K); m^2/(K*s^2)]
S = 60 #[kg/m^3]
Rho_solid = 2680 #[kg/m^3]
Cp_solid_0 = 810 #[J/(kg*K); m^2/(K*s^2)]
%% Cell type:markdown id: tags:
## Perform calculation
Now, with all preparations done, you should be able to smoothly run the calculation by running the following cells.
<div class="alert alert-block alert-warning">
<b>Note:</b> There is an option to generate figures in addition to the output implemented in the script. You can decide whether executing or skipping the option by ticking or unticking the checker box 'Want to generate figures?' below.</div>
%% Cell type:markdown id: tags:
### Load files
Firstly, load your file into the script.
%% Cell type:markdown id: tags:
You have to assign the correct header manually in `line 6`, as this is dependant on your personal input table. Here, we assigned the header of the test-dataset.
Use `line 3` (by uncommenting it) if you want to load a `.txt` or `.dat` file and assign the correct delimiter (test-dataset uses space). Use `line 5` if you want to load a `.csv` file. (If your input file is an excel table, save it as Comma Seperated Values file `.csv` first.)
You should see your input data displayed correctly below, if everything has loaded without errors. Make sure to check for correctness here!
%% Cell type:code id: tags:
``` python
# Load data files
# line 3 is to load data files in .txt or .dat format:
# df_test = pd.read_csv(input_path + input_file, comment='#', delimiter = ' ', header=None)
# line 5 is to load data file is .csv (comma-separated format):
df_test = pd.read_csv(input_path + input_file, comment='#', delimiter = ',')
# Assign headers (depends on the input table)
df_test.columns = ['x', 'y', 'dmean', 'Mtotal', 'Tmean', 'Ztop', 'Zmean']
# Display loaded data
%% Output
%% Cell type:markdown id: tags:
### Heat In Place calculation
This is the part where the data finally get processed.
**Don't change anything in the code cell below! All input and adjustments have been defined in the previous steps!**
%% Cell type:code id: tags:
``` python
# Heat in Place Calculation
df_test['P_hydr'] = Rho_brine_0 * g * df_test['dmean'] + atm
df_test['Cpm_solid']=(0.000006295*df_test['Tmean']**3)-(0.00499*df_test['Tmean']**2)+(1.71*df_test['Tmean'])+ Cp_solid_0
df_test['H']=df_test['Mtotal']*(((1-Phi)* Rho_solid *df_test['Cpm_solid'])+(Phi*df_test['Rho_brine']*Cpm_brine))*(df_test['Tmean']- Tref)
%% Output
0 2.062227e+09
1 2.045876e+09
2 2.013719e+09
3 1.966807e+09
4 1.906610e+09
52495 2.780037e+07
52496 3.017139e+07
52497 3.215381e+07
52498 3.358419e+07
52499 3.433459e+07
Name: H, Length: 52500, dtype: float64
%% Cell type:markdown id: tags:
## Save results
Lastly, save your freshly generated dataset to your output folder, which you have chosen in the beginning.
%% Cell type:markdown id: tags:
<div class="alert alert-block alert-success"><b>Save dataset:</b> Save the data in required format by uncommenting the respective line of which format you'd like to save your data in.
One also needs to assign the correct header (which is the same as in section Load files above) AND the new, just generated columns.</div>
%% Cell type:code id: tags:
``` python
# Save generated dataset to the desired output folder
df_test.to_csv(output_path + 'HIP_test.txt', sep=' ',
columns=['x', 'y', 'dmean', 'Mtotal', 'Tmean', 'Ztop', 'Zmean', 'H'], header = None, index = None)
# df_test.to_csv(output_path + 'HIP_test.dat', sep=' ',
# columns=['x', 'y', 'dmean', 'Mtotal', 'Tmean', 'Ztop', 'Zmean', 'H'], header = None, index = None)
# df_test.to_excel(output_path + 'HIP_test.xlsx', sheet_name='HIP_results',
# columns=['x', 'y', 'dmean', 'Mtotal', 'Tmean', 'Ztop', 'Zmean', 'H'], header=True, index = None)
%% Cell type:markdown id: tags:
## Plot figures (optional)
%% Cell type:markdown id: tags:
<div class="alert alert-block alert-warning">
<b>This is optional:</b> You can decide whether executing or skipping the option by ticking or unticking the checker box 'Want to generate figures?'.</div>
%% Cell type:code id: tags:
``` python
checker1 = widgets.Checkbox(
description='Want to generate figures?',
#hide this
%% Output
%% Cell type:markdown id: tags:
If you wish to have some nice figures of your just generated Heat In Place dataset, after ticking the checkbox above rerun the code cell below by clicking on its triangle (below left) or pressing <kbd>SHIFT</kbd> + <kbd>ENTER</kbd> when the code cell below has a blue or green frame.
Additionally, here is a short explanation of code lines below in which you **might** need to adjust something for your personal data set:
- `line 2` add .reshape(x, y) for reshaping the array into your model extend/format
- `line 5` replace 0.1 with your default thickness value for non-deposited lithology
- `line 18/19` replace 500 by your chosen resolution for interpolation
- `line 21/22` replace linear by your chosen type of interpolation
- `line 25` set the title of your figure by replacing `'Example'`
<div class="alert alert-block alert-success"><b>Save figure:</b> Set the name of the figure file and desired DPI in the last line.</div>
%% Cell type:code id: tags:
``` python
if checker1.value == True:
HIP_array = df_test['H'].values # 210, 250 is this test model extend .reshape(210, 250)
HIP_array = HIP_array / 10**9 # for pic in [PJ]
Mtotal_array = df_test['Mtotal'].values
HIP_array[Mtotal_array <= 0.1] = np.nan # replace not deposited lithology points as NAN values only for pics
# minimum and maximum extent of model
xmin = np.min(df_test['x'])
xmax = np.max(df_test['x'])
ymin = np.min(df_test['y'])
ymax = np.max(df_test['y'])
# variables for interpolation (question: What?)
x = df_test['x']
y = df_test['y']
z = HIP_array
z_M = Mtotal_array
# define axes for interpolation grid (question: Where?):
xi = np.linspace(xmin, xmax, 500) # 500 as chosen resolution
yi = np.linspace(ymin, ymax, 500)
# perform interpolation: ( )
zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='linear')
zi_M = griddata((x, y), z_M, (xi[None,:], yi[:,None]), method='linear')
fig_1, (ax2, ax1) = plt.subplots(1,2, figsize = (10, 5))
fig_1.suptitle('Example', fontsize = 'x-large') # Set the tite of your figure here!
# define color for nan values in chosen color map (
current_cmap = plt.get_cmap('Spectral_r').copy() # ,lut=4 integer giving the number of entries desired in the lookup table
im1 = ax1.imshow(zi, origin = 'lower', extent = [xmin, xmax, ymin, ymax], cmap=current_cmap)
im2 = ax2.imshow(zi_M, origin = 'lower', extent = [xmin, xmax, ymin, ymax], cmap='gist_yarg')
divider = make_axes_locatable(ax1)
cax1 = divider.append_axes("bottom", size="5%", pad=0.42)
plt.colorbar(im1, cax=cax1, orientation='horizontal')
cax1.set_xlabel('heat in place [PJ]')
divider = make_axes_locatable(ax2)
cax2 = divider.append_axes("bottom", size="5%", pad=0.42)
plt.colorbar(im2, cax=cax2, orientation='horizontal')
cax2.set_xlabel('layer thickness [m]')
plt.savefig(output_path_fig + 'HIP_test_dataset.png', dpi=300)
%% Output
%% Cell type:markdown id: tags:
## References
<a id='Baer2012'></a>
Bär, K. M. (2012). Untersuchung der tiefengeothermischen Potenziale von Hessen. (Ph.D. Thesis). TU Darmstadt, Darmstadt.
<a id='BatzleWang1992'></a>
Batzle, M., & Wang, Z. (1992). Seismic properties of pore fluids. Geophysics, 57(11), 1396-1408.
<a id='GargCombs2015'></a>
Garg, S. K., & Combs, J. (2015). A reformulation of USGS volumetric “heat in place” resource estimation method. Geothermics, 55, 150-158.
<a id='MufflerCataldi1978'></a>
Muffler, P., & Cataldi, R. (1978). Methods for regional assessment of geothermal resources. Geothermics, 7(2-4), 53-89.
<a id='Nathenson1975'></a>
Nathenson, M. (1975). Physical factors determining the fraction of stored energy recoverable from hydrothermal convection systems and conduction-dominated areas (No. USGS-OFR-75-525). Geological Survey, Menlo Park, Calif.(USA).
<a id='WeesKronimusEtAl2012'></a>
Van Wees, J.-D., Kronimus, A., Van Putten, M., Pluymaekers, M., Mijnlieff, H., Van Hooff, P., Obdam, A., Kramers, L. (2012). Geothermal aquifer performance
assessment for direct heat production–Methodology and application to Rotliegend aquifers. Netherlands Journal of Geosciences, 91(4), 651-665.
<a id='VosteenSchellschmidt2003'></a>
Vosteen, H.-D., & Schellschmidt, R. (2003). Influence of temperature on thermal conductivity, thermal capacity and thermal diffusivity for different types of rock. Physics and Chemistry of the Earth, Parts A/B/C, 28(9-11), 499-509.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment