trheology.rst - 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
---
trheology.rst (10225B)
---
1 .. include:: ../../../global.txt
2
3 .. _sec-rheology:
4
5 Ice rheology
6 ------------
7
8 .. contents::
9
10 The "rheology" of a viscous fluid refers to the relation between the applied stress and
11 the resulting deformation, the strain rate. The models of ice rheology available in PISM
12 are all isotropic :cite:`Paterson`. A rheology in this class is described by a "flow law",
13 which is, in the most general case in PISM, a function `F(\sigma,T,\omega,P,d)` in the
14 "constitutive relation" form
15
16 .. math::
17 :label: eq-constitutive
18
19 D_{ij} = F(\sigma,T,\omega,P,d)\, \sigma_{ij}'.
20
21 Here `D_{ij}` is the strain rate tensor, `\sigma_{ij}'` is the stress deviator tensor, `T`
22 is the ice temperature, `\omega` is the liquid water fraction, `P` is the pressure, `d` is
23 the grain size, and `\sigma^2 = \frac{1}{2} \|\sigma_{ij}'\|_F = \frac{1}{2} \sigma_{ij}'
24 \sigma_{ij}'` defines the second invariant `\sigma` of the stress deviator tensor.
25
26 Form :eq:`eq-constitutive` of the flow law is used in the SIA, but the "viscosity" form of
27 a flow law, found by inverting the constitutive relation :eq:`eq-constitutive`, is needed
28 for ice shelf and ice stream (SSA) flow :cite:`BBssasliding`:
29
30 .. math::
31 :label: eq-viscosityform
32
33 \sigma_{ij}' = 2 \nu(D,T,\omega,P,d)\,D_{ij}
34
35 Here `\nu(D,T,\omega,P,d)` is the "effective viscosity" and `D^2 = \frac{1}{2}
36 D_{ij} D_{ij}`.
37
38 Most of the flow laws in PISM are of Glen-Nye single-power type. For example,
39
40 .. math::
41 :label: eq-glen
42
43 F(\sigma,T) = A(T) \sigma^{n-1}
44
45 is the common temperature-dependent Glen law :cite:`PatersonBudd`, :cite:`BBL` (which has no
46 dependence on liquid water fraction, pressure, or grain size). If the ice softness
47 `A(T)=A_0` is constant then the law is isothermal, whereas if there is dependence on
48 temperature then `A(T)` is usually a generalization of "Arrhenius" form
49
50 .. math::
51
52 A(T) = A \exp(-Q/(R T)).
53
54 The more elaborate Goldsby-Kohlstedt law :cite:`GoldsbyKohlstedt` is a function
55 `F(\sigma,T,P,d)`, but in this case the function `F` cannot be factored into a
56 product of a function of `T,P,d` and a single power of `\sigma`, as in form
57 :eq:`eq-glen`.
58
59 There is only one choice for the flow law which takes full advantage of the enthalpy mode
60 of PISM, which is the thermodynamical modeling (i.e. conservation of energy) default.
61 Namely the Glen-Paterson-Budd-Lliboutry-Duval flow law
62 :cite:`AschwandenBuelerKhroulevBlatter`, :cite:`LliboutryDuval1985`, :cite:`PatersonBudd`,
63 which is a function `F(\sigma,T,\omega,P)`. This law is the only one in the literature
64 where the ice softness depends on both the temperature and the liquid water fraction, so
65 it parameterizes the (observed) softening of pressure-melting-temperature ice as its
66 liquid fraction increases. One can use this default polythermal law or one may choose
67 among a number of "cold ice" laws listed in :numref:`tab-flowlaw` which do not use the
68 liquid water fraction.
69
70 All flow law parameters can be changed using configuration parameters; see section
71 :ref:`sec-pism-defaults` and the implementation of flow laws in the `Source Code Browser
72 <pism-browser_>`_. Note that different flow laws have different numbers of parameters, but
73 all have at least two parameters (e.g. `A_0` and `n` in ``isothermal_glen``). One can
74 create a new, and reasonably arbitrarily, scalar function `F` by modifying source code;
75 see source files in ``src/base/rheology/``.
76
77 Choosing the flow laws for SIA and SSA stress balances
78 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
79
80 Command-line options :opt:`-sia_flow_law` and :opt:`-ssa_flow_law` choose which flow law
81 is used by the SIA and SSA stress balances, respectively. Allowed arguments are listed in
82 :numref:`tab-flowlaw` below. Viscosity form :eq:`eq-viscosityform` is not known for the
83 Goldsby-Kohlstedt law :cite:`GoldsbyKohlstedt`, so option "``-ssa_flow_law gk``" is an
84 error.
85
86 .. list-table:: Single-power flow laws. Choose the ice rheology using ``-sia_flow_law``
87 and ``-ssa_flow_law`` and one of the names in this table. Flow law choices
88 other than ``gpbld`` do not use the liquid water fraction `\omega`
89 but only the temperature `T`.
90 :name: tab-flowlaw
91 :header-rows: 1
92 :widths: 1,3
93
94 * - Name
95 - Comments and References
96
97 * - ``gpbld``
98 - Glen-Paterson-Budd-Lliboutry-Duval law :cite:`LliboutryDuval1985`, the enthalpy-based
99 default in PISM :cite:`AschwandenBuelerKhroulevBlatter`. Extends the Paterson-Budd law
100 (below) to positive liquid water fraction. If `A_{c}(T)` is from Paterson-Budd then
101 this law returns
102
103 `A(T,\omega) = A_{c}(T) (1 + C \omega),`
104
105 where `\omega` is the liquid water fraction, `C` is a configuration parameter
106 :config:`flow_law.gpbld.water_frac_coeff` [default `C=181.25`\], and `\omega` is
107 capped at level :config:`flow_law.gpbld.water_frac_observed_limit`.
108
109 * - ``pb``
110 - Paterson-Budd law, the cold-mode default. Fixed Glen exponent `n=3`. Has a split
111 "Arrhenius" term `A(T) = A \exp(-Q/RT^*)` where
112
113 `A = 3.615 \times 10^{-13}\, \text{s}^{-1}\, \text{Pa}^{-3},`
114 `Q = 6.0 \times 10^4\, \text{J}\, \text{mol}^{-1}`
115
116 if `T^* < 263` K and
117
118 `A = 1.733 \times 10^{3}\, \text{s}^{-1}\, \text{Pa}^{-3},`
119 `Q = 13.9 \times 10^4\, \text{J}\, \text{mol}^{-1}`
120
121 if `T^* > 263` K.
122
123 Here `T^*` is pressure-adjusted temperature :cite:`PatersonBudd`.
124
125 * - ``arr``
126 - *Cold* part of Paterson-Budd. Regardless of temperature, the `A` and `Q` values for
127 `T^*<263` K in the Paterson-Budd law apply. This is the flow law used in the
128 thermomechanically-coupled exact solutions run by ``pismv -test F`` and
129 ``pismv -test G`` :cite:`BBL`, :cite:`BB`.
130
131 * - ``arrwarm``
132 - *Warm* part of Paterson-Budd. Regardless of temperature, the `A` and `Q` values for
133 `T^*>263` K in Paterson-Budd apply.
134
135 * - ``hooke``
136 - Hooke law with
137
138 `A(T) = A \exp(-Q/(RT^*) + 3C (T_r - T^*)^\kappa).`
139
140 Fixed Glen exponent `n=3` and constants as in :cite:`Hooke`, :cite:`PayneBaldwin`.
141
142 * - ``isothermal_glen``
143 - The isothermal Glen flow law. Here `F(\sigma) = A_0 \sigma^{n-1}` with inverse
144 `\nu(D) = \frac{1}{2} B_0 D^{(1-n)/(2n)}` where `A_0` is the ice softness and
145 `B_0=A_0^{-1/n}` is the ice hardness.
146
147 * - ``gk``
148 - This law has a combination of exponents from `n=1.8` to `n=4`
149 :cite:`GoldsbyKohlstedt`. It can only be used by the SIA stress balance. Because it has
150 more than one power, option ``-sia_n`` has no effect, though ``-sia_e`` works as
151 expected. This law does not use the liquid water fraction, but only the
152 temperature.
153
154 Choose enhancement factor and exponent
155 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
156
157 An enhancement factor can be added to any flow law through a runtime option. Single-power
158 laws also permit control of the flow law exponent through a runtime option.
159
160 Options :opt:`-sia_e` and :opt:`-ssa_e` set flow enhancement factors for the SIA and SSA
161 respectively. Option ``-sia_e`` sets "`e`" in `D_{ij} = e\, F(\sigma,T,\omega,P,d)\,
162 \sigma_{ij}',` in equation :eq:`eq-constitutive`. Option ``-ssa_e`` sets "`e`" in the
163 viscosity form so that `\sigma_{ij}' = e^{-1/n}\, 2\, \nu(D,T,\omega,P,d)\, D_{ij}.`
164
165 Options :opt:`-sia_n` and :opt:`-ssa_n` set the exponent when a single-power flow law is
166 used (see :numref:`tab-flowlaw`). Simply changing to a different value from the default
167 `n=3` is not recommended without a corresponding change to the enhancement factor,
168 however. This is because the coefficient and the power are non-trivially linked when a
169 power law is fit to experimental data :cite:`CuffeyPaterson`, :cite:`PatersonBudd`.
170
171 Here is a possible approach to adjusting both the enhancement factor and the exponent.
172 Suppose `\sigma_0` is preferred as a scale (reference) for the driving stress that
173 appears in both SIA and SSA models. Typically this is on the order of one bar or
174 `10^5` Pa. Suppose one wants the same amount of deformation `D_0` at this
175 reference driving stress as one changes from the old exponent `n_{old}` to the new
176 exponent `n_{new}`. That is, suppose one wants both
177
178 .. math::
179
180 D_0 = E_{old}\, A\, \sigma_0^{n_{old}} \qquad \text{and} \qquad D_0
181 = E_{new}\, A\, \sigma_0^{n_{new}}
182
183 to be true with a new enhancement factor `E_{new}`. Eliminating `D_0` and
184 solving for the new enhancement factor gives
185
186 .. math::
187 :label: eq-renewexponent
188
189 E_{new} = E_{old}\, \sigma_0^{n_{old} - n_{new}}.
190
191 It follows, for example, that if one has a run with values
192
193 .. code-block:: none
194
195 -sia_e 3.0 -sia_n 3.0
196
197 then a new run with exponent `n=6.0` and the same deformation at the reference
198 driving stress of `10^5` Pa will use
199
200 .. code-block:: none
201
202 -sia_e 3.0e-15 -sia_n 6.0
203
204 because `E_{new} = 3.0 \sigma_0^{3-6} = 3.0 \times (10^5)^{-3}` from equation
205 :eq:`eq-renewexponent`.
206
207 A corresponding formula applies to ``-ssa_e`` if the ``-ssa_n`` value changes.
208
209 .. list-table:: For all flow laws, an enhancement factor can be added by a runtime option.
210 For the single-power flow laws in :numref:`tab-flowlaw`, the (Glen)
211 exponent can be controlled by a runtime option.
212 :name: tab-enhancementandexponent
213 :header-rows: 1
214 :widths: 1,2,2
215
216 * - Option
217 - Configuration parameter
218 - Comments
219
220 * - :opt:`-sia_e` (1.0)
221 - ``stress_balance.sia.enhancement_factor``
222 - Note (see the supplement of :cite:`AschwandenAdalgeirsdottirKhroulev`) used `3.0`
223 for Greenland ice sheet simulations while :cite:`Martinetal2011` used `4.5` for
224 simulations of the Antarctic ice sheet with PISM-PIK.
225
226 * - :opt:`-sia_n` (3.0)
227 - ``stress_balance.sia.Glen_exponent``
228 - See text and eqn :eq:`eq-renewexponent` to also set ``-sia_e`` if ``-sia_n`` changes.
229
230 * - :opt:`-ssa_e` (1.0)
231 - ``stress_balance.ssa.enhancement_factor``
232 - Note :cite:`Martinetal2011` used `0.512` for simulations of the Antarctic ice sheet with
233 PISM-PIK.
234
235 * - :opt:`-ssa_n` (3.0)
236 - ``stress_balance.ssa.Glen_exponent``
237 - See text and eqn :eq:`eq-renewexponent` to also set ``-ssa_e`` if ``-ssa_n``
238 changes.