Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EnthalpyConverter.cc
Go to the documentation of this file.
1// Copyright (C) 2009-2017, 2021, 2023, 2024, 2025 Andreas Aschwanden, Ed Bueler and Constantine Khroulev
2//
3// This file is part of PISM.
4//
5// PISM is free software; you can redistribute it and/or modify it under the
6// terms of the GNU General Public License as published by the Free Software
7// Foundation; either version 3 of the License, or (at your option) any later
8// version.
9//
10// PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13// details.
14//
15// You should have received a copy of the GNU General Public License
16// along with PISM; if not, write to the Free Software
17// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19#include "pism/util/EnthalpyConverter.hh"
20#include "pism/util/ConfigInterface.hh"
21
22#include "pism/util/error_handling.hh"
23#include "pism/util/pism_utilities.hh"
24
25namespace pism {
26
27/*! @class EnthalpyConverter
28
29 Maps from @f$ (H,p) @f$ to @f$ (T,\omega,p) @f$ and back.
30
31 Requirements:
32
33 1. A converter has to implement an invertible map @f$ (H,p) \to (T,
34 \omega, p) @f$ *and* its inverse. Both the forward map and the
35 inverse need to be defined for all permissible values of @f$
36 (H,p) @f$ and @f$ (T, \omega, p) @f$, respectively.
37
38 2. A converter has to be consistent with laws and parameterizations
39 used elsewhere in the model. This includes models coupled to
40 PISM.
41
42 3. For a fixed volume of liquid (or solid) water and given two
43 energy states @f$ (H_1, p_1) @f$ and @f$ (H_2, p_2) @f$ , let @f$
44 \Delta U_H @f$ be the difference in internal energy of this
45 volume between the two states *computed using enthalpy*. We require
46 that @f$ \Delta U_T = \Delta U_H @f$ , where @f$ \Delta U_T @f$
47 is the difference in internal energy *computed using corresponding*
48 @f$ (T_1, \omega_1, p_1) @f$ *and* @f$ (T_2, \omega_2, p_2) @f$.
49
50 4. We assume that ice and water are incompressible, so a change in
51 pressure does no work, and @f$ \Diff{H}{p} = 0 @f$. In addition
52 to this, for cold ice and liquid water @f$ \Diff{T}{p} = 0 @f$.
53*/
54
56 m_p_air = config.get_number("surface.pressure"); // Pa
57 m_g = config.get_number("constants.standard_gravity"); // m s-2
58 m_beta = config.get_number("constants.ice.beta_Clausius_Clapeyron"); // K Pa-1
59 m_rho_i = config.get_number("constants.ice.density"); // kg m-3
60 m_c_i = config.get_number("constants.ice.specific_heat_capacity"); // J kg-1 K-1
61 m_c_w = config.get_number("constants.fresh_water.specific_heat_capacity"); // J kg-1 K-1
62 m_L = config.get_number("constants.fresh_water.latent_heat_of_fusion"); // J kg-1
63 m_T_melting = config.get_number("constants.fresh_water.melting_point_temperature"); // K
64 m_T_tolerance = config.get_number("enthalpy_converter.relaxed_is_temperate_tolerance"); // K
65 m_T_0 = config.get_number("enthalpy_converter.T_reference"); // K
66
67 m_cold_mode = member(config.get_string("energy.model"), {"cold", "none"});
68}
69
70//! Return `true` if ice at `(E, P)` is temperate.
71//! Determines if E >= E_s(p), that is, if the ice is at the pressure-melting point.
72bool EnthalpyConverter::is_temperate(double E, double P) const {
73 if (m_cold_mode) {
74 return is_temperate_relaxed(E, P);
75 }
76
77 return (E >= enthalpy_cts(P));
78}
79
80//! A relaxed version of `is_temperate()`.
81/*! Returns `true` if the pressure melting temperature corresponding to `(E, P)` is within
82 `enthalpy_converter.relaxed_is_temperate_tolerance` from `fresh_water.melting_point_temperature`.
83 */
84bool EnthalpyConverter::is_temperate_relaxed(double E, double P) const {
86}
87
88void EnthalpyConverter::validate_T_omega_P(double T, double omega, double P) {
89#if (Pism_DEBUG==1)
90 const double T_melting = melting_temperature(P);
91 if (T <= 0.0) {
92 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T = %f <= 0 is not a valid absolute temperature",T);
93 }
94 const double eps = 1.0e-6;
95 if ((omega < 0.0 - eps) || (1.0 + eps < omega)) {
96 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "water fraction omega=%f not in range [0,1]",omega);
97 }
98 if (T > T_melting + eps) {
99 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T=%f exceeds T_melting=%f; not allowed",T,T_melting);
100 }
101 if ((T < T_melting - eps) && (omega > 0.0 + eps)) {
102 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T < T_melting AND omega > 0 is contradictory;"
103 " got T=%f, T_melting=%f, omega=%f",
104 T, T_melting, omega);
105 }
106#else
107 (void) T;
108 (void) omega;
109 (void) P;
110#endif
111}
112
113void EnthalpyConverter::validate_E_P(double E, double P) {
114#if (Pism_DEBUG==1)
115 if (E >= enthalpy_liquid(P)) {
116 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "E=%f J/kg at P=%f Pa equals or exceeds that of liquid water (%f J/kg)",
117 E, P, enthalpy_liquid(P));
118 }
119#else
120 (void) E;
121 (void) P;
122#endif
123}
124
125
126//! Get pressure in ice from depth below surface using the hydrostatic assumption.
127/*! If \f$d\f$ is the depth then
128 \f[ p = p_{\text{air}} + \rho_i g d. \f]
129Frequently \f$d\f$ is computed from the thickess minus a level in the ice,
130something like `ice_thickness(i, j) - z[k]`. The input depth to this routine is allowed to
131be negative, representing a position above the surface of the ice.
132 */
133double EnthalpyConverter::pressure(double depth) const {
134 if (depth >= 0.0) {
135 return m_p_air + m_rho_i * m_g * depth;
136 }
137
138 return m_p_air; // at or above surface of ice
139}
140
141//! Compute pressure in a column of ice. Does not check validity of `depth`.
142void EnthalpyConverter::pressure(const std::vector<double> &depth,
143 unsigned int ks,
144 std::vector<double> &result) const {
145 for (unsigned int k = 0; k <= ks; ++k) {
146 result[k] = m_p_air + m_rho_i * m_g * depth[k];
147 }
148}
149
150//! Get melting temperature from pressure p.
151/*!
152 \f[ T_m(p) = T_{melting} - \beta p. \f]
153 */
155 return m_T_melting - m_beta * P;
156}
157
158
159//! @brief Compute the maximum allowed value of ice enthalpy
160//! (corresponds to @f$ \omega = 1 @f$).
162 return enthalpy_cts(P) + L(melting_temperature(P));
163}
164
165
166//! Get absolute (not pressure-adjusted) ice temperature (K) from enthalpy and pressure.
167/*! From \ref AschwandenBuelerKhroulevBlatter,
168 \f[ T= T(E,p) = \begin{cases}
169 c_i^{-1} E + T_0, & E < E_s(p), \\
170 T_m(p), & E_s(p) \le E < E_l(p).
171 \end{cases} \f]
172
173We do not allow liquid water (%i.e. water fraction \f$\omega=1.0\f$) so we
174throw an exception if \f$E \ge E_l(p)\f$.
175 */
176double EnthalpyConverter::temperature(double E, double P) const {
177#if (Pism_DEBUG==1)
178 validate_E_P(E, P);
179#endif
180
181 if (E < enthalpy_cts(P)) {
182 return temperature_cold(E);
183 }
184
185 return melting_temperature(P);
186}
187
188
189//! Get pressure-adjusted ice temperature, in kelvin, from enthalpy and pressure.
190/*!
191The pressure-adjusted temperature is:
192 \f[ T_{pa}(E,p) = T(E,p) - T_m(p) + T_{melting}. \f]
193 */
194double EnthalpyConverter::pressure_adjusted_temperature(double E, double P) const {
195 return temperature(E, P) - melting_temperature(P) + m_T_melting;
196}
197
198
199//! Get liquid water fraction from enthalpy and pressure.
200/*!
201 From [@ref AschwandenBuelerKhroulevBlatter],
202 @f[
203 \omega(E,p) =
204 \begin{cases}
205 0.0, & E \le E_s(p), \\
206 (E-E_s(p)) / L, & E_s(p) < E < E_l(p).
207 \end{cases}
208 @f]
209
210 We do not allow liquid water (i.e. water fraction @f$ \omega=1.0 @f$).
211 */
212double EnthalpyConverter::water_fraction(double E, double P) const {
213#if (Pism_DEBUG==1)
214 validate_E_P(E, P);
215#endif
216
217 double E_s = enthalpy_cts(P);
218 if (E <= E_s) {
219 return 0.0;
220 }
221
222 return (E - E_s) / L(melting_temperature(P));
223}
224
225
226//! Compute enthalpy from absolute temperature, liquid water fraction, and pressure.
227/*! This is an inverse function to the functions \f$T(E,p)\f$ and
228\f$\omega(E,p)\f$ [\ref AschwandenBuelerKhroulevBlatter]. It returns:
229 \f[E(T,\omega,p) =
230 \begin{cases}
231 c_i (T - T_0), & T < T_m(p) \quad\text{and}\quad \omega = 0, \\
232 E_s(p) + \omega L, & T = T_m(p) \quad\text{and}\quad \omega \ge 0.
233 \end{cases} \f]
234Certain cases are not allowed and throw exceptions:
235- \f$T<=0\f$ (error code 1)
236- \f$\omega < 0\f$ or \f$\omega > 1\f$ (error code 2)
237- \f$T>T_m(p)\f$ (error code 3)
238- \f$T<T_m(p)\f$ and \f$\omega > 0\f$ (error code 4)
239These inequalities may be violated in the sixth digit or so, however.
240 */
241double EnthalpyConverter::enthalpy(double T, double omega, double P) const {
242#if (Pism_DEBUG==1)
243 validate_T_omega_P(T, omega, P);
244#endif
245
246 const double T_melting = melting_temperature(P);
247
248 if (T < T_melting) {
249 return enthalpy_cold(T);
250 }
251
252 return enthalpy_cts(P) + omega * L(T_melting);
253}
254
255
256//! Compute enthalpy more permissively than enthalpy().
257/*! Computes enthalpy from absolute temperature, liquid water fraction, and
258pressure. Use this form of enthalpy() when outside sources (e.g. information
259from a coupler) might generate a temperature above the pressure melting point or
260generate cold ice with a positive water fraction.
261
262Treats temperatures above pressure-melting point as \e at the pressure-melting
263point. Interprets contradictory case of \f$T < T_m(p)\f$ and \f$\omega > 0\f$
264as cold ice, ignoring the water fraction (\f$\omega\f$) value.
265
266Calls enthalpy(), which validates its inputs.
267
268Computes:
269 @f[
270 E =
271 \begin{cases}
272 E(T,0.0,p), & T < T_m(p) \quad \text{and} \quad \omega \ge 0, \\
273 E(T_m(p),\omega,p), & T \ge T_m(p) \quad \text{and} \quad \omega \ge 0,
274 \end{cases}
275 @f]
276 but ensures @f$ 0 \le \omega \le 1 @f$ in second case. Calls enthalpy() for
277 @f$ E(T,\omega,p) @f$.
278 */
279double EnthalpyConverter::enthalpy_permissive(double T, double omega, double P) const {
280 const double T_m = melting_temperature(P);
281
282 if (T < T_m) {
283 return enthalpy(T, 0.0, P);
284 }
285
286 // T >= T_m(P) replaced with T = T_m(P)
287 return enthalpy(T_m, std::max(0.0, std::min(omega, 1.0)), P);
288}
289
291 : EnthalpyConverter(config) {
292 // turn on the "cold" enthalpy converter mode
293 m_cold_mode = true;
294 // set melting temperature to one million kelvin so that all ice is cold
295 m_T_melting = 1e6;
296 // disable pressure-dependence of the melting temperature by setting Clausius-Clapeyron beta to
297 // zero
298 m_beta = 0.0;
299}
300
301//! Latent heat of fusion of water as a function of pressure melting
302//! temperature.
303/*!
304
305 Following a re-interpretation of [@ref
306 AschwandenBuelerKhroulevBlatter], we require that @f$ \Diff{H}{p} =
307 0 @f$:
308
309 @f[
310 \Diff{H}{p} = \diff{H_w}{p} + \diff{H_w}{p}\Diff{T}{p}
311 @f]
312
313 We assume that water is incompressible, so @f$ \Diff{T}{p} = 0 @f$
314 and the second term vanishes.
315
316 As for the first term, equation (5) of [@ref
317 AschwandenBuelerKhroulevBlatter] defines @f$ H_w @f$ as follows:
318
319 @f[
320 H_w = \int_{T_0}^{T_m(p)} C_i(t) dt + L + \int_{T_m(p)}^T C_w(t)dt
321 @f]
322
323 Using the fundamental theorem of Calculus, we get
324 @f[
325 \diff{H_w}{p} = (C_i(T_m(p)) - C_w(T_m(p))) \diff{T_m(p)}{p} + \diff{L}{p}
326 @f]
327
328 Assuming that @f$ C_i(T) = c_i @f$ and @f$ C_w(T) = c_w @f$ (i.e. specific heat
329 capacities of ice and water do not depend on temperature) and using
330 the Clausius-Clapeyron relation
331 @f[
332 T_m(p) = T_m(p_{\text{air}}) - \beta p,
333 @f]
334
335 we get
336 @f{align}{
337 \Diff{H}{p} &= (c_i - c_w)\diff{T_m(p)}{p} + \diff{L}{p}\\
338 &= \beta(c_w - c_i) + \diff{L}{p}\\
339 @f}
340 Requiring @f$ \Diff{H}{p} = 0 @f$ implies
341 @f[
342 \diff{L}{p} = -\beta(c_w - c_i),
343 @f]
344 and so
345 @f{align}{
346 L(p) &= -\beta p (c_w - c_i) + C\\
347 &= (T_m(p) - T_m(p_{\text{air}})) (c_w - c_i) + C.
348 @f}
349
350 Letting @f$ p = p_{\text{air}} @f$ we find @f$ C = L(p_\text{air}) = L_0 @f$, so
351 @f[
352 L(p) = (T_m(p) - T_m(p_{\text{air}})) (c_w - c_i) + L_0,
353 @f]
354 where @f$ L_0 @f$ is the latent heat of fusion of water at atmospheric
355 pressure.
356
357 Therefore a consistent interpretation of [@ref
358 AschwandenBuelerKhroulevBlatter] requires the temperature-dependent
359 approximation of the latent heat of fusion of water given above.
360
361 Note that this form of @f$ L(p) @f$ also follows from Kirchhoff's
362 law of thermochemistry.
363*/
364double EnthalpyConverter::L(double T_pm) const {
365 return m_L + (m_c_w - m_c_i) * (T_pm - 273.15);
366}
367
368//! Specific heat capacity of ice.
369double EnthalpyConverter::c() const {
370 return m_c_i;
371}
372
373//! Get enthalpy E_s(p) at cold-temperate transition point from pressure p.
374/*! Returns
375 \f[ E_s(p) = c_i (T_m(p) - T_0), \f]
376 */
377double EnthalpyConverter::enthalpy_cts(double P) const {
378 return m_c_i * (melting_temperature(P) - m_T_0);
379}
380
381//! Convert temperature into enthalpy (cold case).
382double EnthalpyConverter::enthalpy_cold(double T) const {
383 return m_c_i * (T - m_T_0);
384}
385
386//! Convert enthalpy into temperature (cold case).
388 return (E / m_c_i) + m_T_0;
389}
390
391} // end of namespace pism
ColdEnthalpyConverter(const Config &config)
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
std::string get_string(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
double enthalpy(double T, double omega, double P) const
Compute enthalpy from absolute temperature, liquid water fraction, and pressure.
EnthalpyConverter(const Config &config)
double m_rho_i
density of ice
bool is_temperate(double E, double P) const
static void validate_T_omega_P(double T, double omega, double P)
double melting_temperature(double P) const
Get melting temperature from pressure p.
double m_T_tolerance
temperature tolerance used in is_temperate() in cold ice mode
double temperature(double E, double P) const
Get absolute (not pressure-adjusted) ice temperature (K) from enthalpy and pressure.
double m_beta
beta in the Clausius-Clapeyron relation ( ).
bool m_cold_mode
if cold ice methods are selected, use is_temperate() check based on temperature, not enthalpy
double enthalpy_cold(double T) const
Convert temperature into enthalpy (cold case).
double enthalpy_liquid(double P) const
Compute the maximum allowed value of ice enthalpy (corresponds to ).
double m_c_i
specific heat capacity of ice
double c() const
Specific heat capacity of ice.
double m_g
acceleration due to gravity
double enthalpy_permissive(double T, double omega, double P) const
Compute enthalpy more permissively than enthalpy().
double m_L
latent heat of fusion of water at atmospheric pressure
double pressure_adjusted_temperature(double E, double P) const
Get pressure-adjusted ice temperature, in kelvin, from enthalpy and pressure.
double enthalpy_cts(double P) const
Get enthalpy E_s(p) at cold-temperate transition point from pressure p.
bool is_temperate_relaxed(double E, double P) const
A relaxed version of is_temperate().
double temperature_cold(double E) const
Convert enthalpy into temperature (cold case).
double pressure(double depth) const
Get pressure in ice from depth below surface using the hydrostatic assumption.
double m_T_melting
melting temperature of pure water at atmospheric pressure
double water_fraction(double E, double P) const
Get liquid water fraction from enthalpy and pressure.
double m_c_w
specific heat capacity of pure water
double m_p_air
atmospheric pressure
static void validate_E_P(double E, double P)
double L(double T_pm) const
double m_T_0
reference temperature in the definition of ice enthalpy
Converts between specific enthalpy and temperature or liquid content.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
#define PISM_ERROR_LOCATION
static const double L
Definition exactTestL.cc:40
static const double k
Definition exactTestP.cc:42
bool member(const std::string &string, const std::set< std::string > &set)