tstress-balance.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
---
tstress-balance.rst (10667B)
---
1 .. include:: ../../../global.txt
2
3 .. _sec-stressbalance:
4
5 Choosing the stress balance
6 ---------------------------
7
8 The basic stress balance used for all grounded ice in PISM is the non-sliding,
9 thermomechanically-coupled SIA :cite:`BBL`. For the vast majority of most ice sheets, as
10 measured by area or volume, this is an appropriate model, which is an `O(\epsilon^2)`
11 approximation to the Stokes model if `\epsilon` is the depth-to-length ratio of the ice
12 sheet :cite:`Fowler`.
13
14 The shallow shelf approximation (SSA) stress balance applies to floating ice. See the Ross
15 ice shelf example in section :ref:`sec-ross` for an example in which the SSA is only
16 applied to floating ice.
17
18 In PISM the SSA is also used to describe the sliding of grounded ice and the formation of
19 ice streams :cite:`BBssasliding`. Specifically for the SSA with "plastic" (Coulomb
20 friction) basal resistance, the locations of ice streams are determined as part of a free
21 boundary problem of Schoof :cite:`SchoofStream`, a model for emergent ice streams within a
22 ice sheet and ice shelf system. This model explains ice streams through a combination of
23 plastic till failure and SSA stress balance.
24
25 This SSA description of ice streams is the preferred "sliding law" for the SIA
26 :cite:`BBssasliding`, :cite:`Winkelmannetal2011`. The SSA should be combined with the SIA,
27 in this way, in preference to classical SIA sliding laws which make the sliding velocity of ice a
28 local function of the basal value of the driving stress. The resulting combination of SIA
29 and SSA is a "hybrid" approximation of the Stokes model :cite:`Winkelmannetal2011`. Option
30 ``-stress_balance ssa+sia`` turns on this "hybrid" model. In this use of the SSA as a
31 sliding law, floating ice is also subject to the SSA.
32
33 :numref:`tab-stress-balance-choice` describes the basic choice of stress balance.
34
35 .. list-table:: The basic choice of stress balance
36 :name: tab-stress-balance-choice
37 :header-rows: 1
38 :widths: 1,2
39
40 * - Option
41 - Description
42
43 * - :opt:`-stress_balance none`
44 - Turn off ice flow completely.
45
46 * - :opt:`-stress_balance sia` (default)
47 - Grounded ice flows by the non-sliding SIA. Floating ice does not flow, so this model
48 is not recommended for marine ice sheets.
49
50 * - :opt:`-stress_balance ssa`
51 - Use the SSA model exclusively. Horizontal ice velocity is constant throughout ice
52 columns.
53
54 * - :opt:`-stress_balance prescribed_sliding`
55 - Use the constant-in-time prescribed sliding velocity field read from a file set
56 using :opt:`-prescribed_sliding_file`, variables ``ubar`` and ``vbar``.
57 Horizontal ice velocity is constant throughout ice columns.
58
59 * - :opt:`-stress_balance ssa+sia`
60 - The recommended sliding law, which gives the SIA+SSA hybrid stress balance.
61 Combines SSA-computed velocity, using pseudo-plastic till, with SIA-computed
62 velocity according to the combination in :cite:`Winkelmannetal2011`; similar to
63 :cite:`BBssasliding`. Floating ice uses SSA only.
64
65 * - :opt:`-stress_balance prescribed_sliding+sia`
66 - Use the constant-in-time prescribed sliding velocity in combination with the
67 non-sliding SIA.
68
69 .. _sec-ssa:
70
71 Controlling the SSA stress balance model
72 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
73
74 If the SSA stress balance is used, a choice of two solvers is available, namely
75 ``-ssa_method fd`` (default) or ``-ssa_method fem``. See :numref:`tab-ssa-usage`, which
76 describes additional controls on the numerical solution of the stress balance equations.
77 If option ``-ssa_method fd`` is chosen then several more controls on numerics are
78 available; see :numref:`tab-ssafd-controls`. If the ice sheet being modeled has any
79 floating ice then the user is advised to read section :ref:`sec-pism-pik` on modeling
80 marine ice sheets.
81
82 When using SSA as a "sliding law" one also needs to model the yield stress, or a
83 pseudo-yield-stress in the case of power law sliding (section :ref:`sec-basestrength`).
84
85 The basal yield stress is normally a function of the amount of water stored in the till
86 and a (generally) spatially-varying till strength. The amount of stored basal water is
87 modeled by the subglacial hydrology model (section :ref:`sec-subhydro`) based on the basal
88 melt rate which is, primarily, thermodynamically-determined (see :ref:`sec-energy`).
89
90
91 .. list-table:: Choice of, and controls on, the numerical SSA stress balance.
92 :name: tab-ssa-usage
93 :header-rows: 1
94 :widths: 1,2
95
96 * - Option
97 - Description
98
99 * - :opt:`-ssa_method` [ ``fd | fem`` ]
100 - Both finite difference (``fd``; the default) and finite element (``fem``) versions
101 of the SSA numerical solver are implemented in PISM. The ``fd`` solver is the only
102 one which allows PIK options (section :ref:`sec-pism-pik`). ``fd`` uses Picard
103 iteration :cite:`BBssasliding`, while ``fem`` uses a Newton method. The ``fem`` solver
104 has surface velocity inversion capability :cite:`Habermannetal2013`.
105
106 * - :opt:`-ssa_eps` (`10^{13}`)
107 - The numerical schemes for the SSA compute an effective viscosity `\nu` which
108 depends on strain rates and ice hardness (thus temperature). The minimum value of
109 the effective viscosity times the thickness (i.e. `\nu H`) largely determines the
110 difficulty of solving the numerical SSA. This constant is added to keep `\nu H`
111 bounded away from zero: `\nu H \to \nu H + \epsilon_{\text{SSA}}`, where
112 `\epsilon_{\text{SSA}}` is set using this option. Units of :opt:`ssa_eps` are
113 `\text{Pa}\,\text{m}\,\text{s}`. Set to zero to turn off this lower bound.
114
115 * - :opt:`-ssa_view_nuh`
116 - View the product `\nu H` for your simulation as a runtime viewer (section
117 :ref:`sec-diagnostic-viewers`). In a typical Greenland run we see a wide range of
118 values for `\nu H` from `\sim 10^{14}` to `\sim 10^{20}`
119 `\text{Pa}\,\text{m}\,\text{s}`.
120
121 .. list-table:: Controls on the numerical iteration of the ``-ssa_method fd`` solver
122 :name: tab-ssafd-controls
123 :header-rows: 1
124 :widths: 1,2
125
126 * - Option
127 - Description
128
129 * - :opt:`-ssafd_picard_maxi` (300)
130 - Set the maximum allowed number of Picard (nonlinear) iterations in solving the
131 shallow shelf approximation.
132
133 * - :opt:`-ssafd_picard_rtol` (`10^{-4}`)
134 - The Picard iteration computes a vertically-averaged effective viscosity which is
135 used to solve the equations for horizontal velocity. Then the new velocities are
136 used to recompute an effective viscosity, and so on. This option sets the relative
137 change tolerance for the effective viscosity. The Picard iteration stops when
138 successive values `\nu^{(k)}` of the vertically-averaged effective viscosity
139 satisfy
140
141 .. math::
142
143 \|(\nu^{(k)} - \nu^{(k-1)}) H\|_1 \le Z \|\nu^{(k)} H\|_1
144
145 where `Z=` ``ssafd_picard_rtol``.
146
147 * - :opt:`-ssafd_ksp_rtol` (`10^{-5}`)
148 - Set the relative change tolerance for the iteration inside the Krylov linear solver
149 used at each Picard iteration.
150
151 * - :opt:`-ssafd_max_speed` (`50 km/yr`)
152 - Limits computed SSA velocities: ice speed is capped at this limit after each Picard
153 iteration of the SSAFD solver. This may allow PISM to take longer time steps by
154 ignoring high velocities at a few troublesome locations.
155
156 .. _sec-sia:
157
158 Controlling the SIA stress balance model
159 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
160
161 The section :ref:`sec-age` describes coupling SIA stress balance to the age of the ice.
162
163 The section :ref:`sec-gradient` covers available methods of computing the surface gradient.
164
165 .. note::
166
167 The explicit time stepping of the mass continuity equation in the case of the SIA flow
168 comes with a severe restriction on time step length:
169
170 .. math::
171
172 \Delta t \le \frac{2 R}{D\left( 1/\Delta x^2 + 1/\Delta y^2 \right)}
173
174 Here `D` is the maximum diffusivity of the SIA flow and `R` is
175 :config:`time_stepping.adaptive_ratio`, a tuning parameter that further reduces the
176 maximum allowed time step length.
177
178 The maximum diffusivity `D` may be achieved at an isolated grid point near the ice
179 margin. In this case it might make sense to limit the diffusivity of the SIA flow,
180 sacrificing accuracy at a few grid points to increase time step length and reduce the
181 computational cost. Set :config:`stress_balance.sia.limit_diffusivity` to enable this
182 mechanism.
183
184 When :config:`stress_balance.sia.limit_diffusivity` is ``false`` PISM stops as soon as
185 the SIA diffusivity at any grid point exceeds
186 :config:`stress_balance.sia.max_diffusivity`. We do this to make it easier to detect
187 problematic model configurations: in many cases it does not make sense to continue a
188 simulation if `D` is very large.
189
190 .. _sec-weertman:
191
192 Weertman-style sliding
193 ^^^^^^^^^^^^^^^^^^^^^^
194
195 .. warning::
196
197 This kind of sliding is, in general, a bad idea. We implement it to simplify
198 comparisons of the "hybrid" model mentioned above to older studies using this
199 parameterization.
200
201 PISM implements equation 5 from :cite:`Tomkin2007`:
202
203 .. math::
204 :label: eq-weertman-sliding
205
206 \mathbf{u}_s = \frac{2 A_s \beta_c (\rho g H)^{n}}{N - P} |\nabla h|^{n-1} \nabla h.
207
208 .. list-table:: Notation used in :eq:`eq-weertman-sliding`
209 :name: tab-weertman-notation
210 :header-rows: 1
211 :widths: 1,9
212
213 * - Variable
214 - Meaning
215
216 * - `H`
217 - ice thickness
218
219 * - `h`
220 - ice surface elevation
221
222 * - `n`
223 - flow law exponent
224
225 * - `g`
226 - acceleration due to gravity
227
228 * - `\rho`
229 - ice density
230
231 * - `N`
232 - ice overburden pressure, `N = \rho g H`
233
234 * - `P`
235 - basal water pressure
236
237 * - `A_s`
238 - sliding parameter
239
240 * - `\beta_c`
241 - "constriction parameter" capturing the effect of valley walls on the flow;
242 set to `1` in this implementation
243
244 We assume that the basal water pressure is a given constant fraction of the overburden
245 pressure: `P = k N`. This simplifies :eq:`eq-weertman-sliding` to
246
247 .. math::
248
249 \mathbf{u}_s = \frac{2 A_s}{1 - k} ( \rho g H\, |\nabla h| )^{n-1} \nabla h.
250
251 This parameterization is used for grounded ice *where the base of the ice is temperate*.
252
253 To enable, use :opt:`-stress_balance weertman_sliding` (this results in constant-in-depth
254 ice velocity) or :opt:`-stress_balance weertman_sliding+sia` to use this parameterization
255 as a sliding law with the deformational flow modeled using the SIA model.
256
257 Use configuration parameters :config:`stress_balance.weertman_sliding.k` and
258 :config:`stress_balance.weertman_sliding.A` tot set `k` and `A_s`, respectively. Default
259 values come from :cite:`Tomkin2007`.