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.