1"""Functions to calculate salt cavern volumes and storage capacities.
2
3References
4----------
5.. [#Jannel22] Jannel, H. and Torquet, M. (2022). Conceptual design of salt
6 cavern and porous media underground storage site. Hystories deliverable
7 D7.1-1. Hystories. Available at:
8 https://hystories.eu/wp-content/uploads/2022/05/Hystories_D7.1-1-Conceptual-design-of-salt-cavern-and-porous-media-underground-storage-site.pdf
9 (Accessed: 9 October 2023).
10.. [#Williams22] Williams, J. D. O., Williamson, J. P., Parkes, D., Evans, D.
11 J., Kirk, K. L., Sunny, N., Hough, E., Vosper, H., and Akhurst, M. C.
12 (2022). ‘Does the United Kingdom have sufficient geological storage
13 capacity to support a hydrogen economy? Estimating the salt cavern storage
14 potential of bedded halite formations’, Journal of Energy Storage, 53, p.
15 105109. https://doi.org/10.1016/j.est.2022.105109.
16.. [#Caglayan20] Caglayan, D. G., Weber, N., Heinrichs, H. U., Linßen, J.,
17 Robinius, M., Kukla, P. A., and Stolten, D. (2020). ‘Technical potential
18 of salt caverns for hydrogen storage in Europe’, International Journal of
19 Hydrogen Energy, 45(11), pp. 6793–6805.
20 https://doi.org/10.1016/j.ijhydene.2019.12.161.
21.. [#Allsop23] Allsop, C. et al. (2023) ‘Utilizing publicly available datasets
22 for identifying offshore salt strata and developing salt caverns for
23 hydrogen storage’, Geological Society, London, Special Publications,
24 528(1), pp. 139–169. Available at: https://doi.org/10.1144/SP528-2022-82.
25.. [#English23] English, J. M., English, K. L., Dunphy, R. B., Blake, S.,
26 Walsh, J., Raine, R., Vafeas, N. A., and Salgado, P. R. (2023). ‘An
27 Overview of Deep Geothermal Energy and Its Potential on the Island of
28 Ireland’, First Break, 41(2), pp. 33–43.
29 https://doi.org/10.3997/1365-2397.fb2023009.
30.. [#Nayar23] Nayar, K.G., Sharqawy, M.H. and Lienhard V, J.H. (2023)
31 Thermophysical properties of seawater. Available at:
32 https://web.mit.edu/seawater (Accessed: 10 April 2024).
33.. [#Nayar16] Nayar, K.G. et al. (2016) ‘Thermophysical properties of seawater:
34 A review and new correlations that include pressure dependence’,
35 Desalination, 390, pp. 1–24. Available at:
36 https://doi.org/10.1016/j.desal.2016.02.024.
37.. [#Bell14] Bell, I. H., Wronski, J., Quoilin, S., and Lemort, V. (2014).
38 ‘Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the
39 Open-Source Thermophysical Property Library CoolProp’, Industrial &
40 Engineering Chemistry Research, 53(6), pp. 2498–2508.
41 https://doi.org/10.1021/ie4033999.
42.. [#PyFluids] Portyanikhin, V. (2023). ‘portyanikhin/PyFluids’. Available at:
43 https://github.com/portyanikhin/PyFluids (Accessed: 1 January 2024).
44.. [#CoolProp] Bell, I. H. and the CoolProp Team (n.d.). Hydrogen - CoolProp
45 documentation. Available at:
46 http://www.coolprop.org/fluid_properties/fluids/Hydrogen.html
47 (Accessed: 10 December 2023).
48.. [#HydrogenTools] Hydrogen Tools (no date) Lower and Higher Heating Values
49 of Fuels. Available at:
50 https://h2tools.org/hyarc/calculator-tools/lower-and-higher-heating-values-fuels
51 (Accessed: 25 April 2024).
52"""
53
54import numpy as np
55import pandas as pd
56from pyfluids import Fluid, FluidsList, Input
57
58from h2ss import functions as fns
59
60
[docs]
61def cavern_volume(height, diameter=fns.CAVERN_DIAMETER, theta=20):
62 """Calculate the cavern volume.
63
64 Parameters
65 ----------
66 height : float
67 Cavern height [m]
68 diameter : float
69 Cavern diameter [m]
70 theta : float
71 Cavern roof angle [°]
72
73 Returns
74 -------
75 float
76 Ideal cavern volume [m³]
77
78 Notes
79 -----
80 The cavern's ideal shape is a cylinder of height :math:`h_{cylinder}` [m]
81 with two conical ends each with a height of :math:`h_{cone}` [m] based on
82 [#Jannel22]_. Since the volume of a cylinder is defined as
83 :math:`\\pi r^2 h_{cylinder}` and the volume of a cone is defined as
84 :math:`\\frac{\\pi r^2 h_{cone}}{3}` (i.e. one third of a cylinder's
85 volume), the bulk cavern volume, :math:`V_{bulk}` [m³] can therefore be
86 deduced by calculating the volume of a cylinder of the cavern height,
87 :math:`h` [m] and deducting 4/3 of the volume of a cone of height
88 :math:`h_{cone}`. :math:`r` [m] is the radius of the cavern and
89 :math:`h_{cone}` is therefore equivalent to :math:`r \\tan(\\theta)`,
90 where :math:`\\theta` [°] is the cavern roof angle.
91
92 .. math::
93 V_{bulk} = \\pi \\, r^2 \\, h
94 - \\frac{4}{3} \\, \\pi \\, r^2 \\, h_{cone}
95 .. math::
96 V_{bulk} = \\pi \\, r^2
97 \\left(h - \\frac{4}{3} \\, h_{cone}\\right)
98 .. math::
99 V_{bulk} = \\pi \\, r^2 \\left(h - \\frac{4}{3} \\,
100 r \\, \\tan(\\theta)\\right)
101 """
102 r = diameter / 2
103 v_cavern = (
104 np.pi * np.square(r) * (height - 4 / 3 * r * np.tan(np.deg2rad(theta)))
105 )
106 return v_cavern
107
108
[docs]
109def corrected_cavern_volume(v_cavern, scf=0.7):
110 """Apply correction factors to the cavern volume.
111
112 Parameters
113 ----------
114 v_cavern : float
115 Ideal cavern volume [m³]
116 scf : float
117 Shape correction factor
118
119 Returns
120 -------
121 float
122 Corrected cavern volume [m³]
123
124 Notes
125 -----
126 The corrected cavern volume, :math:`V_{cavern}` [m³] is defined as
127
128 .. math::
129 V_{cavern} = V_{bulk} \\times SCF
130
131 where :math:`V_{bulk}` [m³] is the bulk cavern volume and :math:`SCF` is
132 the shape correction factor (deviation of the cavern's shape due to
133 geological differences). A :math:`SCF` of 0.7 is used as the default, as
134 defined in [#Jannel22]_, [#Williams22]_, [#Caglayan20]_, and [#Allsop23]_.
135 """
136 return v_cavern * scf
137
138
[docs]
139def temperature_cavern_mid_point(height, depth_top, t_0=10, delta_t=37.5):
140 """Cavern mid-point temperature.
141
142 Parameters
143 ----------
144 height : float
145 Cavern height [m]
146 depth_top : float
147 Cavern top depth [m]
148 t_0 : float
149 Mean annual seabed surface temperature [°C]
150 delta_t : float
151 Geothermal gradient; change in temperature with depth [°C km⁻¹]
152
153 Returns
154 -------
155 float
156 Mid-point temperature [K]
157
158 Notes
159 -----
160 See [#Williams22]_, Eqn. (2).
161
162 .. math::
163 T_{midpoint} = T_0 + \\Delta_T \\,
164 \\frac{z + 0.5 \\, h}{1,000} + 273.15
165
166 where :math:`T_{midpoint}` is the cavern mid-point temperature [K],
167 :math:`T_0` is the mean annual surface temperature [°C], :math:`\\Delta_T`
168 is the change in temperature with depth, i.e. geothermal gradient
169 [°C km⁻¹], :math:`z` is the cavern top depth [m], and :math:`h` is the
170 cavern height [m].
171
172 A :math:`\\Delta_T` of 37.5 °C km⁻¹ is used based on the geothermal
173 gradient of Kish Basin wells in the Mercia Mudstone Group reported in
174 [#English23]_.
175 Note that the cavern top depth is used instead of the casing shoe depth,
176 as the casing shoe is 30 m above the cavern top depth in this study,
177 rather than at the cavern top as modelled in [#Williams22]_.
178 """
179 return t_0 + delta_t * (depth_top + height / 2) / 1000 + 273.15
180
181
[docs]
182def pressure_operating(
183 depth_water,
184 thickness_overburden,
185 thickness_salt=50,
186 rho_overburden=2400,
187 rho_salt=2200,
188 rho_water=1027,
189 minf=0.3,
190 maxf=0.8,
191):
192 """Cavern operating pressures.
193
194 Parameters
195 ----------
196 depth_water : float
197 Sea water depth [m]
198 thickness_overburden : float
199 Overburden thickness / halite top depth [m]
200 thickness_salt : float
201 Thickness of the salt above the casing shoe [m]
202 rho_overburden : float
203 Density of the overburden [kg m⁻³]
204 rho_salt : float
205 Density of the halite [kg m⁻³]
206 rho_water : float
207 Density of the sea water [kg m⁻³]
208 minf : float
209 Factor of lithostatic pressure for the minimum operating pressure
210 maxf : float
211 Factor of lithostatic pressure for the maximum operating pressure
212
213 Returns
214 -------
215 tuple[float, float]
216 Min and max operating pressures [Pa]
217
218 Notes
219 -----
220 .. math::
221 p_{casing} = (\\rho_{water} \\, t_{water} + \\rho_{overburden}
222 \\, t_{overburden} + \\rho_{salt} \\, t_{salt}) \\, g
223 .. math::
224 p_{H_2min} = 0.3 \\, p_{casing}
225 .. math::
226 p_{H_2max} = 0.8 \\, p_{casing}
227
228 :math:`p_{casing}` is the lithostatic pressure at the casing shoe [Pa].
229 The thickness of the overburden, :math:`t_{overburden}` [m] is the same
230 as the depth to top of salt. :math:`t_{salt}` [m] is the thickness of
231 the salt above the casing shoe. :math:`t_{water}` [m] is the sea water
232 depth. :math:`\\rho_{water}`, :math:`\\rho_{overburden}`, and
233 :math:`\\rho_{salt}` are the densities of the sea water, overburden, and
234 salt, respectively [kg m⁻³]. :math:`g` is the acceleration due to gravity
235 [m s⁻²].
236
237 :math:`p_{H_2min}` and :math:`p_{H_2max}` are the minimum and
238 maximum cavern operating pressures, respectively [Pa].
239
240 See [#Williams22]_, Eqn. (3) and (4). This function has been modified to
241 include sea water depth and density to suit offshore conditions.
242 The sea water density is assumed as the mean value at a temperature of 10
243 °C, a salinity of 35 g kg⁻¹, and a surface pressure of 1 atm [#Nayar23]_
244 [#Nayar16]_.
245 """
246 p_casing = (
247 rho_overburden * thickness_overburden
248 + rho_salt * thickness_salt
249 + rho_water * depth_water
250 ) * 9.81
251 p_operating_min = minf * p_casing
252 p_operating_max = maxf * p_casing
253 return p_operating_min, p_operating_max
254
255
[docs]
256def density_hydrogen_gas(p_operating_min, p_operating_max, t_mid_point):
257 """Density of hydrogen at cavern conditions.
258
259 Parameters
260 ----------
261 p_operating_min : float
262 Minimum operating pressure [Pa]
263 p_operating_max : float
264 Maximum operating pressure [Pa]
265 t_mid_point : float
266 Mid-point temperature [K]
267
268 Returns
269 -------
270 tuple[float, float]
271 Min and max hydrogen gas densities [kg m⁻³]
272
273 Notes
274 -----
275 This function uses the CoolProp [#Bell14]_ wrapper called PyFluids
276 [#PyFluids]_. See [#Williams22]_, Section 3.4.2 and also the CoolProp
277 documentation for useful information on hydrogen [#CoolProp]_.
278 The ``pyproject.toml`` configuration file has been set such that the
279 default units used by PyFluids are SI units. PyFluids can also be used to
280 derive the compressibility factor.
281
282 Based on [#Caglayan20]_, Eqn. (3), the density can be approximated using
283 the following equation.
284
285 .. math::
286 \\rho = \\frac{p \\, M}{Z \\, R \\, T}
287
288 where :math:`M` is the molar mass of hydrogen gas [kg mol⁻¹], :math:`R` is
289 the universal gas constant [J K⁻¹ mol⁻¹], :math:`p` is the pressure [Pa],
290 :math:`T` is the temperature [K], :math:`Z` is the compressibility factor,
291 and :math:`\\rho` is the hydrogen gas density [kg m⁻³].
292 """
293 rho_h2 = []
294 for p_min, p_max, t in zip(p_operating_min, p_operating_max, t_mid_point):
295 h2_min = Fluid(FluidsList.Hydrogen).with_state(
296 Input.pressure(p_min), Input.temperature(t)
297 )
298 h2_max = Fluid(FluidsList.Hydrogen).with_state(
299 Input.pressure(p_max), Input.temperature(t)
300 )
301 rho_h2.append((h2_min.density, h2_max.density))
302 rho_h2 = pd.DataFrame(rho_h2)
303 return rho_h2[0], rho_h2[1]
304
305
[docs]
306def mass_hydrogen_working(rho_h2_min, rho_h2_max, v_cavern):
307 """The working mass of hydrogen.
308
309 Parameters
310 ----------
311 rho_h2_min : float
312 Minimum hydrogen density [kg m⁻³]
313 rho_h2_max : float
314 Maximum hydrogen density [kg m⁻³]
315 v_cavern : float
316 Cavern volume [m³]
317
318 Returns
319 -------
320 tuple[float, float, float]
321 Working mass and the minimum and maximum operating mass of hydrogen
322 [kg]
323
324 Notes
325 -----
326 See [#Williams22]_, Eqn. (5) and (6).
327 The mass of hydrogen at the minimum operating pressure represents the
328 cushion gas requirement.
329
330 .. math::
331 m_{min} = \\rho_{H_2min} \\, V_{cavern}
332 .. math::
333 m_{max} = \\rho_{H_2max} \\, V_{cavern}
334 .. math::
335 m_{working} = m_{max} - m_{min}
336
337 The working mass of hydrogen, :math:`m_{working}` [kg] is the difference
338 between the stored mass of hydrogen at maximum, :math:`m_{max}` and
339 minimum, :math:`m_{min}` operating pressures [kg], which were derived
340 using the minimum, :math:`\\rho_{H_2min}` and maximum,
341 :math:`\\rho_{H_2max}` hydrogen densities [kg m⁻³], respectively, and the
342 cavern volume, :math:`V_{cavern}` [m³].
343 """
344 m_min_operating = rho_h2_min * v_cavern
345 m_max_operating = rho_h2_max * v_cavern
346 m_working = m_max_operating - m_min_operating
347 return m_working, m_min_operating, m_max_operating
348
349
[docs]
350def energy_storage_capacity(m_working):
351 """Cavern energy storage capacity.
352
353 Parameters
354 ----------
355 m_working : float
356 Working mass of hydrogen [kg]
357
358 Returns
359 -------
360 float
361 Energy storage capacity of cavern [GWh]
362
363 Notes
364 -----
365 See [#Williams22]_, Eqn. (7).
366 See also [#HydrogenTools]_ for the heating values.
367 The energy storage capacity of the cavern, :math:`E_{cavern}` [GWh] is a
368 function of the working mass of hydrogen, :math:`m_{working}` [kg] and the
369 lower heating value of hydrogen, :math:`LHV` [MJ kg⁻¹].
370
371 .. math::
372 E_{cavern} = \\frac{m_{working} \\times LHV}{3,600,000}
373 """
374 return m_working * 119.96 / 3.6e6