Ice rheology

The “rheology” of a viscous fluid refers to the relation between the applied stress and the resulting deformation, the strain rate. The models of ice rheology available in PISM are all isotropic [43]. A rheology in this class is described by a “flow law”, which is, in the most general case in PISM, a function \(F(\sigma,T,\omega,P,d)\) in the “constitutive relation” form

(8)\[D_{ij} = F(\sigma,T,\omega,P,d)\, \sigma_{ij}'.\]

Here \(D_{ij}\) is the strain rate tensor, \(\sigma_{ij}'\) is the stress deviator tensor, \(T\) is the ice temperature, \(\omega\) is the liquid water fraction, \(P\) is the pressure, \(d\) is the grain size, and \(\sigma^2 = \frac{1}{2} \|\sigma_{ij}'\|_F = \frac{1}{2} \sigma_{ij}' \sigma_{ij}'\) defines the second invariant \(\sigma\) of the stress deviator tensor.

Form (8) of the flow law is used in the SIA, but the “viscosity form” of a flow law, found by inverting the constitutive relation (8), is needed for ice shelf and ice stream (SSA) flow and the first-order stress balance approximation [10]:

(9)\[\sigma_{ij}' = 2 \nu(D,T,\omega,P,d)\,D_{ij}\]

Here \(\nu(D,T,\omega,P,d)\) is the “effective viscosity” and \(D^2 = \frac{1}{2} D_{ij} D_{ij}\).

Most of the flow laws in PISM are of Glen-Nye single-power type. For example,

(10)\[F(\sigma,T) = A(T) \sigma^{n-1}\]

is the common temperature-dependent Glen law [36], [16] (which has no dependence on liquid water fraction, pressure, or grain size). If the ice softness \(A(T)=A_0\) is constant then the law is isothermal, whereas if there is dependence on temperature then \(A(T)\) is usually a generalization of “Arrhenius” form

\[A(T) = A \exp(-Q/(R T)).\]

The more elaborate Goldsby-Kohlstedt law [70] is a function \(F(\sigma,T,P,d)\), but in this case the function \(F\) cannot be factored into a product of a function of \(T,P,d\) and a single power of \(\sigma\), as in form (10).

There is only one choice for the flow law which takes full advantage of the enthalpy mode of PISM, which is the thermodynamical modeling (i.e. conservation of energy) default. Namely the Glen-Paterson-Budd-Lliboutry-Duval flow law [22], [35], [36], which is a function \(F(\sigma,T,\omega,P)\). This law is the only one in the literature where the ice softness depends on both the temperature and the liquid water fraction, so it parameterizes the (observed) softening of pressure-melting-temperature ice as its liquid fraction increases. One can use this default polythermal law or one may choose among a number of “cold ice” laws listed below which do not use the liquid water fraction.

Flow law choices

Configuration parameters

choose which flow law is used by the SIA, SSA, and the Blatter stress balances models, respectively. Allowed arguments are listed below.

  1. gpbld: Glen-Paterson-Budd-Lliboutry-Duval law [35], the enthalpy-based default in PISM [22]. Extends the Paterson-Budd law (below) to positive liquid water fraction. If \(A_{c}(T)\) is from Paterson-Budd then this law returns

    \[A(T,\omega) = A_{c}(T) (1 + C \omega),\]

    where \(\omega\) is the liquid water fraction, \(C\) is a configuration parameter flow_law­.gpbld­.water_frac_coeff, and \(\omega\) is capped at level flow_law­.gpbld­.water_frac_observed_limit.


    This flow law uses all the parameters controlling the Paterson-Budd law, plus the ones listed below.

    Prefix: flow_law.gpbld.

    1. water_frac_coeff (181.25) coefficient in Glen-Paterson-Budd flow law for extra dependence of softness on liquid water fraction (omega) [55], [35]

    2. water_frac_observed_limit (0.01) maximum value of liquid water fraction omega for which softness values are parameterized by [35]; used in Glen-Paterson-Budd-Lliboutry-Duval flow law; compare [22]

  2. pb: Paterson-Budd law, the cold-mode default. Fixed Glen exponent \(n=3\). Has a split “Arrhenius” term \(A(T) = A \exp(-Q/RT^*)\) where

    \[ \begin{align}\begin{aligned}A &= 3.615 \times 10^{-13}\, \text{s}^{-1}\, \text{Pa}^{-3},\\Q &= 6.0 \times 10^4\, \text{J}\, \text{mol}^{-1}\end{aligned}\end{align} \]

    if \(T^* < T_{\text{critical}}\) and

    \[ \begin{align}\begin{aligned}A &= 1.733 \times 10^{3}\, \text{s}^{-1}\, \text{Pa}^{-3},\\Q &= 13.9 \times 10^4\, \text{J}\, \text{mol}^{-1}\end{aligned}\end{align} \]

    if \(T^* > T_{\text{critical}}\). Here \(T^*\) is pressure-adjusted temperature [36].


    Prefix: flow_law.Paterson_Budd.

    1. A_cold (3.61e-13 Pascal-3 / second) Paterson-Budd \(A_\text{cold}\), see [36]

    2. A_warm (1730 Pascal-3 / second) Paterson-Budd \(A_\text{warm}\), see [36]

    3. Q_cold (60000 Joule / mol) Paterson-Budd \(Q_\text{cold}\), see [36]

    4. Q_warm (139000 Joule / mol) Paterson-Budd \(Q_\text{warm}\), see [36]

    5. T_critical (263.15 Kelvin) Paterson-Budd critical temperature, see [36]

  3. arr: Cold part of Paterson-Budd. Regardless of temperature, the \(A\) and \(Q\) values for \(T^* < T_{\text{critical}}\) in the Paterson-Budd law apply. This is the flow law used in the thermomechanically-coupled exact solutions run by pismv -test F and pismv -test G [16], [78].

  4. arrwarm: Warm part of Paterson-Budd. Regardless of temperature, the \(A\) and \(Q\) values for \(T^* > T_{\text{critical}}\) in Paterson-Budd apply.

  5. hooke: Hooke law with

    \[A(T) = A \exp\left(-\frac{Q}{RT^*} + 3C (T_r - T^*)^\kappa\right).\]

    Fixed Glen exponent \(n=3\) and constants as in [79], [46].


    Prefix: flow_law.Hooke.

    1. A (4.42165e-09 Pascal-3 second-1) \(A_{\text{Hooke}} = (1/B_0)^n\) where n=3 and \(B_0\) = 1.928 \(a^{1/3}\) Pa. See [79]

    2. C (0.16612 Kelvin^k) See [79]

    3. Q (78800 Joule / mol) Activation energy, see [79]

    4. Tr (273.39 Kelvin) See [79]

    5. k (1.17) See [79]

  6. isothermal_glen: The isothermal Glen flow law.


    (11)\[ \begin{align}\begin{aligned}F(\sigma) &= A_0 \sigma^{n-1},\\\nu(D) &= \frac{1}{2} B_0 D^{(1-n)/(2n)},\end{aligned}\end{align} \]

    where \(A_0\) is the ice softness and \(B_0=A_0^{-1/n}\) is the ice hardness.


    Prefix: flow_law.isothermal_Glen.

    1. ice_softness (3.1689e-24 Pascal-3 second-1) ice softness used by the isothermal Glen flow law [59]

  7. gk: The Goldsby-Kohlstedt flow law. This law has a combination of exponents from \(n=1.8\) to \(n=4\) [70].


    The viscosity form (9) of this flow law is not known, so it can only be used by the SIA stress balance.

    Because it has more than one power, stress_balance­.sia­.Glen_exponent has no effect, though stress_balance­.sia­.enhancement_factor works as expected. This law does not use the liquid water fraction, but only the temperature.

    Constants defining this flow law are hard-wired in the implementation. Please see the source code for details.

Enhancement factor and exponent

An enhancement factor can be added to any flow law. Single-power laws also permit control of the flow law exponent.

The parameter stress_balance­.sia­.enhancement_factor sets \(e\) in

(12)\[D_{ij} = e\, F(\sigma,T,\omega,P,d)\, \sigma_{ij}',\]

see (8).

Parameters stress_balance­.ssa­.enhancement_factor and stress_balance­.blatter­.enhancement_factor set \(e\) in

(13)\[\sigma_{ij}' = e^{-1/n}\, 2\, \nu(D,T,\omega,P,d)\, D_{ij},\]

see (9).

Parameters stress_balance­.sia­.Glen_exponent, stress_balance­.ssa­.Glen_exponent, stress_balance­.blatter­.Glen_exponent set the exponent when a single-power flow law is used.

Simply changing to a different value from the default \(n=3\) is not recommended without a corresponding change to the enhancement factor, however. This is because the coefficient and the power are non-trivially linked when a power law is fit to experimental data [34], [36].

Here is a possible approach to adjusting both the enhancement factor and the exponent. Suppose \(\sigma_0\) is preferred as a scale (reference) for the driving stress that appears in both SIA and SSA models. Typically this is on the order of one bar or \(10^5\) Pa. Suppose one wants the same amount of deformation \(D_0\) at this reference driving stress as one changes from the old exponent \(n_{old}\) to the new exponent \(n_{new}\). That is, suppose one wants

\[ \begin{align}\begin{aligned}D_0 &= E_{old}\, A\, \sigma_0^{n_{old}},\\D_0 &= E_{new}\, A\, \sigma_0^{n_{new}}\end{aligned}\end{align} \]

to be true with a new enhancement factor \(E_{new}\). Eliminating \(D_0\) and solving for the new enhancement factor gives

(14)\[E_{new} = E_{old}\, \sigma_0^{n_{old} - n_{new}}.\]

It follows, for example, that if one has a run with values

-sia_e 3.0 -sia_n 3.0

then a new run with exponent \(n=6.0\) and the same deformation at the reference driving stress of \(10^5\) Pa will use

-sia_e 3.0e-15 -sia_n 6.0

because \(E_{new} = 3.0 \sigma_0^{3-6} = 3.0 \times (10^5)^{-3}\) from equation (14).

A corresponding formula applies to changing the enhancement factor for the SSA and Blatter stress balance models.


  1. [2] used \(e_{\text{SIA}} = 3.0\) for Greenland ice sheet simulations (see the supplement) while [52] used \(e_{\text{SIA}} = 4.5\) for simulations of the Antarctic ice sheet with PISM-PIK.

  2. [52] used \(e_{\text{SSA}} =0.512\) for simulations of the Antarctic ice sheet with PISM-PIK.

Previous Up Next