EnthalpyConverter Class Reference

Converts ice specific enthalpy to-and-from temperature and water content. More...

#include <enthalpyConverter.hh>

Inheritance diagram for EnthalpyConverter:

Public Member Functions
	EnthalpyConverter (const NCConfigVariable &config)
virtual	~EnthalpyConverter ()
virtual PetscErrorCode	viewConstants (PetscViewer viewer) const
	Simple view of state of EnthalpyConverter. viewer==NULL sends to stdout.
virtual double	getPressureFromDepth (double depth) const
	Get pressure in ice from depth below surface using the hydrostatic assumption.
virtual double	getMeltingTemp (double p) const
	Get melting temperature from pressure p.
virtual double	getEnthalpyCTS (double p) const
	Get enthalpy E_s(p) at cold-temperate transition point, which normalizes enthalpy, from pressure p.
virtual PetscErrorCode	getEnthalpyInterval (double p, double &E_s, double &E_l) const
	Get enthalpies E_s(p) and E_l(p) (endpoints of temperate ice enthalpy range) from pressure p.
virtual double	getCTS (double E, double p) const
	Computes the ratio CTS = E / E_s(p). The cold-temperate transition surface (CTS) is the level surface CTS = 1.
virtual bool	isTemperate (double E, double p) const
	Determines if E >= E_s(p), that is, if the ice is at the pressure-melting point.
virtual bool	isLiquified (double E, double p) const
	Is ice at given enthalpy and pressure actually liquid water?
virtual PetscErrorCode	getAbsTemp (double E, double p, double &T) const
	Get absolute ice temperature (K) from enthalpy and pressure.
virtual PetscErrorCode	getPATemp (double E, double p, double &T_pa) const
	Get pressure-adjusted ice temperature (K) from enthalpy and pressure.
virtual PetscErrorCode	getWaterFraction (double E, double p, double &omega) const
	Get liquid water fraction from enthalpy and pressure.
virtual double	getWaterFractionLimited (double E, double p) const
	Get liquid water fraction from enthalpy and pressure, but return omega=1 if high enthalpy.
virtual PetscErrorCode	getEnth (double T, double omega, double p, double &E) const
	Compute enthalpy from absolute temperature, liquid water fraction, and pressure.
virtual PetscErrorCode	getEnthPermissive (double T, double omega, double p, double &E) const
	Compute enthalpy more permissively from temperature and water fraction. Compare getEnth().
virtual PetscErrorCode	getEnthAtWaterFraction (double omega, double p, double &E) const
	Returns enthalpy for temperate ice with a given liquid fraction.
Protected Attributes
double	T_0
double	L
double	c_i
double	rho_i
double	g
double	p_air
double	beta

Detailed Description

Converts ice specific enthalpy to-and-from temperature and water content.

ONCE THIS CLASS IS CONFIRMED TO BE CORRECT, IT *CAN* GO IN src/base/pism_const.{hh|cc}

Use this way, for example within IceModel with NCConfigVariable config member:

  #include "enthalpyConverter.hh"
  EnthalpyConverter EC(&config);  // this runs constructor, so make sure this
                                  //   is done after the creation/initialization of
                                  //   NCConfigVariable config
  ...
  for (...) {
    ...
    E_s = EC.getEnthalpyCTS(p);
    ... etc ...
  }

Some methods are functions which return the computed values or a boolean. These do no error checking. Others return PetscErrorCode (and set arguments). These check either that the enthalpy is below that of liquid water, or that the temperature (in K) is positive.

Specifically, getAbsTemp() gives return value 1 if the input enthalpy exceeded that of liquid water, but puts the temperature of maximum-liquid-content temperate ice into its computed value for T. And getWaterFraction() gives return value of 1 under the same condition, but puts the maximum-liquid-content into its computed value for omega.

The three methods that get the enthalpy from temperatures and liquid fractions, namely getEnth, getEnthPermissive, getEnthAtWaterFraction, are more strict about error checking. They call SETERRQ() if there arguments are invalid.

Definition at line 57 of file enthalpyConverter.hh.

Constructor & Destructor Documentation

EnthalpyConverter ( const NCConfigVariable & config )

Definition at line 24 of file enthalpyConverter.cc.

References beta, c_i, g, NCConfigVariable.get(), L, p_air, rho_i, and T_0.

virtual ~EnthalpyConverter ( ) [virtual]

Definition at line 60 of file enthalpyConverter.hh.

Member Function Documentation

PetscErrorCode getAbsTemp	(	double	E,
		double	p,
		double &	T
	)			const `[virtual]`

Get absolute ice temperature (K) from enthalpy and pressure.

From AschwandenBlatter, equation (12)

$T= T(E,p) = \begin{cases} c_i^{-1} (E-E_s(p)) + T_m(p), & E < E_s(p), \\ T_m(p), & E_s(p) \le E < E_l(p). \end{cases}$

But the first case simplifies if we expand $E_s$ :

$c_i^{-1} (E-E_s(p)) + T_m(p) = c_i^{-1} (E-c_i T_m(p)) + T_m(p) = c_i^{-1} E.$

We do not allow liquid water (i.e. water fraction $\omega=1.0$ ) so we return a 1 as output PetscErrorCode if $E \ge E_l(p)$ .

Reimplemented in ICMEnthalpyConverter.

Definition at line 164 of file enthalpyConverter.cc.

References c_i, getEnthalpyInterval(), and getMeltingTemp().

Referenced by IceModel.basalSlidingHeatingSIA(), IceModel.compute_temp(), IceModel.computeEffectiveViscosity(), IceModel.correctSigma(), IceModel.energyStats(), IceModel.enthalpyAndDrainageStep(), getPATemp(), and IceModel.velocitySIAStaggered().

double getCTS	(	double	E,
		double	p
	)			const `[virtual]`

Computes the ratio CTS = E / E_s(p). The cold-temperate transition surface (CTS) is the level surface CTS = 1.

If $E$ and $E_s$ are the enthalpy and the CTS enthalpy, respectively, then

$CTS = \frac{E}{E_s}.$

where CTS = 1 is the CTS contour line. Thus CTS is greater than one for temperate ice and less than one for cold ice. Presumably only used for postprocessing.

Definition at line 129 of file enthalpyConverter.cc.

References getEnthalpyCTS().

Referenced by IceModel.setCTSFromEnthalpy().

PetscErrorCode getEnth	(	double	T,
		double	omega,
		double	p,
		double &	E
	)			const `[virtual]`

Compute enthalpy from absolute temperature, liquid water fraction, and pressure.

This is an inverse function to the functions $T(E,p)$ and $\omega(E,p)$ . It returns this enthalpy value:

$E(T,\omega,p) = \begin{cases} E_s(p) + c_i (T-T_m(p)), & T \le 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}$

Certain cases are not allowed and return errors:

is not allowed,
and $\omega > 0$ is not allowed,
$\omega < 0$ or $\omega > 1$ is not allowed.

Because of these not-allowed cases, the following expression is also valid:

$E(T,\omega,p) = E_s(p) + c_i (T-T_m(p)) + \omega L.$

Reimplemented in ICMEnthalpyConverter.

Definition at line 292 of file enthalpyConverter.cc.

References c_i, e, getEnthalpyCTS(), getMeltingTemp(), and L.

Referenced by getEnthPermissive().

double getEnthalpyCTS ( double p ) const [virtual]

Get enthalpy E_s(p) at cold-temperate transition point, which normalizes enthalpy, from pressure p.

Returns

$E_s(p) = c_i T_m(p),$

In particular,

$E_s( p_{\text{air}}) = c_i (T_0 - \beta p_{\text{air}}) = 548743.22\, \frac{\text{J}}{\text{kg}}$

is the enthalpy for surface ice.

Note that $p_{\text{air}} = 10^5$ Pa at temperature $T_0=273.15$ K, with standard choices $c_i=2009$ J kg-1 K-1 and $\beta=7.53\times 10^{-8}$ K Pa-1.

Definition at line 104 of file enthalpyConverter.cc.

References c_i, and getMeltingTemp().

Referenced by getCTS(), getEnth(), IceModel.getEnthalpyCTSColumn(), getEnthalpyInterval(), getEnthAtWaterFraction(), and isTemperate().

PetscErrorCode getEnthalpyInterval	(	double	p,
		double &	E_s,
		double &	E_l
	)			const `[virtual]`

Get enthalpies E_s(p) and E_l(p) (endpoints of temperate ice enthalpy range) from pressure p.

Ice at enthalpy $E$ is temperate if $E_s(p) < E < E_l(p)$ :

$E_s(p) = c_i T_m(p),$

$E_l(p) = H_s(p) + L.$

Definition at line 114 of file enthalpyConverter.cc.

References getEnthalpyCTS(), and L.

Referenced by getAbsTemp(), getWaterFraction(), getWaterFractionLimited(), isLiquified(), and PolyThermalGPBLDIce.softnessParameterFromEnth().

PetscErrorCode getEnthAtWaterFraction	(	double	omega,
		double	p,
		double &	E
	)			const `[virtual]`

Returns enthalpy for temperate ice with a given liquid fraction.

Computes

$E = E_s(p) + \omega L.$

Definition at line 356 of file enthalpyConverter.cc.

References e, getEnthalpyCTS(), and L.

Referenced by IceModel.enthalpyAndDrainageStep().

PetscErrorCode getEnthPermissive	(	double	T,
		double	omega,
		double	p,
		double &	E
	)			const `[virtual]`

Compute enthalpy more permissively from temperature and water fraction. Compare getEnth().

Computes enthalpy from absolute temperature, liquid water fraction, and pressure as before. Use this form of getEnth() when outside sources (e.g. information from a coupler) might generate a temperature above the pressure melting point or cold ice with a positive water fraction.

Treats temperatures above pressure-melting point as at the pressure-melting point. Interprets contradictory case of $T < T_m(p)$ and $\omega > 0$ as cold ice, ignoring water fraction $\omega > 0$ .

Computes:

$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}$

but ensures $0\le \omega \le 1$ in second case. Calls getEnth() for $E(T,\omega,p)$ .

Checks that $T > 0$ and returns error if not.

Reimplemented in ICMEnthalpyConverter.

Definition at line 334 of file enthalpyConverter.cc.

References getEnth(), and getMeltingTemp().

Referenced by IceModel.check_maximum_thickness(), IceModel.enthalpyAndDrainageStep(), IceModel.putTempAtDepth(), IceModel.setEnth3FromT3_ColdIce(), and IceModel.setEnth3FromT3AndLiqfrac3().

double getMeltingTemp ( double p ) const [virtual]

Get melting temperature from pressure p.

$T_m(p) = T_0 - \beta p.$

Definition at line 88 of file enthalpyConverter.cc.

References beta, and T_0.

Referenced by IceModel.enthalpyAndDrainageStep(), getAbsTemp(), getEnth(), getEnthalpyCTS(), getEnthPermissive(), and getPATemp().

PetscErrorCode getPATemp	(	double	E,
		double	p,
		double &	T_pa
	)			const `[virtual]`

Get pressure-adjusted ice temperature (K) from enthalpy and pressure.

Calls getAbsTemp(), which computes $T(E,p)$ , and getMeltingTemp(), which computes $T_m(p)$ .

We define the pressure-adjusted temperature to be:

$T_{pa} = T_{pa}(E,p) = T(E,p) - T_m(p) + T_0.$

Definition at line 193 of file enthalpyConverter.cc.

References getAbsTemp(), getMeltingTemp(), and T_0.

Referenced by IceModel.compute_temp_pa(), and PolyThermalGPBLDIce.softnessParameterFromEnth().

double getPressureFromDepth ( double depth ) const [virtual]

Get pressure in ice from depth below surface using the hydrostatic assumption.

If $d$ is the depth then

$p = p_{\text{air}} + \rho_i g d.$

Frequently $d$ is computed from the thickess minus a level in the ice, something like "H[i][j] - z[k]". The input depth to this routine is allowed to be negative, representing a position above the surface of the ice.

Definition at line 75 of file enthalpyConverter.cc.

References g, p_air, and rho_i.

Referenced by IceModel.basalSlidingHeatingSIA(), IceModel.check_maximum_thickness(), IceModel.compute_temp(), IceModel.compute_temp_pa(), IceModel.computeEffectiveViscosity(), IceModel.correctSigma(), PolyThermalGPBLDIce.effectiveViscosityColumnFromEnth(), IceModel.energyStats(), IceModel.enthalpyAndDrainageStep(), IceModel.getEnthalpyCTSColumn(), IceModel.putTempAtDepth(), IceModel.setCTSFromEnthalpy(), IceModel.setEnth3FromT3_ColdIce(), IceModel.setEnth3FromT3AndLiqfrac3(), and IceModel.setLiquidFracFromEnthalpy().