tbed_roughness.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
---
tbed_roughness.rst (13749B)
---
1 .. include:: ../global.txt
2
3 .. only:: html
4
5 .. in LaTeX I can just \usepackage{txfonts}...
6
7 `\require{cancel} \newcommand{\fint}{\cancel{\phantom{x}}\kern{-1.1em}\int}`
8
9 .. _sec-bed-roughness:
10
11 Using Schoof's parameterized bed roughness technique in PISM
12 ============================================================
13
14 .. contents::
15
16 An explanation
17 --------------
18
19 If the bed elevation ``topg`` is smoothed by preprocessing then we observe a reduction in
20 the peak values of the SIA diffusivities. From such smoothing there is (generically) also
21 a reduction in the peak magnitudes of horizontal velocities from both the SIA and SSA
22 models. The major consequence of these reductions, through the adaptive time-stepping
23 mechanism, is that PISM can take longer time steps and thus that it can complete model
24 runs in shorter time.
25
26 Large peak diffusivities coming from bed roughness are located (generically) at margin
27 locations where the ice is on, or has flowed onto, fjord-like bed topography. At coarser
28 resolutions (e.g. 20km and up), it appears that the effect of increasing bed roughness is
29 not as severe as at finer resolutions (e.g. 10km, 5km and finer). Of course it is true
30 that the shallow models we use, namely the SIA and SSA models, are missing significant
31 stress gradients at the same margin locations which have large bed slopes.
32
33 Here we are emphasizing the performance "hit" which the whole model experiences if some
34 small part of the ice sheet is on a rough bed. That part therefore is not well-modeled
35 anyway, compared to the majority of the ice sheet. (Switching to full Stokes or Blatter
36 higher order models without major spatial adaptivity would probably imply a gain in the
37 balanced stress components *and* a loss of the ability to model the ice sheet at high
38 resolution. There is a tradeoff between completeness of the continuum model and usable
39 resolution needed to resolve the features of the real ice sheet.)
40
41 There exists a theory which addresses exactly this situation for the SIA model, and
42 explains rigorously that one should use a smoothed bed in that model, but with an
43 associated reduction in diffusivity. This theory explains how to improve the SIA model to
44 handle bed roughness more correctly, because it parameterizes the effects of
45 "higher-order" stresses which act on the ice as it flows over bed topography. Specifically
46 it shows the way to a double performance boost for PISM:
47
48 - smoothed beds give longer time steps directly, and
49 - the parameterized effect of the local bed roughness is to further reduce the
50 diffusivity, giving even longer time-steps.
51
52
53 Theory
54 ------
55
56 The theory is in Christian Schoof's (2003) *The effect of basal topography on ice sheet
57 dynamics* :cite:`Schoofbasaltopg2003`. His mathematical technique is to expand the
58 Stokes equations in two levels of horizontal scales, one for the entire ice sheet (denoted
59 `[L]`) and one for the horizontal scale (wavelength) of bed topography (`[S]`).
60 The "inner" scaling assumes that the typical ice sheet thickness `[D]` is small
61 compared to `[S]`, while the "outer" scaling assumes that `[S]` is small compared
62 to `[L]`:
63
64 .. math::
65
66 \nu = \frac{[D]}{[S]} \ll 1, \qquad \delta = \frac{[S]}{[L]} \ll 1.
67
68 Specifically, there is an "inner" horizontal variable `x` describing the local
69 topography on scales comparable to `[S]` or smaller, and an "outer" horizontal
70 variable `X` describing the large scale bed topography and ice sheet flow on scales
71 larger than `[S]`.
72
73 In order to describe the Schoof scheme using PISM notation, we start by recalling the mass
74 continuity equation which is fundamental to any shallow ice theory:
75
76 .. math::
77
78 \frac{\partial H}{\partial t} = (M - S) - \nabla\cdot \mathbf{q}.
79
80 Within PISM this equation is handled by GeometryEvolution. Recall that `M-S` is the mass
81 balance added to the ice column per time. (It plays no further role here.) In the SIA case
82 with zero basal sliding, the horizontal mass flux is
83
84 .. math::
85
86 \mathbf{q} = - D_{SIA} \nabla h,
87
88 where `D_{SIA}\ge 0` is given next. Thus the mass continuity equation is *diffusive*. The
89 diffusivity `D_{SIA}` is a function of the ice geometry and the ice flow law. In the
90 isothermal Glen power law (power `= n`) case we recall
91
92 .. math::
93 :label: eq-siadiffusivity
94
95 D_{SIA} = \Gamma H^{n+2} |\nabla h|^{n-1}
96
97 where `\Gamma = 2 A (\rho g)^n / (n+2)` (c.f. details in :cite:`BLKCB`\).
98
99 Consider now the "original" bed topography `b_0(x_1,x_2)`, which we assume for the moment
100 is independent of time. (Time-independence is not actually critical, and such a
101 restriction can be removed.) We will use `x_1,x_2` to denote the horizontal model
102 coordinates, though they are denoted `x,y` elsewhere in these PISM docs. Suppose a
103 locally-smoothed bed is computed by averaging `b_0` over a rectangular region of sides
104 `2\lambda_1` by `2\lambda_2`, namely:
105
106 .. math::
107
108 b_s(x_1,x_2) = \fint b_0(x_1+\xi_1,x_2+\xi_2)\,d\xi_1\,d\xi_2
109
110 where the slashed integral symbol is defined as
111
112 .. math::
113
114 \fint F(\xi_1,\xi_2)\,d\xi_1\,d\xi_2 = \frac{1}{(2\lambda_1)(2\lambda_2)} \int_{-\lambda_1}^{\lambda_1} \int_{-\lambda_2}^{\lambda_2} F(\xi_1,\xi_2)\,d\xi_1\,d\xi_2.
115
116 Consider also the "local bed topography"
117
118 .. math::
119
120 \tilde b(x_1,x_2,\xi_1,\xi_2) = b_0(x_1+\xi_1,x_2+\xi_2) - b_s(x_1,x_2).
121
122 As a function of the local coordinates `\xi_1,\xi_2`, the local bed topography `\tilde b`
123 is the amount by which the bed deviates from the "local average" `b_s(x_1,x_2)`. Generally
124 we will use `-\lambda_1 \le \xi_1 \le \lambda_1`, `-\lambda_2 \le \xi_2 \le \lambda_2` as
125 the smoothing domain, but these specific ranges are not required by the formulas above.
126 Note that the average of the local bed topgraphy is zero by definition:
127
128 .. math::
129
130 \fint \tilde b(x_1,x_2,\xi_1,\xi_2)\,d\xi_1\,d\xi_2 = 0.
131
132 The result of Schoof's scaling arguments (:cite:`Schoofbasaltopg2003`, equation (49)) is to
133 modify the diffusivity by multiplying by a factor `0 \le \theta \le 1`:
134
135 .. math::
136
137 D = \theta(h(x_1,x_2),x_1,x_2) \, D_{SIA}.
138
139 where `D_{SIA}` is defined by :eq:`eq-siadiffusivity` earlier, and
140
141 .. math::
142 :label: eq-thetadefn
143
144 \theta(h,x_1,x_2) = \left[ \fint \left(1 - \frac{\tilde b(x_1,x_2,\xi_1,\xi_2)}{H}\right)^{-(n+2)/n}\,d\xi_1\,d\xi_2 \right]^{-n}
145
146 Here the ice thickness and ice surface elevation are related to the *smoothed* bed
147 topography, so that in PISM notation
148
149 .. math::
150
151 H(t,x_1,x_2) = h(t,x_1,x_2) - b_s(x_1,x_2).
152
153 This can be treated as the *definition* of the ice thickness `H` in the above formula for
154 `\theta`.
155
156 The formula for `\theta` has additional terms if there is basal sliding, but we consider
157 only the non-sliding SIA here.
158
159 The very important fact that `0 \le \theta \le 1` is proven in appendix A of
160 :cite:`Schoofbasaltopg2003` by a Jensen's inequality argument. (See also the convexity
161 argument at the bottom of this page.)
162
163
164 Practical application, and Taylor approximation
165 -----------------------------------------------
166
167 The above formulas already reflect the recommendations Schoof gives on how to apply his
168 formulas (:cite:`Schoofbasaltopg2003`, subsection 4.2). The rest of this page is devoted to
169 how the class ``stressbalance::BedSmoother`` implements a practical version of this
170 theory, based on these recommendations plus some additional approximation.
171
172 The averages appearing in his scaling arguments are over an infinite domain, e.g.
173
174 .. math::
175
176 f_s(x) = \lim_{R\to\infty} \frac{1}{2R} \int_{-R}^R f(\xi,x)\,d\xi.
177
178 For practical modeling use, Schoof specifically recommends averaging over some finite
179 length scale which should be "considerably greater than the grid spacing, but much smaller
180 than the size of the ice sheet." Furthermore he recommends that, because of the typical
181 aspect ratio of ice sheets, "Bed topography on much larger length scales than 10 km should
182 then be resolved explicitly through the smoothed bed height `b_s` rather than the
183 correction factor `\theta`." Thus in PISM we use `\lambda_1 = \lambda_2 = 5` km as the
184 default (set :config:`stress_balance.sia.bed_smoother.range` to change this value).
185
186 It is, of course, possible to have bed roughness of significant magnitude at essentially
187 any wavelength. We make no claim that PISM results are good models of ice flow over
188 *arbitrary* geometry; clearly the current models cannot come close to the non-shallow
189 solution (Stokes) in such cases. Rather, the goal right now is to improve on the existing
190 shallow models, the diffusive SIA specifically, while maintaining or increasing
191 high-resolution performance and comprehensive model quality, which necessarily includes
192 many other modeled physical processes like ice thermal state, basal lubrication, and so
193 on. The desirable properties of the Schoof scheme are accepted not because the resulting
194 model is perfect, but because we gain both a physical modeling improvement *and* a
195 computational performance improvement from its use.
196
197 How do we actually compute expression :eq:`eq-thetadefn` quickly? Schoof has this suggestion,
198 which we follow: "To find `\theta` for values of [surface elevation for which `\theta` has
199 not already been computed], some interpolation scheme should then be used. `\theta` is
200 then represented at each grid point by some locally-defined interpolating function [of the
201 surface elevation]."
202
203 We need a "locally-defined interpolating function". As with any approximation scheme,
204 higher accuracy is achieved by doing "more work", which here is an increase in memory used
205 for storing spatially-dependent coefficients. Pade rational approximations, for example,
206 were considered but are excluded because of the appearance of uncontrolled poles in the
207 domain of approximation. The 4th degree Taylor polynomial is chosen here because it shares
208 the same convexity as the rational function it approximates; this is proven below.
209
210 Use of Taylor polynomial `P_4(x)` only requires the storage of three fields (below),
211 but it has been demonstrated to be reasonably accurate by trying beds of various
212 roughnesses in Matlab/Octave scripts. The derivation of the Taylor polynomial is most
213 easily understood by starting with an abstract rational function like the one which
214 appears in :eq:`eq-thetadefn`, as follows.
215
216 The fourth-order Taylor polynomial of the function `F(s)=(1-s)^{-k}` around `s=0` is
217
218 .. math::
219
220 P_4(s) = 1 + k s + \frac{k(k+1)}{2} s^2 + \dots + \frac{k(k+1)(k+2)(k+3)}{4!} s^4,
221
222 so `F(s) = P_4(s) + O(s^5)`. Let
223
224 .. math::
225
226 \omega(f,z) = \fint (1 - f(\xi) z)^{-k}\,d\xi
227
228 where `f` is some function and `z` a scalar. Then
229
230 .. math::
231
232 \omega(f,z) &= \fint (1 - f(\xi) z)^{-k}\,d\xi = \fint P_4(f(\xi) z)\,d\xi + O((\max |f(\xi)|)^5 |z|^5)
233
234 &\approx 1 + k z \fint f(\xi)\,d\xi + \frac{k(k+1)}{2} z^2 \fint f(\xi)^2\,d\xi + \dots + \frac{k(k+1)(k+2)(k+3)}{4!} z^5 \fint f(\xi)^4\,d\xi.
235
236
237 Now, `\theta` can be written
238
239 .. math::
240
241 \theta(h,x_1,x_2) = \left[ \omega(\tilde b(x_1,x_2,\cdot,\cdot), H^{-1}) \right]^{-n}.
242
243 So our strategy should be clear, to approximate `\omega(\tilde b(x_1,x2,\cdot,\cdot),
244 H^{-1})` by the Taylor polynomial as a function of `H^{-1}`, whose the coefficients depend
245 on `x_1,x_2`. We thereby get a rapidly-computable approximation for `\theta` using stored
246 coefficients which depend on `x_1,x_2`. In fact, let `f(\xi) = \tilde
247 b(x_1,x_2,\xi_1,\xi_2)` for fixed `x_1,x_2`, and let `z=H^{-1}`. Recall that the mean of
248 this `f(\xi)` is zero, so that the first-order term drops out in the above expansion of
249 `\omega`. We have the following approximation of `\theta`:
250
251 .. math::
252 :label: eq-thetaapprox
253
254 \theta(h,x_1,x_2) \approx \left[ 1 + C_2(x_1,x_2) H^{-2} + C_3(x_1,x_2) H^{-3} + C_4(x_1,x_2) H^{-4} \right]^{-n}
255
256 where
257
258 .. math::
259
260 C_q(x_1,x_2) = \frac{k(k+1)\dots(k+q-1)}{q!} \fint \tilde b(x_1,x_2,\xi_1,\xi_2)^q\,d\xi_1\,d\xi_2
261
262 for `q=2,3,4` and `k = (n+2)/n`.
263
264 Note that the coefficients `C_q` depend only on the bed topography, and not on the ice
265 geometry. Thus we will pre-process the original bed elevations `b_0` to compute and store
266 the fields `b_s,C_2,C_3,C_4`. The computation of equation :eq:`eq-thetaapprox` is fast and
267 easily-parallelized if these fields are pre-stored. The computation of the coefficients
268 `C_q` and the smoothed bed `b_s` at the pre-processing stage is more involved, especially
269 when done in parallel.
270
271 The parameters `\lambda_1,\lambda_2` must be set, but as noted above we use a default
272 value of 5 km based on Schoof's recommendation. This physical distance may be less than or
273 more than the grid spacing. In the case that the grid spacing is 1 km, for example, we see
274 that there is a large smoothing domain in terms of the number of grid points. Generally,
275 the ghosting width (in PETSc sense) is unbounded. We therefore move the unsmoothed
276 topography to processor zero and do the smoothing and the coefficient-computing there. The
277 class ``stressbalance::BedSmoother`` implements these details.
278
279 Convexity of `P_4`
280 ^^^^^^^^^^^^^^^^^^
281
282 The approximation :eq:`eq-thetaapprox` given above relates to the Jensen's inequality argument
283 used by Schoof in Appendix A of :cite:`Schoofbasaltopg2003`. For the nonsliding case, his
284 argument shows that because `F(s)=(1-s)^{-k}` is convex on `[-1,1]` for `k>0`, therefore
285 `0 \le \theta \le 1`.
286
287 Thus it is desirable for the approximation `P_4(z)` to be convex on the same interval,
288 and this is true. In fact,
289
290 .. math::
291
292 P_4''(s) = k(k+1) \left[1 + (k+2) s + \frac{(k+2)(k+3)}{2} s^2\right],
293
294 and this function turns out to be positive for all `s`. In fact we will show that the
295 minimum of `P_4''(s)` is positive. That minimum occurs where `P_4'''(s)=0` or
296 `s=s_{min}=-1/(k+3)`. It is a minimum because `P_4^{(4)}(s)` is a positive
297 constant. And
298
299 .. math::
300
301 P_4''(s_{min}) = \frac{k(k+1)(k+4)}{2(k+3)} > 0.