tHydrology.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
---
tHydrology.hh (8507B)
---
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 _PISMHYDROLOGY_H_
20 #define _PISMHYDROLOGY_H_
21
22 #include "pism/util/iceModelVec.hh"
23 #include "pism/util/Component.hh"
24 #include "pism/util/iceModelVec2T.hh"
25
26 namespace pism {
27
28 class IceModelVec2CellType;
29 class Geometry;
30
31 //! @brief Sub-glacial hydrology models and related diagnostics.
32 namespace hydrology {
33
34 class Inputs {
35 public:
36 Inputs();
37
38 // modeling domain (set to NULL in whole-ice-sheet configurations)
39 const IceModelVec2Int *no_model_mask;
40 // geometry
41 const Geometry *geometry;
42 // hydrological inputs
43 const IceModelVec2S *surface_input_rate;
44 const IceModelVec2S *basal_melt_rate;
45 const IceModelVec2S *ice_sliding_speed;
46 };
47
48 //! \brief The PISM subglacial hydrology model interface.
49 /*!
50 This is a virtual base class.
51
52 The purpose of this class and its derived classes is to provide
53 \code
54 subglacial_water_thickness()
55 subglacial_water_pressure()
56 till_water_thickness()
57 \endcode
58 These correspond to state variables \f$W\f$, \f$P\f$, and \f$W_{\text{till}}\f$
59 in [\ref BuelervanPeltDRAFT], though not all derived classes of Hydrology
60 have all of them as state variables.
61
62 Additional modeled fields, for diagnostic purposes, are
63 \code
64 overburden_pressure(IceModelVec2S &result)
65 wall_melt(IceModelVec2S &result)
66 \endcode
67
68 This interface is appropriate to subglacial hydrology models which track a
69 two-dimensional water layer with a well-defined thickness and pressure at each
70 map-plane location. The methods subglacial_water_thickness() and
71 subglacial_water_pressure() return amount and pressure of the *transportable*
72 water, that is, the subglacial water which moves along a modeled hydraulic
73 head gradient, in contrast to the water stored in the till.
74
75 The transportable water moves through a subglacial morphology which is not
76 determined in this base class.
77
78 The Hydrology models have separate, but potentially-coupled, water
79 which is held in local till storage. Thus the
80 transportable water (bwat) and till water (tillwat) thicknesses are different.
81 Published models with till storage include [\ref BBssasliding, \ref SchoofTill,
82 \ref TrufferEchelmeyerHarrison2001, \ref Tulaczyketal2000b, \ref vanderWeletal2013].
83
84 The till water thickness is can be used, via the theory of
85 [\ref Tulaczyketal2000], to compute an effective pressure for the water in the
86 pore spaces of the till, which can then be used by the Mohr-Coulomb criterion
87 to provide a yield stress. Class MohrCoulombYieldStress does this
88 calculation. Here in Hydrology only the till water thickness tillwat is
89 computed.
90
91 Hydrology is a timestepping component. Because of the
92 short physical timescales associated to liquid water moving under a glacier,
93 Hydrology (and derived) classes generally take many substeps in PISM's major
94 ice dynamics time steps. Thus when an update() method in a Hydrology
95 class is called it will advance its internal time to the new goal t+dt
96 using its own internal time steps.
97
98 Generally Hydrology classes use the ice geometry, the basal melt
99 rate, and the basal sliding velocity in determining the evolution of the
100 hydrology state variables. Note that the basal melt rate is an
101 energy-conservation-derived field and the basal-sliding velocity is derived
102 from the solution of a stress balance. The basal melt rate and
103 sliding velocity fields therefore generally come from IceModel and
104 StressBalance, respectively.
105
106 Additional, time-dependent and spatially-variable water input to the basal
107 layer, taken directly from a file, is possible too.
108
109 Ice geometry and energy fields are normally treated as constant in time
110 during the update() call for the interval [t,t+dt]. Thus the coupling is
111 one-way during the update() call.
112 */
113 class Hydrology : public Component {
114 public:
115 Hydrology(IceGrid::ConstPtr g);
116 virtual ~Hydrology();
117
118 void restart(const File &input_file, int record);
119
120 void bootstrap(const File &input_file,
121 const IceModelVec2S &ice_thickness);
122
123 void init(const IceModelVec2S &W_till,
124 const IceModelVec2S &W,
125 const IceModelVec2S &P);
126
127 void update(double t, double dt, const Inputs& inputs);
128
129 const IceModelVec2S& till_water_thickness() const;
130 const IceModelVec2S& subglacial_water_thickness() const;
131 const IceModelVec2S& overburden_pressure() const;
132 const IceModelVec2S& surface_input_rate() const;
133 const IceModelVec2V& flux() const;
134
135 const IceModelVec2S& mass_change() const;
136 const IceModelVec2S& mass_change_at_grounded_margin() const;
137 const IceModelVec2S& mass_change_at_grounding_line() const;
138 const IceModelVec2S& mass_change_at_domain_boundary() const;
139 const IceModelVec2S& mass_change_due_to_conservation_error() const;
140 const IceModelVec2S& mass_change_due_to_input() const;
141 const IceModelVec2S& mass_change_due_to_lateral_flow() const;
142
143 protected:
144 virtual void restart_impl(const File &input_file, int record);
145
146 virtual void bootstrap_impl(const File &input_file,
147 const IceModelVec2S &ice_thickness);
148
149 virtual void init_impl(const IceModelVec2S &W_till,
150 const IceModelVec2S &W,
151 const IceModelVec2S &P);
152
153 virtual void update_impl(double t, double dt, const Inputs& inputs) = 0;
154 virtual std::map<std::string, Diagnostic::Ptr> diagnostics_impl() const;
155
156 virtual void define_model_state_impl(const File &output) const;
157 virtual void write_model_state_impl(const File &output) const;
158
159 void compute_overburden_pressure(const IceModelVec2S &ice_thickness,
160 IceModelVec2S &result) const;
161
162 void compute_surface_input_rate(const IceModelVec2CellType &mask,
163 const IceModelVec2S *surface_input_rate,
164 IceModelVec2S &result);
165
166 void compute_basal_melt_rate(const IceModelVec2CellType &mask,
167 const IceModelVec2S &basal_melt_rate,
168 IceModelVec2S &result);
169 protected:
170 // water flux on the regular grid
171 IceModelVec2V m_Q;
172
173 //! effective thickness of basal water stored in till
174 IceModelVec2S m_Wtill;
175
176 //! effective thickness of transportable basal water
177 IceModelVec2S m_W;
178
179 //! overburden pressure
180 IceModelVec2S m_Pover;
181
182 // surface input rate
183 IceModelVec2S m_surface_input_rate;
184
185 // input rate due to basal melt
186 IceModelVec2S m_basal_melt_rate;
187
188 // change due to flow for the current hydrology time step
189 IceModelVec2S m_flow_change_incremental;
190
191 // changes in water thickness
192 //
193 // these quantities are re-set to zero at the beginning of the PISM time step
194 IceModelVec2S m_conservation_error_change;
195 IceModelVec2S m_grounded_margin_change;
196 IceModelVec2S m_grounding_line_change;
197 IceModelVec2S m_input_change;
198 IceModelVec2S m_no_model_mask_change;
199 IceModelVec2S m_total_change;
200 IceModelVec2S m_flow_change;
201
202 // when we update the water amounts, careful mass accounting at the boundary
203 // is needed
204 void enforce_bounds(const IceModelVec2CellType &cell_type,
205 const IceModelVec2Int *no_model_mask,
206 double max_thickness,
207 IceModelVec2S &water_thickness,
208 IceModelVec2S &grounded_margin_change,
209 IceModelVec2S &grounding_line_change,
210 IceModelVec2S &conservation_error_change,
211 IceModelVec2S &no_model_mask_change);
212 private:
213 virtual void initialization_message() const = 0;
214 };
215
216 void check_bounds(const IceModelVec2S& W, double W_max);
217
218 } // end of namespace hydrology
219 } // end of namespace pism
220
221 #endif /* _PISMHYDROLOGY_H_ */