tRouting.hh - 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.hh (7735B)
---
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 #ifndef _ROUTING_H_
20 #define _ROUTING_H_
21
22 #include "Hydrology.hh"
23
24 namespace pism {
25
26
27 namespace hydrology {
28
29 //! \brief A subglacial hydrology model which assumes water pressure
30 //! equals overburden pressure.
31 /*!
32 This PISM hydrology model has lateral motion of subglacial water and which
33 conserves the water mass. Further documentation is in [\ref BuelervanPeltDRAFT].
34
35 The water velocity is along the steepest descent route for the hydraulic
36 potential. This potential is (mostly) a function of ice sheet geometry,
37 because the water pressure is set to the overburden pressure, a simplified but
38 well-established model [\ref Shreve1972]. However, the water layer thickness
39 is also a part of the hydraulic potential because it is actually the potential
40 of the *top* of the water layer.
41
42 This (essential) model has been used for finding locations of subglacial lakes
43 [\ref Siegertetal2007, \ref Livingstoneetal2013]. Subglacial lakes occur
44 at local minima of the hydraulic potential. If water builds up significantly
45 (e.g. thickness of 10s of meters or more) then in the model here the resulting
46 lakes diffuse instead of becoming infinitely deep. Thus we avoid delta
47 functions of water thickness at the minima of the hydraulic potential in this
48 well-posed model.
49
50 This model should generally be tested using static ice geometry first.
51
52 The state space includes both the till water effective thickness \f$W_{till}\f$,
53 which is in Hydrology, and the transportable water layer thickness \f$W\f$.
54
55 For more complete modeling where the water pressure is determined by a
56 physical model for the opening and closing of cavities, and where the state
57 space includes a nontrivial pressure variable, see hydrology::Distributed.
58
59 There is an option `-hydrology_null_strip` `X` which produces a strip of
60 `X` km around the edge of the computational domain. In that strip the water flow
61 velocity is set to zero. The water amount is also reset to zero at the end
62 of each time step in this strip (in an accounted way).
63
64 As noted this is the minimal model which has a lateral water flux. This flux is
65 \f[ \mathbf{q} = - K \nabla \psi = \mathbf{V} W - D \nabla W \f]
66 where \f$\psi\f$ is the hydraulic potential
67 \f[ \psi = P + \rho_w g (b + W). \f]
68 The generalized conductivity \f$K\f$ is nontrivial and it generally also
69 depends on the water thickness:
70 \f[ K = k W^{\alpha-1} |\nabla (P+\rho_w g b)|^{\beta-2}. \f]
71
72 This model contains enough information (enough modeled fields) so that we can compute
73 the wall melt generated by dissipating the gravitational potential energy in the moving,
74 presumably turbulent, subglacial water. If we suppose that this heat is dissipated
75 immediately as melt on the cavity/conduit walls then we get a formula for a wall melt
76 contribution. (This is in addition to the `basal_melt_rate_grounded` field coming from
77 conserving energy in the flowing ice.) See wall_melt(). At this time the wall melt is
78 diagnostic only and does not add to the water amount W; such an addition is generally
79 unstable.
80 */
81 class Routing : public Hydrology {
82 public:
83
84 Routing(IceGrid::ConstPtr g);
85 virtual ~Routing();
86
87 const IceModelVec2S& subglacial_water_pressure() const;
88 const IceModelVec2Stag& velocity_staggered() const;
89
90 protected:
91 virtual void restart_impl(const File &input_file, int record);
92
93 virtual void bootstrap_impl(const File &input_file,
94 const IceModelVec2S &ice_thickness);
95
96 virtual void init_impl(const IceModelVec2S &W_till,
97 const IceModelVec2S &W,
98 const IceModelVec2S &P);
99
100 virtual void update_impl(double t, double dt, const Inputs& inputs);
101
102 virtual std::map<std::string, Diagnostic::Ptr> diagnostics_impl() const;
103 virtual std::map<std::string, TSDiagnostic::Ptr> ts_diagnostics_impl() const;
104
105 virtual void define_model_state_impl(const File &output) const;
106 virtual void write_model_state_impl(const File &output) const;
107
108 double max_timestep_W_diff(double KW_max) const;
109 double max_timestep_W_cfl() const;
110 protected:
111
112 // edge-centered (staggered) advection flux
113 IceModelVec2Stag m_Qstag;
114
115 IceModelVec2Stag m_Qstag_average;
116
117 // edge-centered (staggered) water velocity
118 IceModelVec2Stag m_Vstag;
119
120 // edge-centered (staggered) W values (averaged from regular)
121 IceModelVec2Stag m_Wstag;
122
123 // edge-centered (staggered) values of nonlinear conductivity
124 IceModelVec2Stag m_Kstag;
125
126 // work space
127 IceModelVec2S m_Wnew, m_Wtillnew;
128
129 // ghosted temporary storage; modified in compute_conductivity and compute_velocity
130 mutable IceModelVec2S m_R;
131
132 double m_dx, m_dy;
133 double m_rg;
134
135 IceModelVec2S m_bottom_surface;
136
137 void water_thickness_staggered(const IceModelVec2S &W,
138 const IceModelVec2CellType &mask,
139 IceModelVec2Stag &result);
140
141 void compute_conductivity(const IceModelVec2Stag &W,
142 const IceModelVec2S &P,
143 const IceModelVec2S &bed,
144 IceModelVec2Stag &result,
145 double &maxKW) const;
146
147 void compute_velocity(const IceModelVec2Stag &W,
148 const IceModelVec2S &P,
149 const IceModelVec2S &bed,
150 const IceModelVec2Stag &K,
151 const IceModelVec2Int *no_model_mask,
152 IceModelVec2Stag &result) const;
153
154 void advective_fluxes(const IceModelVec2Stag &V,
155 const IceModelVec2S &W,
156 IceModelVec2Stag &result) const;
157
158 void W_change_due_to_flow(double dt,
159 const IceModelVec2S &W,
160 const IceModelVec2Stag &Wstag,
161 const IceModelVec2Stag &K,
162 const IceModelVec2Stag &Q,
163 IceModelVec2S &result);
164 void update_W(double dt,
165 const IceModelVec2S &surface_input_rate,
166 const IceModelVec2S &basal_melt_rate,
167 const IceModelVec2S &W,
168 const IceModelVec2Stag &Wstag,
169 const IceModelVec2S &Wtill,
170 const IceModelVec2S &Wtill_new,
171 const IceModelVec2Stag &K,
172 const IceModelVec2Stag &Q,
173 IceModelVec2S &W_new);
174
175 void update_Wtill(double dt,
176 const IceModelVec2S &Wtill,
177 const IceModelVec2S &surface_input_rate,
178 const IceModelVec2S &basal_melt_rate,
179 IceModelVec2S &Wtill_new);
180
181 private:
182 virtual void initialization_message() const;
183 };
184
185 void wall_melt(const Routing &model,
186 const IceModelVec2S &bed_elevation,
187 IceModelVec2S &result);
188
189 } // end of namespace hydrology
190 } // end of namespace pism
191
192 #endif /* _ROUTING_H_ */