tsteady-hydrology.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
---
tsteady-hydrology.rst (7076B)
---
1 .. include:: ../global.txt
2 .. include:: ../math-definitions.txt
3
4 .. _sec-steady-hydro:
5
6 Computing steady-state subglacial water flux
7 ============================================
8
9 The "routing" subglacial hydrology model is described by equations
10
11 .. math::
12 :label: eq-hydrology-routing
13
14 \diff{W}{t} + \diff{\Wtill}{t} + \div \bq = \frac{m}{\rho_w}
15
16 \diff{\Wtill}{t} = \frac{m}{\rho_w} - C_d
17
18 \bq = -k W^{\alpha} |\nabla \psi|^{\beta - 2} \nabla \psi
19
20 \psi = P_o + \rho_w g (b + W)
21
22 on a an ice covered area `\Omega`. We assume zero flux boundary conditions on the inflow
23 part of the boundary and no boundary condition on the outflow boundary. See
24 :cite:`BuelervanPelt2015` (equations 1, 2, 6, 16, 26) for details and the notation. Here
25 we also assume that `m \ge 0`.
26
27 Our goal is to estimate `Q = \bq \cdot \n`, the flux through the outflow part of the
28 boundary of `\Omega` corresponding to the steady state of :eq:`eq-hydrology-routing` using
29 a method that is computationally cheaper than using the explicit in time approximation
30 of :eq:`eq-hydrology-routing` described by :cite:`BuelervanPelt2015`.
31
32 Pick a contiguous section `\omega` of `\partial \Omega` (the terminus of an outlet
33 glacier, for example). Let `B` be the union of all the trajectories of the vector field
34 `\bq` in `\Omega` that pass through `\omega`. The area `B` is the "drainage basin"
35 corresponding to `\omega`.
36
37 Let `\gamma = \partial B \setminus \omega`. Note that if a point `P` is in `\gamma` then
38 one of the following conditions is satisfied.
39
40 #. `|\bq| = 0` (it is the origin of a trajectory that starts within `\Omega`) or
41 #. `P \in \partial \Omega` (specifically, `P` is a part of the inflow part of the boundary
42 of `\Omega`)
43 #. `\bq \cdot \n = 0` (`P` is not at the end of a trajectory, and so the normal to the
44 boundary is orthogonal to `\bq`).
45
46 Therefore `\bq \cdot \n = 0` on `\gamma` and
47
48 .. math::
49
50 \oint_{\partial B} \bq \cdot \n\; ds &= \int_{\omega} \bq \cdot \n\; ds + \int_{\gamma} \bq \cdot \n\; ds
51
52 &= \int_{\omega} \bq \cdot \n ds.
53
54 Assuming the steady state (and setting time derivatives in :eq:`eq-hydrology-routing` to
55 zero), integrating over `B`, and applying the divergence theorem gives
56
57 .. math::
58 :label: eq-steady-hydro-1
59
60 \int_{\omega} \bq \cdot \n\; ds = \int_{B} \frac{m}{\rho_w},
61
62 i.e. *in a steady state the flux through a terminus is equal to the total rate at which
63 water is added to the corresponding drainage basin due to the source term*.
64
65 Next, consider a related initial boundary value problem
66
67 .. math::
68 :label: eq-emptying-problem
69
70 \diff{u}{t} = -\div (\V u)
71
72 on `B` with `u(x, y, 0) = u_0(x, y)` (`u_0 \ge 0`), `\V = -k(x, y) \nabla \psi`, zero flux
73 on the inflow boundary, and no boundary condition on the outflow boundary.
74
75 Here `\psi` is the hydraulic potential corresponding to the steady state of
76 :eq:`eq-hydrology-routing` and `k(x, y)` is a strictly positive but otherwise arbitrary
77 conductivity function.
78
79 Note that since `\psi` is a steady state hydraulic potential all trajectories of the
80 vector field `\V` leave `B` and for `\epsilon > 0` there is a time `T > 0` such that
81
82 .. math::
83
84 \int_B u(T) = \epsilon \int_B u_0.
85
86 Integrating over time from `0` to `T`, we get
87
88 .. math::
89
90 \int_0^T \diff{u}{t}\, dt = - \int_0^T \div (\V u),\, \text{or}
91
92 .. math::
93
94 u_0 = u(T) + \int_0^T \div (\V u).
95
96 Integrating over `B` and using the divergence theorem gives
97
98 .. math::
99
100 \int_B u_0 &= \int_B u(T) + \int_B \int_0^T \div (\V u)
101
102 &= \epsilon \int_B u_0 + \int_0^T \int_B \div (\V u)
103
104 &= \epsilon \int_B u_0 + \int_0^T \oint_{\partial B} (\V u) \cdot \n
105
106 &= \epsilon \int_B u_0 + \int_0^T \int_{\omega} (\V u) \cdot \n.
107
108 Finally,
109
110 .. math::
111 :label: eq-steady-hydro-2
112
113 \int_B u_0 = \frac{1}{1 - \epsilon} \int_0^T \int_{\omega} (\V u) \cdot \n
114
115 Combining :eq:`eq-steady-hydro-1` and :eq:`eq-steady-hydro-2` and choosing `u_0 = \tau m\,
116 /\, \rho_w` for some `\tau > 0`\ [#]_ gives us a way to estimate the flux through the
117 outflow boundary *if we know the direction of the steady state flux*:
118
119 .. math::
120 :label: eq-steady-hydro-3
121
122 \int_{\omega} \bq \cdot \n\; ds = \frac{1}{\tau(1 - \epsilon)} \int_0^T \int_{\omega} (\V u) \cdot \n\; ds.
123
124 Here the right hand side of :eq:`eq-steady-hydro-3` can be estimated by advancing an
125 explicit-in-time approximation of :eq:`eq-emptying-problem` until `\int_B u` drops below
126 a chosen threshold.
127
128 However, the direction of the steady state flux `\bq` depends on steady state
129 distributions of `W` and `\Wtill` and these quantities are expensive to compute.
130
131 To avoid this issue we note that `W \ll H` and so `\psi` is well approximated by `\psi_0 =
132 P_o + \rho_w g b` everywhere except the vicinity of subglacial lakes. Moreover, if
133 `|\nabla W|` is small then `\nabla \psi_0` is a reasonable approximation of `\nabla \psi`.
134
135 We approximate `\psi` by `\tilde \psi = P_o + \rho_w g b + \delta` where `\delta > 0` is
136 an adjustment needed to ensure that `\tilde \psi` has no local minima in the interior of
137 `\Omega` and `|\nabla \tilde \psi| > 0` everywhere on `\Omega` except possibly on a set of
138 measure zero (no "plateaus").
139
140 The approximation of `\tilde \psi` on a computational grid is computed as follows.
141
142 1. Set `k = 0`, `\tilde \psi_{0} = \psi`.
143 2. Iterate over all grid points. If a grid point `(i, j)` is at a local minimum, set
144 `\tilde \psi_{k + 1}(i, j)` to the average of neighboring values of `\tilde \psi_{k}`
145 plus a small increment `\Delta \psi`, otherwise set `\tilde \psi_{k + 1}(i, j)` to
146 `\tilde \psi_{k}(i, j)`.
147 3. If step 2 found no local minima, stop. Otherwise increment `k` and proceed to step 2.
148
149 Next, note that it is not necessary to identify the drainage basin `B` for a terminus
150 `\omega`: it is defined by `\psi` and therefore an approximation of
151 :eq:`eq-emptying-problem` will automatically distribute water inputs from the ice surface
152 (or melting) along the ice margin.
153
154 .. _sec-steady-hydro-algorithm:
155
156 The algorithm
157 ^^^^^^^^^^^^^
158
159 Using an explicit time stepping approximation of :eq:`eq-emptying-problem` we can estimate
160 `\int_{\omega} \bq \cdot \n \; ds` as follows.
161
162 #. Given ice thickness `H` and bed elevation `b` compute `\tilde \psi` by filling "dips"
163 as described above.
164
165 #. Choose the stopping criterion `\epsilon > 0` and the scaling for the source term `\tau > 0`.
166
167 #. Set
168
169 .. math::
170
171 u &\leftarrow \frac{\tau m}{\rho_w},
172
173 t &\leftarrow 0,
174
175 Q &\leftarrow (0, 0).
176
177 #. Compute the CFL time step `\Delta t` using `u` and `\V`.
178
179 #. Perform an explicit step from `t` to `t + \Delta t`, updating `u`.
180
181 #. Accumulate this step's contribution to `Q`:
182
183 .. math::
184
185 Q \leftarrow Q + \Delta t \cdot \V u.
186
187 #. Set `t \leftarrow t + \Delta t`
188
189 #. If `\int_{\Omega} u\; dx\, dy > \epsilon`, go to 4.
190
191 #. Set
192
193 .. math::
194
195 Q \leftarrow \frac{1}{t (1 - \epsilon^{*})}\; Q,
196
197 where
198
199 .. math::
200
201 \epsilon^{*} = \frac{\int_{\Omega} u}{\int_{\Omega} u_{0}}.
202
203 .. rubric:: Footnotes
204
205 .. [#] The constant `\tau` is needed to get appropriate units, but its value is irrelevant.