tbasal-strength.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
---
tbasal-strength.rst (17222B)
---
1 .. include:: ../../../global.txt
2
3 .. include:: ../../../math-definitions.txt
4
5 .. _sec-basestrength:
6
7 Controlling basal strength
8 --------------------------
9
10 .. contents::
11
12 When using option :opt:`-stress_balance ssa+sia`, the SIA+SSA hybrid stress balance, a
13 model for basal resistance is required. This model for basal resistance is based, at least
14 conceptually, on the hypothesis that the ice sheet is underlain by a layer of till
15 :cite:`Clarke05`. The user can control the parts of this model:
16
17 - the so-called sliding law, typically a power law, which relates the ice base (sliding)
18 velocity to the basal shear stress, and which has a coefficient which is or has the
19 units of a yield stress,
20 - the model relating the effective pressure on the till layer to the yield stress of that
21 layer, and
22 - the model for relating the amount of water stored in the till to the effective pressure
23 on the till.
24
25 This subsection explains the relevant options.
26
27 The primary example of ``-stress_balance ssa+sia`` usage is in section :ref:`sec-start` of
28 this Manual, but the option is also used in sections :ref:`sec-MISMIP`,
29 :ref:`sec-MISMIP3d`, and :ref:`sec-jako`.
30
31 In PISM the key coefficient in the sliding is always denoted as yield stress `\tau_c`,
32 which is :var:`tauc` in PISM output files. This parameter represents the strength of the
33 aggregate material at the base of an ice sheet, a poorly-observed mixture of pressurized
34 liquid water, ice, granular till, and bedrock bumps. The yield stress concept also extends
35 to the power law form, and thus most standard sliding laws can be chosen by user options
36 (below). One reason that the yield stress is a useful parameter is that it can be
37 compared, when looking at PISM output files, to the driving stress (:var:`taud_mag` in
38 PISM output files). Specifically, where :var:`tauc` `<` :var:`taud_mag` you are likely to
39 see sliding if option ``-stress_balance ssa+sia`` is used.
40
41 A historical note on modeling basal sliding is in order. Sliding can be added directly to
42 a SIA stress balance model by making the sliding velocity a local function of the basal
43 value of the driving stress. Such an SIA sliding mechanism appears in ISMIP-HEINO
44 :cite:`Calovetal2009HEINOfinal` and in EISMINT II experiment H :cite:`EISMINT00`, among
45 other places. This kind of sliding is *not* recommended, as it does not make sense to
46 regard the driving stress as the local generator of flow if the bed is not holding all of
47 that stress :cite:`BBssasliding`, :cite:`Fowler01`. Within PISM, for historical reasons,
48 there is an implementation of SIA-based sliding only for verification test E; see section
49 :ref:`sec-verif`. PISM does *not* support this SIA-based sliding mode in other contexts.
50
51 Choosing the sliding law
52 ^^^^^^^^^^^^^^^^^^^^^^^^
53
54 In PISM the sliding law can be chosen to be a purely-plastic (Coulomb) model, namely,
55
56 .. math::
57 :label: eq-plastic
58
59 |\boldsymbol{\tau}_b| \le \tau_c \quad \text{and} \quad \boldsymbol{\tau}_b =
60 - \tau_c \frac{\mathbf{u}}{|\mathbf{u}|} \quad\text{if and only if}\quad |\mathbf{u}| > 0.
61
62 Equation :eq:`eq-plastic` says that the (vector) basal shear stress `\boldsymbol{\tau}_b`
63 is at most the yield stress `\tau_c`, and that only when the shear stress reaches the
64 yield value can there be sliding. The sliding law can, however, also be chosen to be the
65 power law
66
67 .. math::
68 :label: eq-pseudoplastic
69
70 \boldsymbol{\tau}_b = - \tau_c \frac{\mathbf{u}}{u_{\text{threshold}}^q |\mathbf{u}|^{1-q}},
71
72 where `u_{\text{threshold}}` is a parameter with units of velocity (see below). Condition
73 :eq:`eq-plastic` is studied in :cite:`SchoofStream` and :cite:`SchoofTill` in particular,
74 while power laws for sliding are common across the glaciological literature (e.g. see
75 :cite:`CuffeyPaterson`, :cite:`GreveBlatter2009`). Notice that the coefficient `\tau_c` in
76 :eq:`eq-pseudoplastic` has units of stress, regardless of the power `q`.
77
78 In both of the above equations :eq:`eq-plastic` and :eq:`eq-pseudoplastic` we call
79 `\tau_c` the *yield stress*. It corresponds to the variable :var:`tauc` in PISM output files.
80 We call the power law :eq:`eq-pseudoplastic` a "pseudo-plastic" law with power `q` and
81 threshold velocity `u_{\text{threshold}}`. At the threshold velocity the basal shear
82 stress `\boldsymbol{\tau}_b` has exact magnitude `\tau_c`. In equation
83 :eq:`eq-pseudoplastic`, `q` is the power controlled by ``-pseudo_plastic_q``, and the
84 threshold velocity `u_{\text{threshold}}` is controlled by ``-pseudo_plastic_uthreshold``.
85 The plastic model :eq:`eq-plastic` is the `q=0` case of :eq:`eq-pseudoplastic`.
86
87 See :numref:`tab-sliding-power-law` for options controlling the choice of sliding law. The
88 purely plastic case is the default; just use ``-stress_balance ssa+sia`` to turn it on.
89 (Or use ``-stress_balance ssa`` if a model with no vertical shear is desired.)
90
91 .. warning::
92
93 Options ``-pseudo_plastic_q`` and ``-pseudo_plastic_uthreshold`` have no effect if
94 ``-pseudo_plastic`` is not set.
95
96 .. list-table:: Sliding law command-line options
97 :name: tab-sliding-power-law
98 :header-rows: 1
99 :widths: 1,2
100
101 * - Option
102 - Description
103 * - :opt:`-pseudo_plastic`
104 - Enables the pseudo-plastic power law model. If this is not set the sliding law is
105 purely-plastic, so ``pseudo_plastic_q`` and ``pseudo_plastic_uthreshold`` are
106 inactive.
107 * - :opt:`-plastic_reg` (m/a)
108 - Set the value of `\epsilon` regularization of the plastic law, in the formula
109 `\boldsymbol{\tau}_b = - \tau_c \mathbf{u}/\sqrt{|\mathbf{u}|^2 + \epsilon^2}`. The
110 default is `0.01` m/a. This parameter is inactive if ``-pseudo_plastic`` is set.
111 * - :opt:`-pseudo_plastic_q`
112 - Set the exponent `q` in :eq:`eq-pseudoplastic`. The default is `0.25`.
113 * - :opt:`-pseudo_plastic_uthreshold` (m/a)
114 - Set `u_{\text{threshold}}` in :eq:`eq-pseudoplastic`. The default is `100` m/a.
115
116 Equation :eq:`eq-pseudoplastic` is a very flexible power law form. For example, the linear
117 case is `q=1`, in which case if `\beta=\tau_c/u_{\text{threshold}}` then the law is of the
118 form
119
120 .. math::
121
122 \boldsymbol{\tau}_b = - \beta \mathbf{u}
123
124 (The "`\beta`" coefficient is also called `\beta^2` in some sources (see :cite:`MacAyeal`,
125 for example).) If you want such a linear sliding law, and you have a value
126 `\beta=` ``beta`` in `\text{Pa}\,\text{s}\,\text{m}^{-1}`, then use this option
127 combination:
128
129 .. code-block:: none
130
131 -pseudo_plastic \
132 -pseudo_plastic_q 1.0 \
133 -pseudo_plastic_uthreshold 3.1556926e7 \
134 -yield_stress constant -tauc beta
135
136 This sets `u_{\text{threshold}}` to 1 `\text{m}\,\text{s}^{-1}` but using units
137 `\text{m}\,\text{a}^{-1}`.
138
139 More generally, it is common in the literature to see power-law sliding relations in the
140 form
141
142 .. math::
143
144 \boldsymbol{\tau}_b = - C |\mathbf{u}|^{m-1} \mathbf{u},
145
146 where `C` is a constant, as for example in sections :ref:`sec-MISMIP` and
147 :ref:`sec-MISMIP3d`. In that case, use this option combination:
148
149 .. code-block:: none
150
151 -pseudo_plastic \
152 -pseudo_plastic_q m \
153 -pseudo_plastic_uthreshold 3.1556926e7 \
154 -yield_stress constant \
155 -tauc C
156
157 .. _sec-lateral-drag:
158
159 Lateral drag
160 ~~~~~~~~~~~~
161
162 PISM prescribes lateral drag at ice margins next to ground with elevation above the ice.
163 (This is relevant in outlet glaciers flowing through fjords, valley glaciers, and next to
164 nunataks.) Set :config:`basal_resistance.beta_lateral_margin` to control the amount of
165 additional drag at these margins.
166
167 Determining the yield stress
168 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
169
170 Other than setting it to a constant, which only applies in some special cases, the above
171 discussion does not determine the yield stress `\tau_c`. As shown in
172 :numref:`tab-yieldstress`, there are two schemes for determining `\tau_c` in a
173 spatially-variable manner:
174
175 - ``-yield_stress mohr_coulomb`` (the default) determines the yields stress by models of
176 till material property (the till friction angle) and of the effective pressure on the
177 saturated till, or
178 - ``-yield_stress constant`` allows the yield stress to be supplied as time-independent
179 data, read from the input file.
180
181 In normal modelling cases, variations in yield stress are part of the explanation of the
182 locations of ice streams :cite:`SchoofStream`. The default model ``-yield_stress
183 mohr_coulomb`` determines these variations in time and space. The value of `\tau_c` is
184 determined in part by a subglacial hydrology model, including the modeled till-pore water
185 amount (`W_{till}`, NetCDF variable :var:`tillwat`; see :ref:`sec-subhydro`), which then
186 determines the effective pressure `N_{till}` (see below). The value of `\tau_c` is also
187 determined in part by a material property field `\phi`, the "till friction angle" (NetCDF
188 variable :var:`tillphi`). These quantities are related by the Mohr-Coulomb criterion
189 :cite:`CuffeyPaterson`:
190
191 .. math::
192 :label: eq-mohrcoulomb
193
194 \tau_c = c_{0} + (\tan\phi)\,N_{till}.
195
196 Here `c_0` is called the "till cohesion", whose default value in PISM is zero (see
197 :cite:`SchoofStream`, formula (2.4)) but which can be set by option :opt:`-till_cohesion`.
198
199 Option combination ``-yield_stress constant -tauc X`` can be used to fix the yield stress
200 to have value `\tau_c = X` at all grounded locations and all times if desired. This is
201 unlikely to be a good modelling choice for real ice sheets.
202
203 .. list-table:: Command-line options controlling how yield stress is determined
204 :name: tab-yieldstress
205 :header-rows: 1
206 :widths: 3,5
207
208 * - Option
209 - Description
210 * - :opt:`-yield_stress mohr_coulomb`
211 - The default. Use equation :eq:`eq-mohrcoulomb` to determine `\tau_c`. Only
212 effective if ``-stress_balance ssa`` or ``-stress_balance ssa+sia`` is also set.
213 * - :opt:`-till_cohesion`
214 - Set the value of the till cohesion (`c_{0}`) in the plastic till model. The value
215 is a pressure, given in Pa.
216 * - :opt:`-tauc_slippery_grounding_lines`
217 - If set, reduces the basal yield stress at grounded-below-sea-level grid points one
218 cell away from floating ice or ocean. Specifically, it replaces the
219 normally-computed `\tau_c` from the Mohr-Coulomb relation, which uses the effective
220 pressure from the modeled amount of water in the till, by the minimum value of
221 `\tau_c` from Mohr-Coulomb, i.e. using the effective pressure corresponding to the
222 maximum amount of till-stored water. Does not alter the reported amount of till
223 water, nor does this mechanism affect water conservation.
224 * - :opt:`-plastic_phi` (degrees)
225 - Use a constant till friction angle. The default is `30^{\circ}`.
226 * - :opt:`-topg_to_phi` (*list of 4 numbers*)
227 - Compute `\phi` using equation :eq:`eq-phipiecewise`.
228 * - :opt:`-yield_stress constant`
229 - Keep the current values of the till yield stress `\tau_c`. That is, do not update
230 them by the default model using the stored basal melt water. Only effective if
231 ``-stress_balance ssa`` or ``-stress_balance ssa+sia`` is also set.
232 * - :opt:`-tauc`
233 - Directly set the till yield stress `\tau_c`, in units Pa, at all grounded locations
234 and all times. Only effective if used with ``-yield_stress constant``, because
235 otherwise `\tau_c` is updated dynamically.
236
237 We find that an effective, though heuristic, way to determine `\phi` in
238 :eq:`eq-mohrcoulomb` is to make it a function of bed elevation
239 :cite:`AschwandenAdalgeirsdottirKhroulev`, :cite:`Martinetal2011`,
240 :cite:`Winkelmannetal2011`. This heuristic is motivated by hypothesis that basal material
241 with a marine history should be weak :cite:`HuybrechtsdeWolde`. PISM has a mechanism
242 setting `\phi` to be a *piecewise-linear* function of bed elevation. The
243 option is
244
245 .. code-block:: none
246
247 -topg_to_phi phimin,phimax,bmin,bmax
248
249 Thus the user supplies 4 parameters: `\phimin`, `\phimax`, `\bmin`, `\bmax`, where `b`
250 stands for the bed elevation. To explain these, we define `M = (\phimax - \phimin) /
251 (\bmax - \bmin)`. Then
252
253 .. math::
254 :label: eq-phipiecewise
255
256 \phi(x,y) =
257 \begin{cases}
258 \phimin, & b(x,y) \le \bmin, \\
259 \phimin + (b(x,y) - \bmin) \,M, & \bmin < b(x,y) < \bmax, \\
260 \phimax, & \bmax \le b(x,y).
261 \end{cases}
262
263 It is worth noting that an earth deformation model (see section :ref:`sec-beddef`) changes
264 `b(x,y)` (NetCDF variable :var:`topg`) used in :eq:`eq-phipiecewise`, so that a sequence
265 of runs such as
266
267 .. code-block:: none
268
269 pismr -i foo.nc -bed_def lc -stress_balance ssa+sia -topg_to_phi 10,30,-50,0 ... -o bar.nc
270 pismr -i bar.nc -bed_def lc -stress_balance ssa+sia -topg_to_phi 10,30,-50,0 ... -o baz.nc
271
272 will use *different* :var:`tillphi` fields in the first and second runs. PISM will print a
273 warning during initialization of the second run:
274
275 .. code-block:: none
276
277 * Initializing the default basal yield stress model...
278 option -topg_to_phi seen; creating tillphi map from bed elev ...
279 PISM WARNING: -topg_to_phi computation will override the 'tillphi' field
280 present in the input file 'bar.nc'!
281
282 Omitting :opt:`-topg_to_phi` in the second run would make PISM continue with the
283 same :var:`tillphi` field which was set in the first run.
284
285 .. _sec-effective-pressure:
286
287 Determining the effective pressure
288 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
289
290 When using the default option ``-yield_stress mohr_coulomb``, the effective pressure on
291 the till `N_{till}` is determined by the modeled amount of water in the till. Lower
292 effective pressure means that more of the weight of the ice is carried by the pressurized
293 water in the till and thus the ice can slide more easily. That is, equation
294 :eq:`eq-mohrcoulomb` sets the value of `\tau_c` proportionately to `N_{till}`. The amount
295 of water in the till is, however, a nontrivial output of the hydrology (see
296 :ref:`sec-subhydro`) and conservation-of-energy (see :ref:`sec-energy`) submodels in PISM.
297
298 Following :cite:`Tulaczyketal2000`, based on laboratory experiments with till extracted
299 from an ice stream in Antarctica, :cite:`BuelervanPelt2015` propose the following
300 parameterization which is used in PISM. It is based on the ratio
301 `s=W_{till}/W_{till}^{max}` where `W_{till}` is the effective thickness of water in the
302 till and `W_{till}^{max}` (:config:`hydrology.tillwat_max`) is the maximum amount of water
303 in the till (see section :ref:`sec-subhydro`):
304
305 .. math::
306 :label: eq-computeNtill
307
308 N_{till} = \min\left\{P_o, N_0 \left(\frac{\delta P_o}{N_0}\right)^s \, 10^{(e_0/C_c) \left(1 - s\right).}\right\}
309
310 Here `P_o` is the ice overburden pressure, which is determined entirely by the ice
311 thickness and density, and the remaining parameters are set by options in
312 :numref:`tab-effective-pressure`.
313
314 .. note::
315
316 While there is experimental support for the default values of `C_c`, `e_0`, and `N_0`,
317 the value of `\delta` should be regarded as uncertain, important, and subject to
318 parameter studies to assess its effect.
319
320 .. list-table:: Parameters controlling how till effective pressure `N_{till}` in equation
321 :eq:`eq-mohrcoulomb` is determined. All these have prefix
322 ``basal_yield_stress.mohr_coulomb.``.
323 :name: tab-effective-pressure
324 :header-rows: 1
325
326 * - Parameter
327 - Description
328 * - :config:`till_reference_void_ratio`
329 - `e_0` in :eq:`eq-computeNtill`, dimensionless, with default value 0.69
330 :cite:`Tulaczyketal2000`
331 * - :config:`till_compressibility_coefficient`
332 - `C_c` in :eq:`eq-computeNtill`, dimensionless, with default value 0.12
333 :cite:`Tulaczyketal2000`
334 * - :config:`till_effective_fraction_overburden`
335 - `\delta` in :eq:`eq-computeNtill`, dimensionless, with default value 0.02
336 :cite:`BuelervanPelt2015`
337 * - :config:`till_reference_effective_pressure`
338 - `N_0` in :eq:`eq-computeNtill`, in Pa, with default value 1000.0
339 :cite:`Tulaczyketal2000`
340
341 .. _sec-min-effective-pressure:
342
343 Controlling minimum effective pressure
344 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
345
346 The effective pressure `N_{till}` above satisfies (see equation 20 in
347 :cite:`BuelervanPelt2015`)
348
349 .. math::
350 :label: eq-mohr-coulomb-Ntill-bounds
351
352 \delta P_o \le N_{till} \le P_o.
353
354 In other words, `\delta` controls the lower bound of the effective pressure. In addition
355 to setting it using a configuration parameter (see :numref:`tab-effective-pressure`) one
356 can use a space- and time-dependent field. Set
357 :config:`basal_yield_stress.mohr_coulomb.delta.file` to the name of the file containing
358 the variable :var:`mohr_coulomb_delta` (dimensionless, `units=1`). Just like when using
359 other 2D time-dependent forcing, set
360 :config:`basal_yield_stress.mohr_coulomb.delta.period` and
361 :config:`basal_yield_stress.mohr_coulomb.delta.reference_year` to use periodic data. (See
362 :ref:`sec-periodic-forcing` for details.)
363
364 .. note::
365
366 PISM restricts the time step to capture changes in `\delta`. For example, if you
367 provide monthly records of `\delta` PISM will make sure no time step spans more than
368 one month.
369
370 PISM uses piecewise-linear interpolation in time for model times between records of
371 `\delta`.
372
373 ..
374 FIXME: EVOLVING CODE: If the :config:`basal_yield_stress.add_transportable_water`
375 configuration flag is set then the above formula becomes...