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