tRouting.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tRouting.cc (33984B)
---
1 // Copyright (C) 2012-2019 PISM Authors
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include <cassert>
20
21 #include "Routing.hh"
22 #include "pism/util/IceModelVec2CellType.hh"
23 #include "pism/util/Mask.hh"
24 #include "pism/util/MaxTimestep.hh"
25
26 #include "pism/util/error_handling.hh"
27
28 #include "pism/util/pism_options.hh"
29 #include "pism/util/pism_utilities.hh"
30 #include "pism/util/Vars.hh"
31 #include "pism/geometry/Geometry.hh"
32 #include "pism/util/Profiling.hh"
33
34 namespace pism {
35 namespace hydrology {
36
37 namespace diagnostics {
38
39 //! \brief Reports the pressure of the transportable water in the subglacial layer.
40 class BasalWaterPressure : public Diag<Routing>
41 {
42 public:
43 BasalWaterPressure(const Routing *m)
44 : Diag<Routing>(m) {
45 m_vars = {SpatialVariableMetadata(m_sys, "bwp")};
46 set_attrs("pressure of transportable water in subglacial layer", "", "Pa", "Pa", 0);
47 }
48
49 protected:
50 virtual IceModelVec::Ptr compute_impl() const {
51 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "bwp", WITHOUT_GHOSTS));
52 result->metadata() = m_vars[0];
53 result->copy_from(model->subglacial_water_pressure());
54 return result;
55 }
56 };
57
58
59 //! \brief Reports the pressure of the transportable water in the subglacial layer as a
60 //! fraction of the overburden pressure.
61 class RelativeBasalWaterPressure : public Diag<Routing>
62 {
63 public:
64 RelativeBasalWaterPressure(const Routing *m)
65 : Diag<Routing>(m) {
66 m_vars = {SpatialVariableMetadata(m_sys, "bwprel")};
67 set_attrs("pressure of transportable water in subglacial layer"
68 " as fraction of the overburden pressure", "",
69 "", "", 0);
70 m_vars[0].set_number("_FillValue", m_fill_value);
71 }
72
73 protected:
74 virtual IceModelVec::Ptr compute_impl() const {
75 double fill_value = m_fill_value;
76
77 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "bwprel", WITHOUT_GHOSTS));
78 result->metadata(0) = m_vars[0];
79
80 const IceModelVec2S
81 &P = model->subglacial_water_pressure(),
82 &Po = model->overburden_pressure();
83
84 IceModelVec::AccessList list{result.get(), &Po, &P};
85 for (Points p(*m_grid); p; p.next()) {
86 const int i = p.i(), j = p.j();
87
88 if (Po(i,j) > 0.0) {
89 (*result)(i,j) = P(i, j) / Po(i,j);
90 } else {
91 (*result)(i,j) = fill_value;
92 }
93 }
94
95 return result;
96 }
97 };
98
99
100 //! \brief Reports the effective pressure of the transportable water in the subglacial
101 //! layer, that is, the overburden pressure minus the pressure.
102 class EffectiveBasalWaterPressure : public Diag<Routing>
103 {
104 public:
105 EffectiveBasalWaterPressure(const Routing *m)
106 : Diag<Routing>(m) {
107 m_vars = {SpatialVariableMetadata(m_sys, "effbwp")};
108 set_attrs("effective pressure of transportable water in subglacial layer"
109 " (overburden pressure minus water pressure)",
110 "", "Pa", "Pa", 0);
111 }
112
113 protected:
114 virtual IceModelVec::Ptr compute_impl() const {
115
116 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "effbwp", WITHOUT_GHOSTS));
117 result->metadata() = m_vars[0];
118
119 const IceModelVec2S
120 &P = model->subglacial_water_pressure(),
121 &Po = model->overburden_pressure();
122
123 IceModelVec::AccessList list{&Po, &P, result.get()};
124
125 for (Points p(*m_grid); p; p.next()) {
126 const int i = p.i(), j = p.j();
127
128 (*result)(i, j) = Po(i, j) - P(i, j);
129 }
130
131 return result;
132 }
133 };
134
135
136 //! \brief Report the wall melt rate from dissipation of the potential energy of the
137 //! transportable water.
138 class WallMelt : public Diag<Routing>
139 {
140 public:
141 WallMelt(const Routing *m)
142 : Diag<Routing>(m) {
143 m_vars = {SpatialVariableMetadata(m_sys, "wallmelt")};
144 set_attrs("wall melt into subglacial hydrology layer"
145 " from (turbulent) dissipation of energy in transportable water",
146 "", "m s-1", "m year-1", 0);
147 }
148
149 protected:
150 virtual IceModelVec::Ptr compute_impl() const {
151 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "wallmelt", WITHOUT_GHOSTS));
152 result->metadata() = m_vars[0];
153
154 const IceModelVec2S &bed_elevation = *m_grid->variables().get_2d_scalar("bedrock_altitude");
155
156 wall_melt(*model, bed_elevation, *result);
157 return result;
158 }
159 };
160
161 //! @brief Diagnostically reports the staggered-grid components of the velocity of the
162 //! water in the subglacial layer.
163 class BasalWaterVelocity : public Diag<Routing>
164 {
165 public:
166 BasalWaterVelocity(const Routing *m)
167 : Diag<Routing>(m) {
168
169 // set metadata:
170 m_vars = {SpatialVariableMetadata(m_sys, "bwatvel[0]"),
171 SpatialVariableMetadata(m_sys, "bwatvel[1]")};
172
173 set_attrs("velocity of water in subglacial layer, i-offset", "",
174 "m s-1", "m s-1", 0);
175 set_attrs("velocity of water in subglacial layer, j-offset", "",
176 "m s-1", "m s-1", 1);
177 }
178 protected:
179 virtual IceModelVec::Ptr compute_impl() const {
180 IceModelVec2Stag::Ptr result(new IceModelVec2Stag(m_grid, "bwatvel", WITHOUT_GHOSTS));
181 result->metadata(0) = m_vars[0];
182 result->metadata(1) = m_vars[1];
183
184 result->copy_from(model->velocity_staggered());
185
186 return result;
187 }
188 };
189
190 //! Compute the hydraulic potential.
191 /*!
192 Computes \f$\psi = P + \rho_w g (b + W)\f$.
193 */
194 void hydraulic_potential(const IceModelVec2S &W,
195 const IceModelVec2S &P,
196 const IceModelVec2S &sea_level,
197 const IceModelVec2S &bed,
198 const IceModelVec2S &ice_thickness,
199 IceModelVec2S &result) {
200
201 IceGrid::ConstPtr grid = result.grid();
202
203 Config::ConstPtr config = grid->ctx()->config();
204
205 double
206 ice_density = config->get_number("constants.ice.density"),
207 sea_water_density = config->get_number("constants.sea_water.density"),
208 C = ice_density / sea_water_density,
209 rg = (config->get_number("constants.fresh_water.density") *
210 config->get_number("constants.standard_gravity"));
211
212 IceModelVec::AccessList list{&P, &W, &sea_level, &ice_thickness, &bed, &result};
213
214 for (Points p(*grid); p; p.next()) {
215 const int i = p.i(), j = p.j();
216
217 double b = std::max(bed(i, j), sea_level(i, j) - C * ice_thickness(i, j));
218
219 result(i, j) = P(i, j) + rg * (b + W(i, j));
220 }
221 }
222
223 /*! @brief Report hydraulic potential in the subglacial hydrology system */
224 class HydraulicPotential : public Diag<Routing>
225 {
226 public:
227 HydraulicPotential(const Routing *m)
228 : Diag<Routing>(m) {
229
230 m_vars = {SpatialVariableMetadata(m_sys, "hydraulic_potential")};
231
232 set_attrs("hydraulic potential in the subglacial hydrology system", "",
233 "Pa", "Pa", 0);
234 }
235
236 protected:
237 IceModelVec::Ptr compute_impl() const {
238
239 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "hydraulic_potential", WITHOUT_GHOSTS));
240 result->metadata(0) = m_vars[0];
241
242 const IceModelVec2S &sea_level = *m_grid->variables().get_2d_scalar("sea_level");
243 const IceModelVec2S &bed_elevation = *m_grid->variables().get_2d_scalar("bedrock_altitude");
244 const IceModelVec2S &ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
245
246 hydraulic_potential(model->subglacial_water_thickness(),
247 model->subglacial_water_pressure(),
248 sea_level,
249 bed_elevation,
250 ice_thickness,
251 *result);
252
253 return result;
254 }
255 };
256
257 } // end of namespace diagnostics
258
259 Routing::Routing(IceGrid::ConstPtr grid)
260 : Hydrology(grid),
261 m_Qstag(grid, "advection_flux", WITH_GHOSTS, 1),
262 m_Qstag_average(grid, "cumulative_advection_flux", WITH_GHOSTS, 1),
263 m_Vstag(grid, "water_velocity", WITHOUT_GHOSTS),
264 m_Wstag(grid, "W_staggered", WITH_GHOSTS, 1),
265 m_Kstag(grid, "K_staggered", WITH_GHOSTS, 1),
266 m_Wnew(grid, "W_new", WITHOUT_GHOSTS),
267 m_Wtillnew(grid, "Wtill_new", WITHOUT_GHOSTS),
268 m_R(grid, "potential_workspace", WITH_GHOSTS, 1), /* box stencil used */
269 m_dx(grid->dx()),
270 m_dy(grid->dy()),
271 m_bottom_surface(grid, "ice_bottom_surface_elevation", WITH_GHOSTS) {
272
273 m_W.metadata().set_string("pism_intent", "model_state");
274
275 m_rg = (m_config->get_number("constants.fresh_water.density") *
276 m_config->get_number("constants.standard_gravity"));
277
278 m_Qstag.set_attrs("internal",
279 "cell face-centered (staggered) components of advective subglacial water flux",
280 "m2 s-1", "m2 s-1", "", 0);
281
282 m_Qstag_average.set_attrs("internal",
283 "average (over time) advection flux on the staggered grid",
284 "m2 s-1", "m2 s-1", "", 0);
285
286 m_Vstag.set_attrs("internal",
287 "cell face-centered (staggered) components of water velocity"
288 " in subglacial water layer",
289 "m s-1", "m s-1", "", 0);
290
291 // auxiliary variables which NEED ghosts
292 m_Wstag.set_attrs("internal",
293 "cell face-centered (staggered) values of water layer thickness",
294 "m", "m", "", 0);
295 m_Wstag.metadata().set_number("valid_min", 0.0);
296
297 m_Kstag.set_attrs("internal",
298 "cell face-centered (staggered) values of nonlinear conductivity",
299 "", "", "", 0);
300 m_Kstag.metadata().set_number("valid_min", 0.0);
301
302 m_R.set_attrs("internal",
303 "work space for modeled subglacial water hydraulic potential",
304 "Pa", "Pa", "", 0);
305
306 // temporaries during update; do not need ghosts
307 m_Wnew.set_attrs("internal",
308 "new thickness of transportable subglacial water layer during update",
309 "m", "m", "", 0);
310 m_Wnew.metadata().set_number("valid_min", 0.0);
311
312 m_Wtillnew.set_attrs("internal",
313 "new thickness of till (subglacial) water layer during update",
314 "m", "m", "", 0);
315 m_Wtillnew.metadata().set_number("valid_min", 0.0);
316
317 {
318 double alpha = m_config->get_number("hydrology.thickness_power_in_flux");
319 if (alpha < 1.0) {
320 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
321 "alpha = %f < 1 which is not allowed", alpha);
322 }
323
324 if (m_config->get_number("hydrology.tillwat_max") < 0.0) {
325 throw RuntimeError(PISM_ERROR_LOCATION,
326 "hydrology::Routing: hydrology.tillwat_max is negative.\n"
327 "This is not allowed.");
328 }
329 }
330 }
331
332 Routing::~Routing() {
333 // empty
334 }
335
336 void Routing::initialization_message() const {
337 m_log->message(2,
338 "* Initializing the routing subglacial hydrology model ...\n");
339
340 if (m_config->get_flag("hydrology.routing.include_floating_ice")) {
341 m_log->message(2, " ... routing subglacial water under grounded and floating ice.\n");
342 } else {
343 m_log->message(2, " ... routing subglacial water under grounded ice only.\n");
344 }
345 }
346
347 void Routing::restart_impl(const File &input_file, int record) {
348 Hydrology::restart_impl(input_file, record);
349
350 m_W.read(input_file, record);
351
352 regrid("Hydrology", m_W);
353 }
354
355 void Routing::bootstrap_impl(const File &input_file,
356 const IceModelVec2S &ice_thickness) {
357 Hydrology::bootstrap_impl(input_file, ice_thickness);
358
359 double bwat_default = m_config->get_number("bootstrapping.defaults.bwat");
360 m_W.regrid(input_file, OPTIONAL, bwat_default);
361
362 regrid("Hydrology", m_W);
363 }
364
365 void Routing::init_impl(const IceModelVec2S &W_till,
366 const IceModelVec2S &W,
367 const IceModelVec2S &P) {
368 Hydrology::init_impl(W_till, W, P);
369
370 m_W.copy_from(W);
371 }
372
373 void Routing::define_model_state_impl(const File &output) const {
374 Hydrology::define_model_state_impl(output);
375 m_W.define(output);
376 }
377
378 void Routing::write_model_state_impl(const File &output) const {
379 Hydrology::write_model_state_impl(output);
380 m_W.write(output);
381 }
382
383 //! Returns the (trivial) overburden pressure as the pressure of the transportable water,
384 //! because this is the model.
385 const IceModelVec2S& Routing::subglacial_water_pressure() const {
386 return m_Pover;
387 }
388
389 const IceModelVec2Stag& Routing::velocity_staggered() const {
390 return m_Vstag;
391 }
392
393
394 //! Average the regular grid water thickness to values at the center of cell edges.
395 /*! Uses mask values to avoid averaging using water thickness values from
396 either ice-free or floating areas. */
397 void Routing::water_thickness_staggered(const IceModelVec2S &W,
398 const IceModelVec2CellType &mask,
399 IceModelVec2Stag &result) {
400
401 bool include_floating = m_config->get_flag("hydrology.routing.include_floating_ice");
402
403 IceModelVec::AccessList list{ &mask, &W, &result };
404
405 for (Points p(*m_grid); p; p.next()) {
406 const int i = p.i(), j = p.j();
407
408 if (include_floating) {
409 // east
410 if (mask.icy(i, j)) {
411 result(i, j, 0) = mask.icy(i + 1, j) ? 0.5 * (W(i, j) + W(i + 1, j)) : W(i, j);
412 } else {
413 result(i, j, 0) = mask.icy(i + 1, j) ? W(i + 1, j) : 0.0;
414 }
415 // north
416 if (mask.icy(i, j)) {
417 result(i, j, 1) = mask.icy(i, j + 1) ? 0.5 * (W(i, j) + W(i, j + 1)) : W(i, j);
418 } else {
419 result(i, j, 1) = mask.icy(i, j + 1) ? W(i, j + 1) : 0.0;
420 }
421 } else {
422 // east
423 if (mask.grounded_ice(i, j)) {
424 result(i, j, 0) = mask.grounded_ice(i + 1, j) ? 0.5 * (W(i, j) + W(i + 1, j)) : W(i, j);
425 } else {
426 result(i, j, 0) = mask.grounded_ice(i + 1, j) ? W(i + 1, j) : 0.0;
427 }
428 // north
429 if (mask.grounded_ice(i, j)) {
430 result(i, j, 1) = mask.grounded_ice(i, j + 1) ? 0.5 * (W(i, j) + W(i, j + 1)) : W(i, j);
431 } else {
432 result(i, j, 1) = mask.grounded_ice(i, j + 1) ? W(i, j + 1) : 0.0;
433 }
434 }
435 }
436
437 result.update_ghosts();
438 }
439
440
441 //! Compute the nonlinear conductivity at the center of cell edges.
442 /*!
443 Computes
444
445 \f[ K = K(W, \nabla P, \nabla b) = k W^{\alpha-1} |\nabla R|^{\beta-2} \f]
446
447 on the staggered grid, where \f$R = P+\rho_w g b\f$. We denote
448
449 \f[ \Pi = |\nabla R|^2 \f]
450
451 internally; this is computed on a staggered grid by a Mahaffy-like ([@ref Mahaffy])
452 scheme. This requires \f$R\f$ to be defined on a box stencil of width 1.
453
454 Also returns the maximum over all staggered points of \f$ K W \f$.
455 */
456 void Routing::compute_conductivity(const IceModelVec2Stag &W,
457 const IceModelVec2S &P,
458 const IceModelVec2S &bed_elevation,
459 IceModelVec2Stag &result,
460 double &KW_max) const {
461 const double
462 k = m_config->get_number("hydrology.hydraulic_conductivity"),
463 alpha = m_config->get_number("hydrology.thickness_power_in_flux"),
464 beta = m_config->get_number("hydrology.gradient_power_in_flux"),
465 betapow = (beta - 2.0) / 2.0;
466
467 IceModelVec::AccessList list({&result, &W});
468
469 KW_max = 0.0;
470
471 if (beta != 2.0) {
472 // Put the squared norm of the gradient of the simplified hydrolic potential (Pi) in
473 // "result"
474 //
475 // FIXME: we don't need to re-compute this during every hydrology time step: the
476 // simplified hydrolic potential does not depend on the water amount and can be
477 // computed *once* in update_impl(), before entering the time-stepping loop
478 {
479 // R <-- P + rhow g b
480 P.add(m_rg, bed_elevation, m_R); // yes, it updates ghosts
481
482 list.add(m_R);
483 for (Points p(*m_grid); p; p.next()) {
484 const int i = p.i(), j = p.j();
485
486 double dRdx, dRdy;
487 dRdx = (m_R(i + 1, j) - m_R(i, j)) / m_dx;
488 dRdy = (m_R(i + 1, j + 1) + m_R(i, j + 1) - m_R(i + 1, j - 1) - m_R(i, j - 1)) / (4.0 * m_dy);
489 result(i, j, 0) = dRdx * dRdx + dRdy * dRdy;
490
491 dRdx = (m_R(i + 1, j + 1) + m_R(i + 1, j) - m_R(i - 1, j + 1) - m_R(i - 1, j)) / (4.0 * m_dx);
492 dRdy = (m_R(i, j + 1) - m_R(i, j)) / m_dy;
493 result(i, j, 1) = dRdx * dRdx + dRdy * dRdy;
494 }
495 }
496
497 // We regularize negative power |\grad psi|^{beta-2} by adding eps because large
498 // head gradient might be 10^7 Pa per 10^4 m or 10^3 Pa/m.
499 const double eps = beta < 2.0 ? 1.0 : 0.0;
500
501 for (Points p(*m_grid); p; p.next()) {
502 const int i = p.i(), j = p.j();
503
504 for (int o = 0; o < 2; ++o) {
505 const double Pi = result(i, j, o);
506
507 // FIXME: same as Pi above: we don't need to re-compute this each time we make a
508 // short hydrology time step
509 double B = pow(Pi + eps * eps, betapow);
510
511 result(i, j, o) = k * pow(W(i, j, o), alpha - 1.0) * B;
512
513 KW_max = std::max(KW_max, result(i, j, o) * W(i, j, o));
514 }
515 }
516 } else {
517 for (Points p(*m_grid); p; p.next()) {
518 const int i = p.i(), j = p.j();
519
520 for (int o = 0; o < 2; ++o) {
521 result(i, j, o) = k * pow(W(i, j, o), alpha - 1.0);
522
523 KW_max = std::max(KW_max, result(i, j, o) * W(i, j, o));
524 }
525 }
526 }
527
528 KW_max = GlobalMax(m_grid->com, KW_max);
529
530 result.update_ghosts();
531 }
532
533
534 //! Compute the wall melt rate which comes from (turbulent) dissipation of flow energy.
535 /*!
536 This code fills `result` with
537 \f[ \frac{m_{wall}}{\rho_w} = - \frac{1}{L \rho_w} \mathbf{q} \cdot \nabla \psi = \left(\frac{k}{L \rho_w}\right) W^\alpha |\nabla R|^\beta \f]
538 where \f$R = P+\rho_w g b\f$.
539
540 Note that conductivity_staggered() computes the related quantity
541 \f$K = k W^{\alpha-1} |\nabla R|^{\beta-2}\f$ on the staggered grid, but
542 contriving to reuse that code would be inefficient because of the
543 staggered-versus-regular change.
544
545 At the current state of the code, this is a diagnostic calculation only.
546 */
547 void wall_melt(const Routing &model,
548 const IceModelVec2S &bed_elevation,
549 IceModelVec2S &result) {
550
551 IceGrid::ConstPtr grid = result.grid();
552
553 Config::ConstPtr config = grid->ctx()->config();
554
555 const double
556 k = config->get_number("hydrology.hydraulic_conductivity"),
557 L = config->get_number("constants.fresh_water.latent_heat_of_fusion"),
558 alpha = config->get_number("hydrology.thickness_power_in_flux"),
559 beta = config->get_number("hydrology.gradient_power_in_flux"),
560 g = config->get_number("constants.standard_gravity"),
561 rhow = config->get_number("constants.fresh_water.density"),
562 rg = rhow * g,
563 CC = k / (L * rhow);
564
565 // FIXME: could be scaled with overall factor hydrology_coefficient_wall_melt ?
566 if (alpha < 1.0) {
567 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
568 "alpha = %f < 1 which is not allowed", alpha);
569 }
570
571 IceModelVec2S R(grid, "R", WITH_GHOSTS);
572
573 // R <-- P + rhow g b
574 model.subglacial_water_pressure().add(rg, bed_elevation, R);
575 // yes, it updates ghosts
576
577 IceModelVec2S W(grid, "W", WITH_GHOSTS);
578 W.copy_from(model.subglacial_water_thickness());
579
580 IceModelVec::AccessList list{&R, &W, &result};
581
582 double dx = grid->dx();
583 double dy = grid->dy();
584
585 for (Points p(*grid); p; p.next()) {
586 const int i = p.i(), j = p.j();
587 double dRdx, dRdy;
588
589 if (W(i, j) > 0.0) {
590 dRdx = 0.0;
591 if (W(i + 1, j) > 0.0) {
592 dRdx = (R(i + 1, j) - R(i, j)) / (2.0 * dx);
593 }
594 if (W(i - 1, j) > 0.0) {
595 dRdx += (R(i, j) - R(i - 1, j)) / (2.0 * dx);
596 }
597 dRdy = 0.0;
598 if (W(i, j + 1) > 0.0) {
599 dRdy = (R(i, j + 1) - R(i, j)) / (2.0 * dy);
600 }
601 if (W(i, j - 1) > 0.0) {
602 dRdy += (R(i, j) - R(i, j - 1)) / (2.0 * dy);
603 }
604 result(i, j) = CC * pow(W(i, j), alpha) * pow(dRdx * dRdx + dRdy * dRdy, beta/2.0);
605 } else {
606 result(i, j) = 0.0;
607 }
608 }
609 }
610
611
612 //! Get the advection velocity V at the center of cell edges.
613 /*!
614 Computes the advection velocity @f$\mathbf{V}@f$ on the staggered
615 (edge-centered) grid. If V = (u, v) in components then we have
616 <code> result(i, j, 0) = u(i+1/2, j) </code> and
617 <code> result(i, j, 1) = v(i, j+1/2) </code>
618
619 The advection velocity is given by the formula
620
621 @f[ \mathbf{V} = - K \left(\nabla P + \rho_w g \nabla b\right) @f]
622
623 where @f$\mathbf{V}@f$ is the water velocity, @f$P@f$ is the water
624 pressure, and @f$b@f$ is the bedrock elevation.
625
626 If the corresponding staggered grid value of the water thickness is zero then that
627 component of V is set to zero. This does not change the flux value (which would be zero
628 anyway) but it does provide the correct max velocity in the CFL calculation. We assume
629 bed has valid ghosts.
630 */
631 void Routing::compute_velocity(const IceModelVec2Stag &W,
632 const IceModelVec2S &pressure,
633 const IceModelVec2S &bed,
634 const IceModelVec2Stag &K,
635 const IceModelVec2Int *no_model_mask,
636 IceModelVec2Stag &result) const {
637 IceModelVec2S &P = m_R;
638 P.copy_from(pressure); // yes, it updates ghosts
639
640 IceModelVec::AccessList list{&P, &W, &K, &bed, &result};
641
642 for (Points p(*m_grid); p; p.next()) {
643 const int i = p.i(), j = p.j();
644
645 if (W(i, j, 0) > 0.0) {
646 double
647 P_x = (P(i + 1, j) - P(i, j)) / m_dx,
648 b_x = (bed(i + 1, j) - bed(i, j)) / m_dx;
649 result(i, j, 0) = - K(i, j, 0) * (P_x + m_rg * b_x);
650 } else {
651 result(i, j, 0) = 0.0;
652 }
653
654 if (W(i, j, 1) > 0.0) {
655 double
656 P_y = (P(i, j + 1) - P(i, j)) / m_dy,
657 b_y = (bed(i, j + 1) - bed(i, j)) / m_dy;
658 result(i, j, 1) = - K(i, j, 1) * (P_y + m_rg * b_y);
659 } else {
660 result(i, j, 1) = 0.0;
661 }
662 }
663
664 if (no_model_mask) {
665 list.add(*no_model_mask);
666
667 for (Points p(*m_grid); p; p.next()) {
668 const int i = p.i(), j = p.j();
669
670 auto M = no_model_mask->int_star(i, j);
671
672 if (M.ij or M.e) {
673 result(i, j, 0) = 0.0;
674 }
675
676 if (M.ij or M.n) {
677 result(i, j, 1) = 0.0;
678 }
679 }
680 }
681 }
682
683
684 //! Compute Q = V W at edge-centers (staggered grid) by first-order upwinding.
685 /*!
686 The field W must have valid ghost values, but V does not need them.
687
688 FIXME: This could be re-implemented using the Koren (1993) flux-limiter.
689 */
690 void Routing::advective_fluxes(const IceModelVec2Stag &V,
691 const IceModelVec2S &W,
692 IceModelVec2Stag &result) const {
693 IceModelVec::AccessList list{&W, &V, &result};
694
695 assert(W.stencil_width() >= 1);
696
697 for (Points p(*m_grid); p; p.next()) {
698 const int i = p.i(), j = p.j();
699
700 result(i, j, 0) = V(i, j, 0) * (V(i, j, 0) >= 0.0 ? W(i, j) : W(i + 1, j));
701 result(i, j, 1) = V(i, j, 1) * (V(i, j, 1) >= 0.0 ? W(i, j) : W(i, j + 1));
702 }
703
704 result.update_ghosts();
705 }
706
707 /*!
708 * See equation (51) in Bueler and van Pelt.
709 */
710 double Routing::max_timestep_W_diff(double KW_max) const {
711 double D_max = m_rg * KW_max;
712 double result = 1.0 / (m_dx * m_dx) + 1.0 / (m_dy * m_dy);
713 return 0.25 / (D_max * result);
714 }
715
716 /*!
717 * See equation (50) in Bueler and van Pelt.
718 */
719 double Routing::max_timestep_W_cfl() const {
720 // V could be zero if P is constant and bed is flat
721 std::vector<double> tmp = m_Vstag.absmaxcomponents();
722
723 // add a safety margin
724 double alpha = 0.95;
725 double eps = 1e-6;
726
727 return alpha * 0.5 / (tmp[0]/m_dx + tmp[1]/m_dy + eps);
728 }
729
730
731 //! The computation of Wtillnew, called by update().
732 /*!
733 Does a step of the trivial integration
734 \f[ \frac{\partial W_{till}}{\partial t} = \frac{m}{\rho_w} - C\f]
735
736 where \f$C=\f$`hydrology_tillwat_decay_rate`. Enforces bounds
737 \f$0 \le W_{till} \le W_{till}^{max}\f$ where the upper bound is
738 `hydrology_tillwat_max`. Here \f$m/\rho_w\f$ is `total_input`.
739
740 Compare hydrology::NullTransport::update_impl().
741
742 The current code is not quite "code duplication" because the code here:
743
744 1. computes `Wtill_new` instead of updating `Wtill` in place;
745 2. uses time steps determined by the rest of the hydrology::Routing model;
746 3. does not check mask because the enforce_bounds() call addresses that.
747
748 Otherwise this is the same physical model with the same configurable parameters.
749 */
750 void Routing::update_Wtill(double dt,
751 const IceModelVec2S &Wtill,
752 const IceModelVec2S &surface_input_rate,
753 const IceModelVec2S &basal_melt_rate,
754 IceModelVec2S &Wtill_new) {
755 const double
756 tillwat_max = m_config->get_number("hydrology.tillwat_max"),
757 C = m_config->get_number("hydrology.tillwat_decay_rate", "m / second");
758
759 IceModelVec::AccessList list{&Wtill, &Wtill_new, &basal_melt_rate};
760
761 bool add_surface_input = m_config->get_flag("hydrology.add_water_input_to_till_storage");
762 if (add_surface_input) {
763 list.add(surface_input_rate);
764 }
765
766 for (Points p(*m_grid); p; p.next()) {
767 const int i = p.i(), j = p.j();
768
769 double input_rate = basal_melt_rate(i, j);
770 if (add_surface_input) {
771 input_rate += surface_input_rate(i, j);
772 }
773
774 Wtill_new(i, j) = clip(Wtill(i, j) + dt * (input_rate - C),
775 0, tillwat_max);
776 }
777 }
778
779 void Routing::W_change_due_to_flow(double dt,
780 const IceModelVec2S &W,
781 const IceModelVec2Stag &Wstag,
782 const IceModelVec2Stag &K,
783 const IceModelVec2Stag &Q,
784 IceModelVec2S &result) {
785 const double
786 wux = 1.0 / (m_dx * m_dx),
787 wuy = 1.0 / (m_dy * m_dy);
788
789 IceModelVec::AccessList list{&W, &Wstag, &K, &Q, &result};
790
791 for (Points p(*m_grid); p; p.next()) {
792 const int i = p.i(), j = p.j();
793
794 auto q = Q.star(i, j);
795 const double divQ = (q.e - q.w) / m_dx + (q.n - q.s) / m_dy;
796
797 auto k = K.star(i, j);
798 auto ws = Wstag.star(i, j);
799
800 const double
801 De = m_rg * k.e * ws.e,
802 Dw = m_rg * k.w * ws.w,
803 Dn = m_rg * k.n * ws.n,
804 Ds = m_rg * k.s * ws.s;
805
806 auto w = W.star(i, j);
807 const double diffW = (wux * (De * (w.e - w.ij) - Dw * (w.ij - w.w)) +
808 wuy * (Dn * (w.n - w.ij) - Ds * (w.ij - w.s)));
809
810 result(i, j) = dt * (- divQ + diffW);
811 }
812 }
813
814
815 //! The computation of Wnew, called by update().
816 void Routing::update_W(double dt,
817 const IceModelVec2S &surface_input_rate,
818 const IceModelVec2S &basal_melt_rate,
819 const IceModelVec2S &W,
820 const IceModelVec2Stag &Wstag,
821 const IceModelVec2S &Wtill,
822 const IceModelVec2S &Wtill_new,
823 const IceModelVec2Stag &K,
824 const IceModelVec2Stag &Q,
825 IceModelVec2S &W_new) {
826
827 W_change_due_to_flow(dt, W, Wstag, K, Q, m_flow_change_incremental);
828
829 IceModelVec::AccessList list{&W, &Wtill, &Wtill_new, &surface_input_rate,
830 &basal_melt_rate, &m_flow_change_incremental, &W_new};
831
832 for (Points p(*m_grid); p; p.next()) {
833 const int i = p.i(), j = p.j();
834
835 double input_rate = surface_input_rate(i, j) + basal_melt_rate(i, j);
836
837 double Wtill_change = Wtill_new(i, j) - Wtill(i, j);
838 W_new(i, j) = (W(i, j) + (dt * input_rate - Wtill_change) + m_flow_change_incremental(i, j));
839 }
840
841 m_flow_change.add(1.0, m_flow_change_incremental);
842 m_input_change.add(dt, surface_input_rate);
843 m_input_change.add(dt, basal_melt_rate);
844 }
845
846 //! Update the model state variables W and Wtill by applying the subglacial hydrology model equations.
847 /*!
848 Runs the hydrology model from time t to time t + dt. Here [t, dt]
849 is generally on the order of months to years. This hydrology model will take its
850 own shorter time steps, perhaps hours to weeks.
851
852 To update W = `bwat` we call update_W(), and to update Wtill = `tillwat` we
853 call update_Wtill().
854 */
855 void Routing::update_impl(double t, double dt, const Inputs& inputs) {
856
857 ice_bottom_surface(*inputs.geometry, m_bottom_surface);
858
859 double
860 ht = t,
861 hdt = 0.0;
862
863 const double
864 t_final = t + dt,
865 dt_max = m_config->get_number("hydrology.maximum_time_step", "seconds");
866
867 m_Qstag_average.set(0.0);
868
869 // make sure W has valid ghosts before starting hydrology steps
870 m_W.update_ghosts();
871
872 unsigned int step_counter = 0;
873 for (; ht < t_final; ht += hdt) {
874 step_counter++;
875
876 #if (Pism_DEBUG==1)
877 double huge_number = 1e6;
878 check_bounds(m_W, huge_number);
879
880 check_bounds(m_Wtill, m_config->get_number("hydrology.tillwat_max"));
881 #endif
882
883 // updates ghosts of m_Wstag
884 water_thickness_staggered(m_W,
885 inputs.geometry->cell_type,
886 m_Wstag);
887
888 double maxKW = 0.0;
889 // updates ghosts of m_Kstag
890 m_grid->ctx()->profiling().begin("routing_conductivity");
891 compute_conductivity(m_Wstag,
892 subglacial_water_pressure(),
893 m_bottom_surface,
894 m_Kstag, maxKW);
895 m_grid->ctx()->profiling().end("routing_conductivity");
896
897 // ghosts of m_Vstag are not updated
898 m_grid->ctx()->profiling().begin("routing_velocity");
899 compute_velocity(m_Wstag,
900 subglacial_water_pressure(),
901 m_bottom_surface,
902 m_Kstag,
903 inputs.no_model_mask,
904 m_Vstag);
905 m_grid->ctx()->profiling().end("routing_velocity");
906
907 // to get Q, W needs valid ghosts (ghosts of m_Vstag are not used)
908 // updates ghosts of m_Qstag
909 m_grid->ctx()->profiling().begin("routing_flux");
910 advective_fluxes(m_Vstag, m_W, m_Qstag);
911 m_grid->ctx()->profiling().end("routing_flux");
912
913 m_Qstag_average.add(hdt, m_Qstag);
914
915 {
916 const double
917 dt_cfl = max_timestep_W_cfl(),
918 dt_diff_w = max_timestep_W_diff(maxKW);
919
920 hdt = std::min(t_final - ht, dt_max);
921 hdt = std::min(hdt, dt_cfl);
922 hdt = std::min(hdt, dt_diff_w);
923 }
924
925 m_log->message(3, " hydrology step %05d, dt = %f s\n", step_counter, hdt);
926
927 // update Wtillnew from Wtill and input_rate
928 {
929 m_grid->ctx()->profiling().begin("routing_Wtill");
930 update_Wtill(hdt,
931 m_Wtill,
932 m_surface_input_rate,
933 m_basal_melt_rate,
934 m_Wtillnew);
935 // remove water in ice-free areas and account for changes
936 enforce_bounds(inputs.geometry->cell_type,
937 inputs.no_model_mask,
938 0.0, // do not limit maximum thickness
939 m_Wtillnew,
940 m_grounded_margin_change,
941 m_grounding_line_change,
942 m_conservation_error_change,
943 m_no_model_mask_change);
944 m_grid->ctx()->profiling().end("routing_Wtill");
945 }
946
947 // update Wnew from W, Wtill, Wtillnew, Wstag, Q, input_rate
948 // uses ghosts of m_W, m_Wstag, m_Qstag, m_Kstag
949 {
950 m_grid->ctx()->profiling().begin("routing_W");
951 update_W(hdt,
952 m_surface_input_rate,
953 m_basal_melt_rate,
954 m_W, m_Wstag,
955 m_Wtill, m_Wtillnew,
956 m_Kstag, m_Qstag,
957 m_Wnew);
958 // remove water in ice-free areas and account for changes
959 enforce_bounds(inputs.geometry->cell_type,
960 inputs.no_model_mask,
961 0.0, // do not limit maximum thickness
962 m_Wnew,
963 m_grounded_margin_change,
964 m_grounding_line_change,
965 m_conservation_error_change,
966 m_no_model_mask_change);
967
968 // transfer new into old (updates ghosts of m_W)
969 m_W.copy_from(m_Wnew);
970 m_grid->ctx()->profiling().end("routing_W");
971 }
972
973 // m_Wtill has no ghosts
974 m_Wtill.copy_from(m_Wtillnew);
975 } // end of the time-stepping loop
976
977 staggered_to_regular(inputs.geometry->cell_type, m_Qstag_average,
978 m_config->get_flag("hydrology.routing.include_floating_ice"),
979 m_Q);
980 m_Q.scale(1.0 / dt);
981
982 m_log->message(2,
983 " took %d hydrology sub-steps with average dt = %.6f years (%.3f s or %.3f hours)\n",
984 step_counter,
985 units::convert(m_sys, dt / step_counter, "seconds", "years"),
986 dt / step_counter,
987 (dt / step_counter) / 3600.0);
988 }
989
990 std::map<std::string, Diagnostic::Ptr> Routing::diagnostics_impl() const {
991 using namespace diagnostics;
992
993 DiagnosticList result = {
994 {"bwatvel", Diagnostic::Ptr(new BasalWaterVelocity(this))},
995 {"bwp", Diagnostic::Ptr(new BasalWaterPressure(this))},
996 {"bwprel", Diagnostic::Ptr(new RelativeBasalWaterPressure(this))},
997 {"effbwp", Diagnostic::Ptr(new EffectiveBasalWaterPressure(this))},
998 {"wallmelt", Diagnostic::Ptr(new WallMelt(this))},
999 {"hydraulic_potential", Diagnostic::Ptr(new HydraulicPotential(this))},
1000 };
1001 return combine(result, Hydrology::diagnostics_impl());
1002 }
1003
1004 std::map<std::string, TSDiagnostic::Ptr> Routing::ts_diagnostics_impl() const {
1005 std::map<std::string, TSDiagnostic::Ptr> result = {
1006 // FIXME: add mass-conservation diagnostics
1007 };
1008 return result;
1009 }
1010
1011 } // end of namespace hydrology
1012 } // end of namespace pism