tEnthalpyConverter.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tEnthalpyConverter.cc (13494B)
---
1 // Copyright (C) 2009-2017 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/pism_utilities.hh"
20 #include "EnthalpyConverter.hh"
21 #include "pism/util/ConfigInterface.hh"
22
23 #include "pism/util/error_handling.hh"
24
25 namespace 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
55 EnthalpyConverter::EnthalpyConverter(const Config &config) {
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_do_cold_ice_methods = config.get_flag("energy.temperature_based");
68 }
69
70 EnthalpyConverter::~EnthalpyConverter() {
71 // empty
72 }
73
74 //! Return `true` if ice at `(E, P)` is temperate.
75 //! Determines if E >= E_s(p), that is, if the ice is at the pressure-melting point.
76 bool EnthalpyConverter::is_temperate(double E, double P) const {
77 if (m_do_cold_ice_methods) {
78 return is_temperate_relaxed(E, P);
79 } else {
80 return (E >= enthalpy_cts(P));
81 }
82 }
83
84 //! A relaxed version of `is_temperate()`.
85 /*! Returns `true` if the pressure melting temperature corresponding to `(E, P)` is within
86 `enthalpy_converter.relaxed_is_temperate_tolerance` from `fresh_water.melting_point_temperature`.
87 */
88 bool EnthalpyConverter::is_temperate_relaxed(double E, double P) const {
89 return (pressure_adjusted_temperature(E, P) >= m_T_melting - m_T_tolerance);
90 }
91
92
93 void EnthalpyConverter::validate_T_omega_P(double T, double omega, double P) const {
94 #if (Pism_DEBUG==1)
95 const double T_melting = melting_temperature(P);
96 if (T <= 0.0) {
97 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T = %f <= 0 is not a valid absolute temperature",T);
98 }
99 if ((omega < 0.0 - 1.0e-6) || (1.0 + 1.0e-6 < omega)) {
100 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "water fraction omega=%f not in range [0,1]",omega);
101 }
102 if (T > T_melting + 1.0e-6) {
103 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T=%f exceeds T_melting=%f; not allowed",T,T_melting);
104 }
105 if ((T < T_melting - 1.0e-6) && (omega > 0.0 + 1.0e-6)) {
106 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T < T_melting AND omega > 0 is contradictory;"
107 " got T=%f, T_melting=%f, omega=%f",
108 T, T_melting, omega);
109 }
110 #else
111 (void) T;
112 (void) omega;
113 (void) P;
114 #endif
115 }
116
117 void EnthalpyConverter::validate_E_P(double E, double P) const {
118 #if (Pism_DEBUG==1)
119 if (E >= enthalpy_liquid(P)) {
120 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "E=%f J/kg at P=%f Pa equals or exceeds that of liquid water (%f J/kg)",
121 E, P, enthalpy_liquid(P));
122 }
123 #else
124 (void) E;
125 (void) P;
126 #endif
127 }
128
129
130 //! Get pressure in ice from depth below surface using the hydrostatic assumption.
131 /*! If \f$d\f$ is the depth then
132 \f[ p = p_{\text{air}} + \rho_i g d. \f]
133 Frequently \f$d\f$ is computed from the thickess minus a level in the ice,
134 something like `ice_thickness(i, j) - z[k]`. The input depth to this routine is allowed to
135 be negative, representing a position above the surface of the ice.
136 */
137 double EnthalpyConverter::pressure(double depth) const {
138 if (depth >= 0.0) {
139 return m_p_air + m_rho_i * m_g * depth;
140 } else {
141 return m_p_air; // at or above surface of ice
142 }
143 }
144
145 //! Compute pressure in a column of ice. Does not check validity of `depth`.
146 void EnthalpyConverter::pressure(const std::vector<double> &depth,
147 unsigned int ks,
148 std::vector<double> &result) const {
149 for (unsigned int k = 0; k <= ks; ++k) {
150 result[k] = m_p_air + m_rho_i * m_g * depth[k];
151 }
152 }
153
154 //! Get melting temperature from pressure p.
155 /*!
156 \f[ T_m(p) = T_{melting} - \beta p. \f]
157 */
158 double EnthalpyConverter::melting_temperature(double P) const {
159 return m_T_melting - m_beta * P;
160 }
161
162
163 //! @brief Compute the maximum allowed value of ice enthalpy
164 //! (corresponds to @f$ \omega = 1 @f$).
165 double EnthalpyConverter::enthalpy_liquid(double P) const {
166 return enthalpy_cts(P) + L(melting_temperature(P));
167 }
168
169
170 //! Get absolute (not pressure-adjusted) ice temperature (K) from enthalpy and pressure.
171 /*! From \ref AschwandenBuelerKhroulevBlatter,
172 \f[ T= T(E,p) = \begin{cases}
173 c_i^{-1} E + T_0, & E < E_s(p), \\
174 T_m(p), & E_s(p) \le E < E_l(p).
175 \end{cases} \f]
176
177 We do not allow liquid water (%i.e. water fraction \f$\omega=1.0\f$) so we
178 throw an exception if \f$E \ge E_l(p)\f$.
179 */
180 double EnthalpyConverter::temperature(double E, double P) const {
181 validate_E_P(E, P);
182
183 if (E < enthalpy_cts(P)) {
184 return temperature_cold(E);
185 } else {
186 return melting_temperature(P);
187 }
188 }
189
190
191 //! Get pressure-adjusted ice temperature, in Kelvin, from enthalpy and pressure.
192 /*!
193 The pressure-adjusted temperature is:
194 \f[ T_{pa}(E,p) = T(E,p) - T_m(p) + T_{melting}. \f]
195 */
196 double EnthalpyConverter::pressure_adjusted_temperature(double E, double P) const {
197 return temperature(E, P) - melting_temperature(P) + m_T_melting;
198 }
199
200
201 //! Get liquid water fraction from enthalpy and pressure.
202 /*!
203 From [@ref AschwandenBuelerKhroulevBlatter],
204 @f[
205 \omega(E,p) =
206 \begin{cases}
207 0.0, & E \le E_s(p), \\
208 (E-E_s(p)) / L, & E_s(p) < E < E_l(p).
209 \end{cases}
210 @f]
211
212 We do not allow liquid water (i.e. water fraction @f$ \omega=1.0 @f$).
213 */
214 double EnthalpyConverter::water_fraction(double E, double P) const {
215 validate_E_P(E, P);
216
217 double E_s = enthalpy_cts(P);
218 if (E <= E_s) {
219 return 0.0;
220 } else {
221 return (E - E_s) / L(melting_temperature(P));
222 }
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]
234 Certain 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)
239 These inequalities may be violated in the sixth digit or so, however.
240 */
241 double EnthalpyConverter::enthalpy(double T, double omega, double P) const {
242 validate_T_omega_P(T, omega, P);
243
244 const double T_melting = melting_temperature(P);
245
246 if (T < T_melting) {
247 return enthalpy_cold(T);
248 } else {
249 return enthalpy_cts(P) + omega * L(T_melting);
250 }
251 }
252
253
254 //! Compute enthalpy more permissively than enthalpy().
255 /*! Computes enthalpy from absolute temperature, liquid water fraction, and
256 pressure. Use this form of enthalpy() when outside sources (e.g. information
257 from a coupler) might generate a temperature above the pressure melting point or
258 generate cold ice with a positive water fraction.
259
260 Treats temperatures above pressure-melting point as \e at the pressure-melting
261 point. Interprets contradictory case of \f$T < T_m(p)\f$ and \f$\omega > 0\f$
262 as cold ice, ignoring the water fraction (\f$\omega\f$) value.
263
264 Calls enthalpy(), which validates its inputs.
265
266 Computes:
267 @f[
268 E =
269 \begin{cases}
270 E(T,0.0,p), & T < T_m(p) \quad \text{and} \quad \omega \ge 0, \\
271 E(T_m(p),\omega,p), & T \ge T_m(p) \quad \text{and} \quad \omega \ge 0,
272 \end{cases}
273 @f]
274 but ensures @f$ 0 \le \omega \le 1 @f$ in second case. Calls enthalpy() for
275 @f$ E(T,\omega,p) @f$.
276 */
277 double EnthalpyConverter::enthalpy_permissive(double T, double omega, double P) const {
278 const double T_m = melting_temperature(P);
279
280 if (T < T_m) {
281 return enthalpy(T, 0.0, P);
282 } else { // T >= T_m(P) replaced with T = T_m(P)
283 return enthalpy(T_m, std::max(0.0, std::min(omega, 1.0)), P);
284 }
285 }
286
287 ColdEnthalpyConverter::ColdEnthalpyConverter(const Config &config)
288 : EnthalpyConverter(config) {
289 // turn on the "cold" enthalpy converter mode
290 m_do_cold_ice_methods = true;
291 // set melting temperature to one million Kelvin so that all ice is cold
292 m_T_melting = 1e6;
293 // disable pressure-dependence of the melting temperature by setting Clausius-Clapeyron beta to
294 // zero
295 m_beta = 0.0;
296 }
297
298 ColdEnthalpyConverter::~ColdEnthalpyConverter() {
299 // empty
300 }
301
302 //! Latent heat of fusion of water as a function of pressure melting
303 //! temperature.
304 /*!
305
306 Following a re-interpretation of [@ref
307 AschwandenBuelerKhroulevBlatter], we require that @f$ \Diff{H}{p} =
308 0 @f$:
309
310 @f[
311 \Diff{H}{p} = \diff{H_w}{p} + \diff{H_w}{p}\Diff{T}{p}
312 @f]
313
314 We assume that water is incompressible, so @f$ \Diff{T}{p} = 0 @f$
315 and the second term vanishes.
316
317 As for the first term, equation (5) of [@ref
318 AschwandenBuelerKhroulevBlatter] defines @f$ H_w @f$ as follows:
319
320 @f[
321 H_w = \int_{T_0}^{T_m(p)} C_i(t) dt + L + \int_{T_m(p)}^T C_w(t)dt
322 @f]
323
324 Using the fundamental theorem of Calculus, we get
325 @f[
326 \diff{H_w}{p} = (C_i(T_m(p)) - C_w(T_m(p))) \diff{T_m(p)}{p} + \diff{L}{p}
327 @f]
328
329 Assuming that @f$ C_i(T) = c_i @f$ and @f$ C_w(T) = c_w @f$ (i.e. specific heat
330 capacities of ice and water do not depend on temperature) and using
331 the Clausius-Clapeyron relation
332 @f[
333 T_m(p) = T_m(p_{\text{air}}) - \beta p,
334 @f]
335
336 we get
337 @f{align}{
338 \Diff{H}{p} &= (c_i - c_w)\diff{T_m(p)}{p} + \diff{L}{p}\\
339 &= \beta(c_w - c_i) + \diff{L}{p}\\
340 @f}
341 Requiring @f$ \Diff{H}{p} = 0 @f$ implies
342 @f[
343 \diff{L}{p} = -\beta(c_w - c_i),
344 @f]
345 and so
346 @f{align}{
347 L(p) &= -\beta p (c_w - c_i) + C\\
348 &= (T_m(p) - T_m(p_{\text{air}})) (c_w - c_i) + C.
349 @f}
350
351 Letting @f$ p = p_{\text{air}} @f$ we find @f$ C = L(p_\text{air}) = L_0 @f$, so
352 @f[
353 L(p) = (T_m(p) - T_m(p_{\text{air}})) (c_w - c_i) + L_0,
354 @f]
355 where @f$ L_0 @f$ is the latent heat of fusion of water at atmospheric
356 pressure.
357
358 Therefore a consistent interpretation of [@ref
359 AschwandenBuelerKhroulevBlatter] requires the temperature-dependent
360 approximation of the latent heat of fusion of water given above.
361
362 Note that this form of @f$ L(p) @f$ also follows from Kirchhoff's
363 law of thermochemistry.
364 */
365 double EnthalpyConverter::L(double T_pm) const {
366 return m_L + (m_c_w - m_c_i) * (T_pm - 273.15);
367 }
368
369 //! Specific heat capacity of ice.
370 double EnthalpyConverter::c() const {
371 return m_c_i;
372 }
373
374 //! Get enthalpy E_s(p) at cold-temperate transition point from pressure p.
375 /*! Returns
376 \f[ E_s(p) = c_i (T_m(p) - T_0), \f]
377 */
378 double EnthalpyConverter::enthalpy_cts(double P) const {
379 return m_c_i * (melting_temperature(P) - m_T_0);
380 }
381
382 //! Convert temperature into enthalpy (cold case).
383 double EnthalpyConverter::enthalpy_cold(double T) const {
384 return m_c_i * (T - m_T_0);
385 }
386
387 //! Convert enthalpy into temperature (cold case).
388 double EnthalpyConverter::temperature_cold(double E) const {
389 return (E / m_c_i) + m_T_0;
390 }
391
392 } // end of namespace pism