thydrology.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
---
thydrology.rst (12323B)
---
1 .. include:: ../../../global.txt
2 .. include:: ../../../math-definitions.txt
3
4 .. _sec-subhydro:
5
6 Subglacial hydrology
7 --------------------
8
9 At the present time, two simple subglacial hydrology models are implemented *and
10 documented* in PISM, namely ``-hydrology null`` (and its modification ``steady``) and
11 ``-hydrology routing``; see :numref:`tab-hydrologychoice` and :cite:`BuelervanPelt2015`.
12 In both models, some of the water in the subglacial layer is stored locally in a layer of
13 subglacial till by the hydrology model. In the ``routing`` model water is conserved by
14 horizontally-transporting the excess water (namely ``bwat``) according to the gradient of
15 the modeled hydraulic potential. In both hydrology models a state variable ``tillwat`` is
16 the effective thickness of the layer of liquid water in the till; it is used to compute
17 the effective pressure on the till (see :ref:`sec-basestrength`). The pressure of the
18 transportable water ``bwat`` in the ``routing`` model does not relate directly to the
19 effective pressure on the till.
20
21 .. note::
22
23 Both ``null`` and ``routing`` described here provide all diagnostic quantities needed
24 for mass accounting, even though ``null`` is not mass-conserving. See
25 :ref:`sec-mass-conservation-hydrology` for details.
26
27 .. list-table:: Command-line options to choose the hydrology model
28 :name: tab-hydrologychoice
29 :header-rows: 1
30 :widths: 2,5
31
32 * - Option
33 - Description
34 * - :opt:`-hydrology null`
35 - The default model with only a layer of water stored in till. Not mass conserving in
36 the map-plane but much faster than ``-hydrology routing``. Based on "undrained
37 plastic bed" model of :cite:`Tulaczyketal2000b`. The only state variable is
38 ``tillwat``.
39
40 * - :opt:`-hydrology steady`
41 - A version of the ``null`` model that computes an approximation of the steady state
42 subglacial water flux that can be used as an input for a frontal melt
43 parameterization such as the one described in :ref:`sec-ismip6-frontal-melt`.
44
45 * - :opt:`-hydrology routing`
46 - A mass-conserving horizontal transport model in which the pressure of transportable
47 water is equal to overburden pressure. The till layer remains in the model, so this
48 is a "drained and conserved plastic bed" model. The state variables are ``bwat``
49 and ``tillwat``.
50
51 See :numref:`tab-hydrology` for options which apply to all hydrology models. Note
52 that the primary water source for these models is the energy conservation model which
53 computes the basal melt rate ``basal_melt_rate_grounded``. There is, however, also option
54 :opt:`-hydrology_input_to_bed_file` which allows the user to *add* water directly into the
55 subglacial layer, in addition to the computed ``basal_melt_rate_grounded`` values. Thus
56 ``-hydrology_input_to_bed_file`` allows the user to model drainage directly to the bed
57 from surface runoff, for example. Also option :opt:`-hydrology_bmelt_file` allows the user
58 to replace the computed ``basal_melt_rate_grounded`` rate by values read from a file,
59 thereby effectively decoupling the hydrology model from the ice dynamics
60 (esp. conservation of energy).
61
62 To use water runoff computed by a PISM surface model, set
63 :config:`hydrology.surface_input_from_runoff`. (The :ref:`sec-surface-pdd` computes runoff
64 using computed melt and the assumption that a fraction of this melt re-freezes, all other
65 models assume that all melt turns into runoff.)
66
67 .. list-table:: Subglacial hydrology command-line options which apply to all hydrology models
68 :name: tab-hydrology
69 :header-rows: 1
70 :widths: 1,1
71
72 * - Option
73 - Description
74 * - :opt:`-hydrology.surface_input.file`
75 - Specifies a NetCDF file which contains a time-dependent field ``water_input_rate``
76 which has units of water thickness per time. This rate is *added to* the basal melt
77 rate computed by the energy conservation code.
78 * - :opt:`-hydrology_tillwat_max` (m)
79 - Maximum effective thickness for water stored in till.
80 * - :opt:`-hydrology_tillwat_decay_rate` (m/a)
81 - Water accumulates in the till at the basal melt rate ``basal_melt_rate_grounded``,
82 minus this rate.
83 * - :opt:`-hydrology_use_const_bmelt`
84 - Replace the conservation-of-energy basal melt rate ``basal_melt_rate_grounded``
85 with a constant.
86
87 .. _sec-hydrology-null:
88
89 The default model: ``-hydrology null``
90 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
91
92 In this model the water is *not* conserved but it is stored locally in the till up to a
93 specified amount; the configuration parameter :config:`hydrology.tillwat_max` sets that
94 amount. The water is not conserved in the sense that water above the
95 :config:`hydrology.tillwat_max` level is lost permanently. This model is based on the
96 "undrained plastic bed" concept of :cite:`Tulaczyketal2000b`; see also
97 :cite:`BBssasliding`.
98
99 In particular, denoting ``tillwat`` by `W_{till}`, the till-stored water layer effective
100 thickness evolves by the simple equation
101
102 .. math::
103 :label: eq-tillwatevolve
104
105 \frac{\partial W_{till}}{\partial t} = \frac{m}{\rho_w} - C
106
107 where `m=` :var:`basal_melt_rate_grounded` (kg `\text{m}^{-2}\,\text{s}^{-1}`), `\rho_w`
108 is the density of fresh water, and `C` :var:`hydrology_tillwat_decay_rate`. At all times
109 bounds `0 \le W_{till} \le W_{till}^{max}` are satisfied.
110
111 This ``-hydrology null`` model has been extensively tested in combination with the
112 Mohr-Coulomb till (section :ref:`sec-basestrength`) for modelling ice streaming (see
113 :cite:`AschwandenAdalgeirsdottirKhroulev` and :cite:`BBssasliding`, among others).
114
115 .. _sec-hydrology-steady:
116
117 "Steady state flow" model
118 ~~~~~~~~~~~~~~~~~~~~~~~~~
119
120 This model (``-hydrology steady``) adds an approximation of the subglacial water flux to
121 the ``null`` model. It *does not* model the steady state distribution of the subglacial
122 water thickness.
123
124 Here we assume that the water input from the surface read from the file specified using
125 :config:`hydrology.surface_input.file` instantaneously percolates to the base of the ice
126 and enters the subglacial water system.
127
128 We also assume that the subglacial drainage system instantaneously reaches its steady
129 state. Under this assumption the water flux can be approximated by using an auxiliary
130 problem that is computationally cheaper to solve, compared to using the mass conserving
131 ``routing`` model (at least at high grid resolutions).
132
133 We solve
134
135 .. math::
136 :label: eq-steady-hydro-aux
137
138 \diff{u}{t} = -\div (\V u)
139
140 on the grounded part of the domain with the initial state `u_0 = \tau M`, where `\tau` is
141 the scaling of the input rate (:config:`hydrology.steady.input_rate_scaling`) and `M` is
142 the input rate read from an input file.
143
144 The velocity field `\V` approximates the steady state flow direction. In this
145 implementation `\V` is the normalized gradient of
146
147 .. math::
148
149 \psi = \rho_i g H + \rho_w g b + \Delta \psi.
150
151 Here `H` is the ice thickness, `b` is the bed elevation, `g` is the acceleration due to
152 gravity and `\rho_i`, `\rho_w` are densities of ice and water, respectively.
153
154 .. note::
155
156 Note that `\psi` ignores subglacial water thickness, producing flow patterns that are
157 more concentrated than would be seen in a model that includes it. This effect is
158 grid-dependent.
159
160 This implies that this code cannot model the distribution of the flow along the
161 grounding line of a particular outlet glacier but it *can* tell us how surface runoff
162 is distributed among the outlets.
163
164 The term `\Delta \psi` is the adjustment needed to remove internal minima from the "raw"
165 potential, filling any "lakes" it might have. This modification of `\psi` is performed
166 using an iterative process; set :config:`hydrology.steady.potential_n_iterations` to
167 control the maximum number of iterations and :config:`hydrology.steady.potential_delta` to
168 affect the amount it is adjusted by at each iteration.
169
170 The equation :eq:`eq-steady-hydro-aux` is advanced forward in time until `\int_{\Omega}u
171 < \epsilon\int_{\Omega} u_0`, where `\epsilon` (:config:`hydrology.steady.volume_ratio`)
172 is a fraction controlling the accuracy of the approximation: lower values of `\epsilon`
173 would result in a better approximation and a higher computational cost (more iterations).
174 Set :config:`hydrology.steady.n_iterations` to control the maximum number of these
175 iterations.
176
177 This model restricts the time step length in order to capture the temporal variability of
178 the forcing: the flux is updated at least once for each time interval in the forcing file.
179
180 In simulations with a slowly-varying surface input rate the flux has to be updated every
181 once in a while to reflect changes in the flow pattern coming from changing geometry. Use
182 :config:`hydrology.steady.flux_update_interval` (years) to set the update frequency.
183
184 See :ref:`sec-steady-hydro` for technical details.
185
186 .. _sec-hydrology-routing:
187
188 The mass-conserving model: ``-hydrology routing``
189 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
190
191 In this model the water *is* conserved in the map-plane. Water does get put into the till,
192 with the same maximum value :config:`hydrology.tillwat_max`, but excess water is
193 horizontally-transported. An additional state variable ``bwat``, the effective thickness
194 of the layer of transportable water, is used by ``routing``. This transportable water will
195 flow in the direction of the negative of the gradient of the modeled hydraulic potential.
196 In the ``routing`` model this potential is calculated by assuming that the transportable
197 subglacial water is at the overburden pressure :cite:`Siegertetal2007`. Ultimately the
198 transportable water will reach the ice sheet grounding line or ice-free-land margin, at
199 which point it will be lost. The amount that is lost this way is reported to the user.
200
201 In this model ``tillwat`` also evolves by equation :eq:`eq-tillwatevolve`, but several
202 additional parameters are used in determining how the transportable water ``bwat`` flows
203 in the model; see :numref:`tab-hydrologyrouting`. Specifically, the horizontal
204 subglacial water flux is determined by a generalized Darcy flux relation :cite:`Clarke05`,
205 :cite:`Schoofetal2012`
206
207 .. math::
208 :label: eq-flux
209
210 \bq = - k\, W^\alpha\, |\nabla \psi|^{\beta-2} \nabla \psi
211
212 where `\bq` is the lateral water flux, `W=` ``bwat`` is the effective thickness of the
213 layer of transportable water, `\psi` is the hydraulic potential, and `k`, `\alpha`,
214 `\beta` are controllable parameters (:numref:`tab-hydrologyrouting`).
215
216 In the ``routing`` model the hydraulic potential `\psi` is determined by
217
218 .. math::
219 :label: eq-hydraulicpotential
220
221 \psi = P_o + \rho_w g (b + W)
222
223 where `P_o=\rho_i g H` is the ice overburden pressure, `g` is gravity, `\rho_i` is ice
224 density, `\rho_w` is fresh water density, `H` is ice thickness, and `b` is the bedrock
225 elevation.
226
227 For most choices of the relevant parameters and most grid spacings, the ``routing`` model
228 is at least two orders of magnitude more expensive computationally than the ``null``
229 model. This follows directly from the CFL-type time-step restriction on lateral flow of a
230 fluid with velocity on the order of centimeters to meters per second (i.e. the subglacial
231 liquid water ``bwat``). (By comparison, much of PISM ice dynamics time-stepping is
232 controlled by the much slower velocity of the flowing ice.) Therefore the user should
233 start with short runs of order a few model years. We also recommend ``daily`` or even
234 ``hourly`` reporting for scalar and spatially-distributed time-series to see hydrology
235 model behavior, especially on fine grids (e.g. `< 1` km).
236
237 .. list-table:: Command-line options specific to hydrology model ``routing``
238 :name: tab-hydrologyrouting
239 :header-rows: 1
240 :widths: 5,4
241
242 * - Option
243 - Description
244 * - :opt:`-hydrology_hydraulic_conductivity` `k`
245 - `=k` in formula :eq:`eq-flux`.
246 * - :opt:`-hydrology_null_strip` (km)
247 - In the boundary strip water is removed and this is reported. This option specifies
248 the width of this strip, which should typically be one or two grid cells.
249 * - :opt:`-hydrology_gradient_power_in_flux` `\beta`
250 - `=\beta` in formula :eq:`eq-flux`.
251 * - :opt:`-hydrology_thickness_power_in_flux` `\alpha`
252 - `=\alpha` in formula :eq:`eq-flux`.
253
254 .. FIXME -hydrology distributed is not documented except by :cite:`BuelervanPelt2015`