Commit 3af116de by Denis Anikiev

### fixes im math

parent 5895b255
 %% 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. *Copyright (Bott et al., 2022). 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: $$H = V ((1 - \phi) \rho_{m} c_{p,m} + \phi \rho_{brine} c_{p, brine}) (T_r - T_{ref})$$
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)))$$ %% 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 matplotlib.cm 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: |
x
|
y
|
$d_{mean}$
|
$M_{total}$
|
$T_{mean}$
| |:------------|:-----|:-----------|:-------------|:------------| | 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: |
Variable
|
Description
|
Unit
|
SI-unit
| |:----------------------------------------|:-------------------------------------------|:-----------------------------------|:---------------------------------------| | $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.
Be aware: There is one option implemented in the skript which is to additionally generate figures to your calculation output. You can decide whether executing or skipping the option by ticking or unticking the checker box 'Want to generate figures?'.
%% 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 display(df_test)  %% 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['Rho_water']=1+0.000001*((-80*df_test['Tmean'])-(3.3*df_test['Tmean']**2)+(0.00175*df_test['Tmean']**3)+(489*(df_test['P_hydr']*10**-6))-(2*df_test['Tmean']*(df_test['P_hydr']*10**-6))+(0.016*df_test['Tmean']**2*(df_test['P_hydr']*10**-6))-(1.3*10**-5*df_test['Tmean']**3*(df_test['P_hydr']*10**-6))-(0.333*(df_test['P_hydr']*10**-6)**2)-(0.002*df_test['Tmean']*(df_test['P_hydr']*10**-6)**2)) df_test['Rho_brine']=df_test['Rho_water']+(S*(0.668+(0.44*S)+(1*10**-6*((300*df_test['P_hydr']*10**-6)-(2400*df_test['P_hydr']*10**-6*S)+(df_test['Tmean']*(80+(3*df_test['Tmean'])-(3300*S)-(13*df_test['P_hydr']*10**-6)+(47*df_test['P_hydr']*10**-6*S))))))) 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) print(df_test['H'])  %% 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:
Save dataset: 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.
%% 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:
This is optional: You can decide whether executing or skipping the option by ticking or unticking the checker box 'Want to generate figures?'.
%% Cell type:code id: tags:  python checker1 = widgets.Checkbox( value=False, description='Want to generate figures?', disabled=False, indent=False ) display(checker1) #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 SHIFT + ENTER 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'
Save figure: Set the name of the figure file and desired DPI in the last line.
%% 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: ( https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html ) 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 (https://matplotlib.org/stable/tutorials/colors/colormaps.html) current_cmap = plt.get_cmap('Spectral_r').copy() # ,lut=4 integer giving the number of entries desired in the lookup table current_cmap.set_bad(color='white') im1 = ax1.imshow(zi, origin = 'lower', extent = [xmin, xmax, ymin, ymax], cmap=current_cmap) ax1.autoscale(tight=True) im2 = ax2.imshow(zi_M, origin = 'lower', extent = [xmin, xmax, ymin, ymax], cmap='gist_yarg') ax2.autoscale(tight=True) 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.tight_layout() plt.savefig(output_path_fig + 'HIP_test_dataset.png', dpi=300)  %% Output %% Cell type:markdown id: tags: ## References Bär, K. M. (2012). Untersuchung der tiefengeothermischen Potenziale von Hessen. (Ph.D. Thesis). TU Darmstadt, Darmstadt. https://tuprints.ulb.tu-darmstadt.de/3067/ Batzle, M., & Wang, Z. (1992). Seismic properties of pore fluids. Geophysics, 57(11), 1396-1408. https://www.doi.org/10.1190/1.1443207 Garg, S. K., & Combs, J. (2015). A reformulation of USGS volumetric “heat in place” resource estimation method. Geothermics, 55, 150-158. https://www.doi.org/10.1016/j.geothermics.2015.02.004 Muffler, P., & Cataldi, R. (1978). Methods for regional assessment of geothermal resources. Geothermics, 7(2-4), 53-89. https://www.doi.org/10.1016/0375-6505(78)90002-0 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). https://www.doi.org/10.2172/7348122 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. https://doi.org/10.1017/S0016774600000433 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. https://doi.org/10.1016/S1474-7065(03)00069-X ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!