EnthalpyConverter.cc

Go to the documentation of this file.

// Copyright (C) 2009-2017, 2021, 2023, 2024, 2025 Andreas Aschwanden, Ed Bueler and Constantine Khroulev
//
// This file is part of PISM.
//
// PISM is free software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the Free Software
// Foundation; either version 3 of the License, or (at your option) any later
// version.
//
// PISM is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License
// along with PISM; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
#include "pism/util/EnthalpyConverter.hh"
#include "pism/util/ConfigInterface.hh"
 
#include "pism/util/error_handling.hh"
#include "pism/util/pism_utilities.hh"
 
namespace pism {
 
/*! @class EnthalpyConverter
 
  Maps from @f$ (H,p) @f$ to @f$ (T,\omega,p) @f$ and back.
 
  Requirements:
 
  1. A converter has to implement an invertible map @f$ (H,p) \to (T,
     \omega, p) @f$ *and* its inverse. Both the forward map and the
     inverse need to be defined for all permissible values of @f$
     (H,p) @f$ and @f$ (T, \omega, p) @f$, respectively.
 
  2. A converter has to be consistent with laws and parameterizations
     used elsewhere in the model. This includes models coupled to
     PISM.
 
  3. For a fixed volume of liquid (or solid) water and given two
     energy states @f$ (H_1, p_1) @f$ and @f$ (H_2, p_2) @f$ , let @f$
     \Delta U_H @f$ be the difference in internal energy of this
     volume between the two states *computed using enthalpy*. We require
     that @f$ \Delta U_T = \Delta U_H @f$ , where @f$ \Delta U_T @f$
     is the difference in internal energy *computed using corresponding*
     @f$ (T_1, \omega_1, p_1) @f$ *and* @f$ (T_2, \omega_2, p_2) @f$.
 
  4. We assume that ice and water are incompressible, so a change in
     pressure does no work, and @f$ \Diff{H}{p} = 0 @f$. In addition
     to this, for cold ice and liquid water @f$ \Diff{T}{p} = 0 @f$.
*/
 
EnthalpyConverter::EnthalpyConverter(const Config &config) {
  m_p_air       = config.get_number("surface.pressure"); // Pa
  m_g           = config.get_number("constants.standard_gravity"); // m s-2
  m_beta        = config.get_number("constants.ice.beta_Clausius_Clapeyron"); // K Pa-1
  m_rho_i       = config.get_number("constants.ice.density"); // kg m-3
  m_c_i         = config.get_number("constants.ice.specific_heat_capacity"); // J kg-1 K-1
  m_c_w         = config.get_number("constants.fresh_water.specific_heat_capacity"); // J kg-1 K-1
  m_L           = config.get_number("constants.fresh_water.latent_heat_of_fusion"); // J kg-1
  m_T_melting   = config.get_number("constants.fresh_water.melting_point_temperature"); // K
  m_T_tolerance = config.get_number("enthalpy_converter.relaxed_is_temperate_tolerance"); // K
  m_T_0         = config.get_number("enthalpy_converter.T_reference"); // K
 
  m_cold_mode = member(config.get_string("energy.model"), {"cold", "none"});
}
EnthalpyConverter::EnthalpyConverter(const Config &config) {…}
 
//! Return `true` if ice at `(E, P)` is temperate.
//! Determines if E >= E_s(p), that is, if the ice is at the pressure-melting point.
bool EnthalpyConverter::is_temperate(double E, double P) const {
  if (m_cold_mode) {
    return is_temperate_relaxed(E, P);
  }
 
  return (E >= enthalpy_cts(P));
}
bool EnthalpyConverter::is_temperate(double E, double P) const  {…}
 
//! A relaxed version of `is_temperate()`.
/*! Returns `true` if the pressure melting temperature corresponding to `(E, P)` is within
    `enthalpy_converter.relaxed_is_temperate_tolerance` from `fresh_water.melting_point_temperature`.
 */
bool EnthalpyConverter::is_temperate_relaxed(double E, double P) const {
  return (pressure_adjusted_temperature(E, P) >= m_T_melting - m_T_tolerance);
}
bool EnthalpyConverter::is_temperate_relaxed(double E, double P) const  {…}
 
void EnthalpyConverter::validate_T_omega_P(double T, double omega, double P) {
#if (Pism_DEBUG==1)
  const double T_melting = melting_temperature(P);
  if (T <= 0.0) {
    throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T = %f <= 0 is not a valid absolute temperature",T);
  }
  const double eps = 1.0e-6;
  if ((omega < 0.0 - eps) || (1.0 + eps < omega)) {
    throw RuntimeError::formatted(PISM_ERROR_LOCATION, "water fraction omega=%f not in range [0,1]",omega);
  }
  if (T > T_melting + eps) {
    throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T=%f exceeds T_melting=%f; not allowed",T,T_melting);
  }
  if ((T < T_melting - eps) && (omega > 0.0 + eps)) {
    throw RuntimeError::formatted(PISM_ERROR_LOCATION, "T < T_melting AND omega > 0 is contradictory;"
                                  " got T=%f, T_melting=%f, omega=%f",
                                  T, T_melting, omega);
  }
#else
  (void) T;
  (void) omega;
  (void) P;
#endif
}
void EnthalpyConverter::validate_T_omega_P(double T, double omega, double P) {…}
 
void EnthalpyConverter::validate_E_P(double E, double P) {
#if (Pism_DEBUG==1)
  if (E >= enthalpy_liquid(P)) {
    throw RuntimeError::formatted(PISM_ERROR_LOCATION, "E=%f J/kg at P=%f Pa equals or exceeds that of liquid water (%f J/kg)",
                                  E, P, enthalpy_liquid(P));
  }
#else
  (void) E;
  (void) P;
#endif
}
void EnthalpyConverter::validate_E_P(double E, double P) {…}
 
 
//! Get pressure in ice from depth below surface using the hydrostatic assumption.
/*! If \f$d\f$ is the depth then
     \f[ p = p_{\text{air}}  + \rho_i g d. \f]
Frequently \f$d\f$ is computed from the thickess minus a level in the ice,
something like `ice_thickness(i, j) - z[k]`.  The input depth to this routine is allowed to
be negative, representing a position above the surface of the ice.
 */
double EnthalpyConverter::pressure(double depth) const {
  if (depth >= 0.0) {
    return m_p_air + m_rho_i * m_g * depth;
  }
 
  return m_p_air; // at or above surface of ice
}
double EnthalpyConverter::pressure(double depth) const  {…}
 
//! Compute pressure in a column of ice. Does not check validity of `depth`.
void EnthalpyConverter::pressure(const std::vector<double> &depth,
                                 unsigned int ks,
                                 std::vector<double> &result) const {
  for (unsigned int k = 0; k <= ks; ++k) {
    result[k] = m_p_air + m_rho_i * m_g * depth[k];
  }
}
void EnthalpyConverter::pressure(const std::vector<double> &depth, {…}
 
//! Get melting temperature from pressure p.
/*!
     \f[ T_m(p) = T_{melting} - \beta p. \f]
 */
double EnthalpyConverter::melting_temperature(double P) const {
  return m_T_melting - m_beta * P;
}
double EnthalpyConverter::melting_temperature(double P) const  {…}
 
 
//! @brief Compute the maximum allowed value of ice enthalpy
//! (corresponds to @f$ \omega = 1 @f$).
double EnthalpyConverter::enthalpy_liquid(double P) const {
  return enthalpy_cts(P) + L(melting_temperature(P));
}
double EnthalpyConverter::enthalpy_liquid(double P) const  {…}
 
 
//! Get absolute (not pressure-adjusted) ice temperature (K) from enthalpy and pressure.
/*! From \ref AschwandenBuelerKhroulevBlatter,
     \f[ T= T(E,p) = \begin{cases}
                       c_i^{-1} E + T_0,  &  E < E_s(p), \\
                       T_m(p),            &  E_s(p) \le E < E_l(p).
                     \end{cases} \f]
 
We do not allow liquid water (%i.e. water fraction \f$\omega=1.0\f$) so we
throw an exception if \f$E \ge E_l(p)\f$.
 */
double EnthalpyConverter::temperature(double E, double P) const {
#if (Pism_DEBUG==1)
  validate_E_P(E, P);
#endif
 
  if (E < enthalpy_cts(P)) {
    return temperature_cold(E);
  }
 
  return melting_temperature(P);
}
double EnthalpyConverter::temperature(double E, double P) const  {…}
 
 
//! Get pressure-adjusted ice temperature, in kelvin, from enthalpy and pressure.
/*!
The pressure-adjusted temperature is:
     \f[ T_{pa}(E,p) = T(E,p) - T_m(p) + T_{melting}. \f]
 */
double EnthalpyConverter::pressure_adjusted_temperature(double E, double P) const {
  return temperature(E, P) - melting_temperature(P) + m_T_melting;
}
double EnthalpyConverter::pressure_adjusted_temperature(double E, double P) const  {…}
 
 
//! Get liquid water fraction from enthalpy and pressure.
/*!
  From [@ref AschwandenBuelerKhroulevBlatter],
   @f[
   \omega(E,p) =
   \begin{cases}
     0.0, & E \le E_s(p), \\
     (E-E_s(p)) / L, & E_s(p) < E < E_l(p).
   \end{cases}
   @f]
 
   We do not allow liquid water (i.e. water fraction @f$ \omega=1.0 @f$).
 */
double EnthalpyConverter::water_fraction(double E, double P) const {
#if (Pism_DEBUG==1)
  validate_E_P(E, P);
#endif
 
  double E_s = enthalpy_cts(P);
  if (E <= E_s) {
    return 0.0;
  }
 
  return (E - E_s) / L(melting_temperature(P));
}
double EnthalpyConverter::water_fraction(double E, double P) const  {…}
 
 
//! Compute enthalpy from absolute temperature, liquid water fraction, and pressure.
/*! This is an inverse function to the functions \f$T(E,p)\f$ and
\f$\omega(E,p)\f$ [\ref AschwandenBuelerKhroulevBlatter].  It returns:
  \f[E(T,\omega,p) =
       \begin{cases}
         c_i (T - T_0),     & T < T_m(p) \quad\text{and}\quad \omega = 0, \\
         E_s(p) + \omega L, & T = T_m(p) \quad\text{and}\quad \omega \ge 0.
       \end{cases} \f]
Certain cases are not allowed and throw exceptions:
- \f$T<=0\f$ (error code 1)
- \f$\omega < 0\f$ or \f$\omega > 1\f$ (error code 2)
- \f$T>T_m(p)\f$ (error code 3)
- \f$T<T_m(p)\f$ and \f$\omega > 0\f$ (error code 4)
These inequalities may be violated in the sixth digit or so, however.
 */
double EnthalpyConverter::enthalpy(double T, double omega, double P) const {
#if (Pism_DEBUG==1)
  validate_T_omega_P(T, omega, P);
#endif
 
  const double T_melting = melting_temperature(P);
 
  if (T < T_melting) {
    return enthalpy_cold(T);
  }
 
  return enthalpy_cts(P) + omega * L(T_melting);
}
double EnthalpyConverter::enthalpy(double T, double omega, double P) const  {…}
 
 
//! Compute enthalpy more permissively than enthalpy().
/*! Computes enthalpy from absolute temperature, liquid water fraction, and
pressure.  Use this form of enthalpy() when outside sources (e.g. information
from a coupler) might generate a temperature above the pressure melting point or
generate cold ice with a positive water fraction.
 
Treats temperatures above pressure-melting point as \e at the pressure-melting
point.  Interprets contradictory case of \f$T < T_m(p)\f$ and \f$\omega > 0\f$
as cold ice, ignoring the water fraction (\f$\omega\f$) value.
 
Calls enthalpy(), which validates its inputs.
 
Computes:
  @f[
  E =
  \begin{cases}
    E(T,0.0,p),         & T < T_m(p) \quad \text{and} \quad \omega \ge 0, \\
    E(T_m(p),\omega,p), & T \ge T_m(p) \quad \text{and} \quad \omega \ge 0,
  \end{cases}
  @f]
  but ensures @f$ 0 \le \omega \le 1 @f$ in second case.  Calls enthalpy() for
  @f$ E(T,\omega,p) @f$.
 */
double EnthalpyConverter::enthalpy_permissive(double T, double omega, double P) const {
  const double T_m = melting_temperature(P);
 
  if (T < T_m) {
    return enthalpy(T, 0.0, P);
  }
 
  // T >= T_m(P) replaced with T = T_m(P)
  return enthalpy(T_m, std::max(0.0, std::min(omega, 1.0)), P);
}
double EnthalpyConverter::enthalpy_permissive(double T, double omega, double P) const  {…}
 
ColdEnthalpyConverter::ColdEnthalpyConverter(const Config &config)
  : EnthalpyConverter(config) {
  // turn on the "cold" enthalpy converter mode
  m_cold_mode = true;
  // set melting temperature to one million kelvin so that all ice is cold
  m_T_melting = 1e6;
  // disable pressure-dependence of the melting temperature by setting Clausius-Clapeyron beta to
  // zero
  m_beta = 0.0;
}
ColdEnthalpyConverter::ColdEnthalpyConverter(const Config &config) {…}
 
//! Latent heat of fusion of water as a function of pressure melting
//! temperature.
/*!
 
  Following a re-interpretation of [@ref
  AschwandenBuelerKhroulevBlatter], we require that @f$ \Diff{H}{p} =
  0 @f$:
 
  @f[
  \Diff{H}{p} = \diff{H_w}{p} + \diff{H_w}{p}\Diff{T}{p}
  @f]
 
  We assume that water is incompressible, so @f$ \Diff{T}{p} = 0 @f$
  and the second term vanishes.
 
  As for the first term, equation (5) of [@ref
  AschwandenBuelerKhroulevBlatter] defines @f$ H_w @f$ as follows:
 
  @f[
  H_w = \int_{T_0}^{T_m(p)} C_i(t) dt + L + \int_{T_m(p)}^T C_w(t)dt
  @f]
 
  Using the fundamental theorem of Calculus, we get
  @f[
  \diff{H_w}{p} = (C_i(T_m(p)) - C_w(T_m(p))) \diff{T_m(p)}{p} + \diff{L}{p}
  @f]
 
  Assuming that @f$ C_i(T) = c_i @f$ and @f$ C_w(T) = c_w @f$ (i.e. specific heat
  capacities of ice and water do not depend on temperature) and using
  the Clausius-Clapeyron relation
  @f[
  T_m(p) = T_m(p_{\text{air}}) - \beta p,
  @f]
 
  we get
  @f{align}{
  \Diff{H}{p} &= (c_i - c_w)\diff{T_m(p)}{p} + \diff{L}{p}\\
  &= \beta(c_w - c_i) + \diff{L}{p}\\
  @f}
  Requiring @f$ \Diff{H}{p} = 0 @f$ implies
  @f[
  \diff{L}{p} = -\beta(c_w - c_i),
  @f]
  and so
  @f{align}{
  L(p) &= -\beta p (c_w - c_i) + C\\
  &= (T_m(p) - T_m(p_{\text{air}})) (c_w - c_i) + C.
  @f}
 
  Letting @f$ p = p_{\text{air}} @f$ we find @f$ C = L(p_\text{air}) = L_0 @f$, so
  @f[
  L(p) = (T_m(p) - T_m(p_{\text{air}})) (c_w - c_i) + L_0,
  @f]
  where @f$ L_0 @f$ is the latent heat of fusion of water at atmospheric
  pressure.
 
  Therefore a consistent interpretation of [@ref
  AschwandenBuelerKhroulevBlatter] requires the temperature-dependent
  approximation of the latent heat of fusion of water given above.
 
  Note that this form of @f$ L(p) @f$ also follows from Kirchhoff's
  law of thermochemistry.
*/
double EnthalpyConverter::L(double T_pm) const {
  return m_L + (m_c_w - m_c_i) * (T_pm - 273.15);
}
double EnthalpyConverter::L(double T_pm) const  {…}
 
//! Specific heat capacity of ice.
double EnthalpyConverter::c() const {
  return m_c_i;
}
double EnthalpyConverter::c() const  {…}
 
//! Get enthalpy E_s(p) at cold-temperate transition point from pressure p.
/*! Returns
     \f[ E_s(p) = c_i (T_m(p) - T_0), \f]
 */
double EnthalpyConverter::enthalpy_cts(double P) const {
  return m_c_i * (melting_temperature(P) - m_T_0);
}
double EnthalpyConverter::enthalpy_cts(double P) const  {…}
 
//! Convert temperature into enthalpy (cold case).
double EnthalpyConverter::enthalpy_cold(double T) const {
  return m_c_i * (T - m_T_0);
}
double EnthalpyConverter::enthalpy_cold(double T) const  {…}
 
//! Convert enthalpy into temperature (cold case).
double EnthalpyConverter::temperature_cold(double E) const {
  return (E / m_c_i) + m_T_0;
}
double EnthalpyConverter::temperature_cold(double E) const  {…}
 
} // end of namespace pism