tvertchange.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
---
tvertchange.rst (7384B)
---
1 .. include:: ../global.txt
2
3 .. _sec-vertchange:
4
5 On the vertical coordinate in PISM, and a critical change of variable
6 =====================================================================
7
8 In PISM all fields in the ice, including enthalpy, age, velocity, and so on, evolve within
9 an ice fluid domain of *changing geometry*. See :numref:`fig-freebdry`. In
10 particular, the upper and lower surfaces of the ice fluid move with respect to the geoid.
11
12 .. figure:: figures/tempbdryfig.png
13 :name: fig-freebdry
14
15 The ice fluid domain evolves, with both the upper and lower surfaces in motion with
16 respect to the geoid.
17
18 .. note:: FIXME: This figure should show the floating case too.
19
20 The `(x,y,z)` coordinates in :numref:`fig-freebdry` are supposed to be from an orthogonal
21 coordinate system with `z` in the direction anti-parallel to gravity, so this is a
22 flat-earth approximation. In practice, the data inputs to PISM are in some particular
23 projection, of course.
24
25 We make a change of the independent variable `z` which simplifies how PISM deals
26 with the changing geometry of the ice, especially in the cases of a non-flat or moving
27 bed. We replace the vertical coordinate relative to the geoid with the vertical coordinate
28 relative to the base of the ice. Let
29
30 .. math::
31
32 s = \begin{cases}
33 z - b(x,y,t), & \text{ice is grounded}, \\
34 z - \frac{\rho_i}{\rho_o} H(x,y,t), & \text{ice is floating;}
35 \end{cases}
36
37 where `H = h - b` is the ice thickness and `\rho_i, \rho_o` are densities of
38 ice and ocean, respectively.
39
40 Now we make the change of variables
41
42 .. math::
43
44 (x,y,z,t) \mapsto (x,y,s,t)
45
46 throughout the PISM code. This replaces `z=b` by `s=0` as the equation of the
47 base surface of the ice. The ice fluid domain in the new coordinates only has a free upper
48 surface as shown in :numref:`fig-sfreebdry`.
49
50 .. figure:: figures/stempbdryfig.png
51 :name: fig-sfreebdry
52
53 In (x,y,s) space the ice fluid domain only has an upper surface which moves,
54 `s=H(x,y,t)`. Compare to :numref:`fig-freebdry`.
55
56 .. note:: FIXME: This figure should show the floating case too, and bedrock."
57
58 In PISM the computational domain (region)
59
60 .. math::
61
62 \mathcal{R}=\left\{(x,y,s)\big| -L_x\le x \le L_x, -L_y\le y \le L_y, -Lb_z \le s \le
63 L_z\right\}
64
65 is divided into a three-dimensional grid. See ``IceGrid``.
66
67 The change of variable `z\to s` used here *is not* the :cite:`Jenssen` change of variable
68 `\tilde s=(z-b)/H` . That change causes the conservation of energy equation to
69 become singular at the boundaries of the ice sheet. Specifically, the Jenssen change
70 replaces the vertical conduction term by a manifestly-singular term at ice sheet margins
71 where `H\to 0`, because
72
73 .. math::
74
75 \frac{\partial^2 E}{\partial z^2} = \frac{1}{H^2}
76 \frac{\partial^2 E}{\partial \tilde s^2}.
77
78 A singular coefficient of this type can be assumed to affect the stability of all
79 time-stepping schemes. The current change `s=z-b` has no such singularizing effect
80 though the change does result in added advection terms in the conservation of energy
81 equation, which we now address. See :ref:`sec-bombproof` for more general
82 considerations about the conservation of energy equation.
83
84 The new coordinates `(x,y,s)` are not orthogonal.
85
86 Recall that if `f=f(x,y,z,t)` is a function written in the old variables and if
87 `\tilde f(x,y,s,t)=f(x,y,z(x,y,s,t),t)` is the "same" function written in the new
88 variables, equivalently `f(x,y,z,t)=\tilde f(x,y,s(x,y,z,t),t)` , then
89
90 .. only:: html
91
92 .. When building PDFs this will be included already.
93
94 .. include:: ../math-definitions.txt
95
96 .. math::
97
98 \diff{f}{x} = \diff{\tilde f}{x} + \diff{\tilde f}{s} \diff{s}{x} = \diff{\tilde f}{x}
99 - \diff{\tilde f}{s} \diff{b}{x}.
100
101 Similarly,
102
103 .. math::
104
105 \diff{f}{y} = \diff{\tilde f}{y} - \diff{\tilde f}{s} \diff{b}{y},
106
107 .. math::
108
109 \diff{f}{t} = \diff{\tilde f}{t} - \diff{\tilde f}{s} \diff{b}{t}.
110
111 On the other hand,
112
113 .. math::
114
115 \diff{f}{z} = \diff{\tilde f}{s}.
116
117 The following table records some important changes to formulae related to conservation of
118 energy:
119
120 .. math::
121
122 \begin{array}{ll}
123 \textbf{old} & \textbf{new} \\
124 P=\rho g(h-z) & P=\rho g(H-s) \\
125 \diff{E}{t} & \diff{E}{t}-\diff{E}{s}\diff{b}{t} \\
126 \nabla E & \nabla E- \diff{E}{s}\nabla b \\
127 \rho_i\left(\diff{E}{t}+\mathbf{U}\cdot\nabla E + w\diff{E}{z}\right)=\frac{k_i}{c_i} \frac{\partial^2 E}{\partial z^2} + Q & \rho_i\left(\diff{E}{t} + \mathbf{U}\cdot\nabla E + \left(w-\diff{b}{t}-\mathbf{U}\cdot\nabla b\right)\diff{E}{s}\right) = \frac{k_i}{c_i} \frac{\partial^2 E}{\partial s^2} + Q
128 \end{array}
129
130 Note `E` is the ice enthalpy and `T` is the ice temperature (which is a
131 function of the enthalpy; see ``EnthalpyConverter``), `P` is the ice pressure
132 (assumed hydrostatic), `\mathbf{U}` is the depth-dependent horizontal velocity, and
133 `Q` is the strain-heating term.
134
135 Now the vertical velocity is computed by
136 ``StressBalance::compute_vertical_velocity(...)``. In the old coordinates
137 `(x,y,z,t)` it has this formula:
138
139 .. math::
140
141 w(z) = -\int_b^z \diff{u}{x}(z') + \diff{v}{y}(z')\,dz' + \diff{b}{t}
142 + \mathbf{U}_b \cdot \nabla b - S.
143
144 Here `S` is the basal melt rate, positive when ice is being melted. We have used the
145 basal kinematical equation and integrated the incompressibility statement
146
147 .. math::
148
149 \diff{u}{x} + \diff{v}{y} + \diff{w}{z} = 0.
150
151 In the new coordinates we have
152
153 .. math::
154
155 w(s) = -\int_0^s \diff{u}{x}(s') + \diff{v}{y}(s')\,ds'
156 + \mathbf{U}(s) \cdot \nabla b + \diff{b}{t} - S.
157
158 (Note that the term `\mathbf{U}(s) \cdot \nabla b` evaluates the horizontal velocity
159 at level `s` and not at the base.)
160
161 Let
162
163 .. math::
164
165 \tilde w(x,y,s,t) = w(s) - \diff{b}{t}-\mathbf{U}(s)\cdot\nabla b.
166
167 This quantity is the vertical velocity of the ice *relative to the location on the bed
168 immediately below it*. In particular, `\tilde w=0` for a slab sliding down a
169 non-moving inclined plane at constant horizontal velocity, if there is no basal melt rate.
170 Also, `\tilde w(s=0)` is nonzero only if there is basal melting or freeze-on, i.e.
171 when `S\ne 0`. Within PISM, `\tilde w` is written with name `wvel_rel` into an
172 input file. Comparing the last two equations, we see how
173 ``StressBalance::compute_vertical_velocity(...)`` computes `\tilde w` :
174
175 .. math::
176
177 \tilde w(s) = -\int_0^s \diff{u}{x}(s') + \diff{v}{y}(s')\,ds' - S.
178
179 The conservation of energy equation is now, in the new coordinate `s` and
180 newly-defined relative vertical velocity,
181
182 .. math::
183
184 \rho_i \left(\diff{E}{t} + \mathbf{U}\cdot\nabla E + \tilde w \diff{E}{s}\right)
185 = \frac{k_i}{c_i} \frac{\partial^2 E}{\partial s^2} + Q.
186
187 Thus it looks just like the conservation of energy equation in the original vertical
188 velocity `z`. This is the form of the equation solved by ``EnthalpyModel`` using
189 ``enthSystemCtx::solve()``.
190
191 Under option ``-o_size big``, all of these vertical velocity fields are available as
192 fields in the output NetCDF file. The vertical velocity relative to the geoid, as a
193 three-dimensional field, is written as the diagnostic variable ``wvel``. This is the
194 "actual" vertical velocity `w = \tilde w + \diff{b}{t} + \mathbf{U}(s)\cdot\nabla b`
195 . Its surface value is written as ``wvelsurf``, and its basal value as ``wvelbase``. The
196 relative vertical velocity `\tilde w` is written to the NetCDF output file as
197 ``wvel_rel``.