tEnthalpyModel.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
---
tEnthalpyModel.cc (15541B)
---
1 /* Copyright (C) 2016, 2017, 2018 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
20 #include "EnthalpyModel.hh"
21
22 #include "DrainageCalculator.hh"
23 #include "pism/util/EnthalpyConverter.hh"
24 #include "pism/energy/enthSystem.hh"
25 #include "pism/util/IceModelVec2CellType.hh"
26 #include "pism/util/io/File.hh"
27 #include "utilities.hh"
28 #include "pism/util/pism_utilities.hh"
29
30 namespace pism {
31 namespace energy {
32
33 EnthalpyModel::EnthalpyModel(IceGrid::ConstPtr grid,
34 stressbalance::StressBalance *stress_balance)
35 : EnergyModel(grid, stress_balance) {
36 // empty
37 }
38
39 void EnthalpyModel::restart_impl(const File &input_file, int record) {
40
41 m_log->message(2, "* Restarting the enthalpy-based energy balance model from %s...\n",
42 input_file.filename().c_str());
43
44 m_basal_melt_rate.read(input_file, record);
45 init_enthalpy(input_file, false, record);
46
47 regrid("Energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
48 regrid_enthalpy();
49 }
50
51 void EnthalpyModel::bootstrap_impl(const File &input_file,
52 const IceModelVec2S &ice_thickness,
53 const IceModelVec2S &surface_temperature,
54 const IceModelVec2S &climatic_mass_balance,
55 const IceModelVec2S &basal_heat_flux) {
56
57 m_log->message(2, "* Bootstrapping the enthalpy-based energy balance model from %s...\n",
58 input_file.filename().c_str());
59
60 m_basal_melt_rate.regrid(input_file, OPTIONAL,
61 m_config->get_number("bootstrapping.defaults.bmelt"));
62
63 regrid("Energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
64
65 int enthalpy_revision = m_ice_enthalpy.state_counter();
66 regrid_enthalpy();
67
68 if (enthalpy_revision == m_ice_enthalpy.state_counter()) {
69 bootstrap_ice_enthalpy(ice_thickness, surface_temperature, climatic_mass_balance,
70 basal_heat_flux, m_ice_enthalpy);
71 }
72 }
73
74 void EnthalpyModel::initialize_impl(const IceModelVec2S &basal_melt_rate,
75 const IceModelVec2S &ice_thickness,
76 const IceModelVec2S &surface_temperature,
77 const IceModelVec2S &climatic_mass_balance,
78 const IceModelVec2S &basal_heat_flux) {
79
80 m_log->message(2, "* Bootstrapping the enthalpy-based energy balance model...\n");
81
82 m_basal_melt_rate.copy_from(basal_melt_rate);
83
84 regrid("Energy balance model", m_basal_melt_rate, REGRID_WITHOUT_REGRID_VARS);
85
86 int enthalpy_revision = m_ice_enthalpy.state_counter();
87 regrid_enthalpy();
88
89 if (enthalpy_revision == m_ice_enthalpy.state_counter()) {
90 bootstrap_ice_enthalpy(ice_thickness, surface_temperature, climatic_mass_balance,
91 basal_heat_flux, m_ice_enthalpy);
92 }
93 }
94
95 //! Update ice enthalpy field based on conservation of energy.
96 /*!
97 This method is documented by the page \ref bombproofenth and by [\ref
98 AschwandenBuelerKhroulevBlatter].
99
100 This method updates IceModelVec3 m_work and IceModelVec2S basal_melt_rate.
101 No communication of ghosts is done for any of these fields.
102
103 We use an instance of enthSystemCtx.
104
105 Regarding drainage, see [\ref AschwandenBuelerKhroulevBlatter] and references therein.
106 */
107
108 void EnthalpyModel::update_impl(double t, double dt, const Inputs &inputs) {
109 // current time does not matter here
110 (void) t;
111
112 EnthalpyConverter::Ptr EC = m_grid->ctx()->enthalpy_converter();
113
114 const double
115 ice_density = m_config->get_number("constants.ice.density"), // kg m-3
116 bulgeEnthMax = m_config->get_number("energy.enthalpy.cold_bulge_max"), // J kg-1
117 target_water_fraction = m_config->get_number("energy.drainage_target_water_fraction");
118
119 energy::DrainageCalculator dc(*m_config);
120
121 inputs.check();
122
123 // give them names that are a bit shorter...
124 const IceModelVec3
125 &strain_heating3 = *inputs.volumetric_heating_rate,
126 &u3 = *inputs.u3,
127 &v3 = *inputs.v3,
128 &w3 = *inputs.w3;
129
130 const IceModelVec2CellType &cell_type = *inputs.cell_type;
131
132 const IceModelVec2S
133 &basal_frictional_heating = *inputs.basal_frictional_heating,
134 &basal_heat_flux = *inputs.basal_heat_flux,
135 &ice_thickness = *inputs.ice_thickness,
136 &surface_liquid_fraction = *inputs.surface_liquid_fraction,
137 &shelf_base_temp = *inputs.shelf_base_temp,
138 &ice_surface_temp = *inputs.surface_temp,
139 &till_water_thickness = *inputs.till_water_thickness;
140
141 energy::enthSystemCtx system(m_grid->z(), "energy.enthalpy", m_grid->dx(), m_grid->dy(), dt,
142 *m_config, m_ice_enthalpy, u3, v3, w3, strain_heating3, EC);
143
144 const size_t Mz_fine = system.z().size();
145 const double dz = system.dz();
146 std::vector<double> Enthnew(Mz_fine); // new enthalpy in column
147
148 IceModelVec::AccessList list{&ice_surface_temp, &shelf_base_temp, &surface_liquid_fraction,
149 &ice_thickness, &basal_frictional_heating, &basal_heat_flux, &till_water_thickness,
150 &cell_type, &u3, &v3, &w3, &strain_heating3, &m_basal_melt_rate, &m_ice_enthalpy,
151 &m_work};
152
153 double margin_threshold = m_config->get_number("energy.margin_ice_thickness_limit");
154
155 unsigned int liquifiedCount = 0;
156
157 ParallelSection loop(m_grid->com);
158 try {
159 for (Points pt(*m_grid); pt; pt.next()) {
160 const int i = pt.i(), j = pt.j();
161
162 const double H = ice_thickness(i, j);
163
164 system.init(i, j,
165 marginal(ice_thickness, i, j, margin_threshold),
166 H);
167
168 // enthalpy and pressures at top of ice
169 const double
170 depth_ks = H - system.ks() * dz,
171 p_ks = EC->pressure(depth_ks); // FIXME issue #15
172
173 const double Enth_ks = EC->enthalpy_permissive(ice_surface_temp(i, j),
174 surface_liquid_fraction(i, j), p_ks);
175
176 const bool ice_free_column = (system.ks() == 0);
177
178 // deal completely with columns with no ice; enthalpy and basal_melt_rate need setting
179 if (ice_free_column) {
180 m_work.set_column(i, j, Enth_ks);
181 // The floating basal melt rate will be set later; cover this
182 // case and set to zero for now. Also, there is no basal melt
183 // rate on ice free land and ice free ocean
184 m_basal_melt_rate(i, j) = 0.0;
185 continue;
186 } // end of if (ice_free_column)
187
188 if (system.lambda() < 1.0) {
189 m_stats.reduced_accuracy_counter += 1; // count columns with lambda < 1
190 }
191
192 const bool
193 is_floating = cell_type.ocean(i, j),
194 base_is_warm = system.Enth(0) >= system.Enth_s(0),
195 above_base_is_warm = system.Enth(1) >= system.Enth_s(1);
196
197 // set boundary conditions and update enthalpy
198 {
199 system.set_surface_dirichlet_bc(Enth_ks);
200
201 // determine lowest-level equation at bottom of ice; see
202 // decision chart in the source code browser and page
203 // documenting BOMBPROOF
204 if (is_floating) {
205 // floating base: Dirichlet application of known temperature from ocean
206 // coupler; assumes base of ice shelf has zero liquid fraction
207 double Enth0 = EC->enthalpy_permissive(shelf_base_temp(i, j), 0.0, EC->pressure(H));
208
209 system.set_basal_dirichlet_bc(Enth0);
210 } else {
211 // grounded ice warm and wet
212 if (base_is_warm && (till_water_thickness(i, j) > 0.0)) {
213 if (above_base_is_warm) {
214 // temperate layer at base (Neumann) case: q . n = 0 (K0 grad E . n = 0)
215 system.set_basal_heat_flux(0.0);
216 } else {
217 // only the base is warm: E = E_s(p) (Dirichlet)
218 // ( Assumes ice has zero liquid fraction. Is this a valid assumption here?
219 system.set_basal_dirichlet_bc(system.Enth_s(0));
220 }
221 } else {
222 // (Neumann) case: q . n = q_lith . n + F_b
223 // a) cold and dry base, or
224 // b) base that is still warm from the last time step, but without basal water
225 system.set_basal_heat_flux(basal_heat_flux(i, j) + basal_frictional_heating(i, j));
226 }
227 }
228
229 // solve the system
230 system.solve(Enthnew);
231
232 }
233
234 // post-process (drainage and bulge-limiting)
235 double Hdrainedtotal = 0.0;
236 double Hfrozen = 0.0;
237 {
238 // drain ice segments by mechanism in [\ref AschwandenBuelerKhroulevBlatter],
239 // using DrainageCalculator dc
240 for (unsigned int k=0; k < system.ks(); k++) {
241 if (Enthnew[k] > system.Enth_s(k)) { // avoid doing any more work if cold
242
243 const double
244 depth = H - k * dz,
245 p = EC->pressure(depth), // FIXME issue #15
246 T_m = EC->melting_temperature(p),
247 L = EC->L(T_m);
248
249 if (Enthnew[k] >= system.Enth_s(k) + 0.5 * L) {
250 liquifiedCount++; // count these rare events...
251 Enthnew[k] = system.Enth_s(k) + 0.5 * L; // but lose the energy
252 }
253
254 double omega = EC->water_fraction(Enthnew[k], p);
255
256 if (omega > target_water_fraction) {
257 double fractiondrained = dc.get_drainage_rate(omega) * dt; // pure number
258
259 fractiondrained = std::min(fractiondrained,
260 omega - target_water_fraction);
261 Hdrainedtotal += fractiondrained * dz; // always a positive contribution
262 Enthnew[k] -= fractiondrained * L;
263 }
264 }
265 }
266
267 // apply bulge limiter
268 const double lowerEnthLimit = Enth_ks - bulgeEnthMax;
269 for (unsigned int k=0; k < system.ks(); k++) {
270 if (Enthnew[k] < lowerEnthLimit) {
271 // Count grid points which have very large cold limit advection bulge... enthalpy not
272 // too low.
273 m_stats.bulge_counter += 1;
274 Enthnew[k] = lowerEnthLimit;
275 }
276 }
277
278 // if there is subglacial water, don't allow ice base enthalpy to be below
279 // pressure-melting; that is, assume subglacial water is at the pressure-
280 // melting temperature and enforce continuity of temperature
281 {
282 if (Enthnew[0] < system.Enth_s(0) && till_water_thickness(i,j) > 0.0) {
283 const double E_difference = system.Enth_s(0) - Enthnew[0];
284
285 const double depth = H,
286 pressure = EC->pressure(depth),
287 T_m = EC->melting_temperature(pressure);
288
289 Enthnew[0] = system.Enth_s(0);
290 // This adjustment creates energy out of nothing. We will
291 // freeze some basal water, subtracting an equal amount of
292 // energy, to make up for it.
293 //
294 // Note that [E_difference] = J/kg, so
295 //
296 // U_difference = E_difference * ice_density * dx * dy * (0.5*dz)
297 //
298 // is the amount of energy created (we changed enthalpy of
299 // a block of ice with the volume equal to
300 // dx*dy*(0.5*dz); note that the control volume
301 // corresponding to the grid point at the base of the
302 // column has thickness 0.5*dz, not dz).
303 //
304 // Also, [L] = J/kg, so
305 //
306 // U_freeze_on = L * ice_density * dx * dy * Hfrozen,
307 //
308 // is the amount of energy created by freezing a water
309 // layer of thickness Hfrozen (using units of ice
310 // equivalent thickness).
311 //
312 // Setting U_difference = U_freeze_on and solving for
313 // Hfrozen, we find the thickness of the basal water layer
314 // we need to freeze co restore energy conservation.
315
316 Hfrozen = E_difference * (0.5*dz) / EC->L(T_m);
317 }
318 }
319
320 } // end of post-processing
321
322 // compute basal melt rate
323 {
324 bool base_is_cold = (Enthnew[0] < system.Enth_s(0)) && (till_water_thickness(i,j) == 0.0);
325 // Determine melt rate, but only preliminarily because of
326 // drainage, from heat flux out of bedrock, heat flux into
327 // ice, and frictional heating
328 if (is_floating) {
329 // The floating basal melt rate will be set later; cover
330 // this case and set to zero for now. Note that
331 // Hdrainedtotal is discarded (the ocean model determines
332 // the basal melt).
333 m_basal_melt_rate(i, j) = 0.0;
334 } else {
335 if (base_is_cold) {
336 m_basal_melt_rate(i, j) = 0.0; // zero melt rate if cold base
337 } else {
338 const double
339 p_0 = EC->pressure(H),
340 p_1 = EC->pressure(H - dz), // FIXME issue #15
341 Tpmp_0 = EC->melting_temperature(p_0);
342
343 const bool k1_istemperate = EC->is_temperate(Enthnew[1], p_1); // level z = + \Delta z
344 double hf_up = 0.0;
345 if (k1_istemperate) {
346 const double
347 Tpmp_1 = EC->melting_temperature(p_1);
348
349 hf_up = -system.k_from_T(Tpmp_0) * (Tpmp_1 - Tpmp_0) / dz;
350 } else {
351 double T_0 = EC->temperature(Enthnew[0], p_0);
352 const double K_0 = system.k_from_T(T_0) / EC->c();
353
354 hf_up = -K_0 * (Enthnew[1] - Enthnew[0]) / dz;
355 }
356
357 // compute basal melt rate from flux balance:
358 //
359 // basal_melt_rate = - Mb / rho in [\ref AschwandenBuelerKhroulevBlatter];
360 //
361 // after we compute it we make sure there is no refreeze if
362 // there is no available basal water
363 m_basal_melt_rate(i, j) = (basal_frictional_heating(i, j) + basal_heat_flux(i, j) - hf_up) / (ice_density * EC->L(Tpmp_0));
364
365 if (till_water_thickness(i, j) <= 0 && m_basal_melt_rate(i, j) < 0) {
366 m_basal_melt_rate(i, j) = 0.0;
367 }
368 }
369
370 // Add drained water from the column to basal melt rate.
371 m_basal_melt_rate(i, j) += (Hdrainedtotal - Hfrozen) / dt;
372 } // end of the grounded case
373 } // end of the basal melt rate computation
374
375 system.fine_to_coarse(Enthnew, i, j, m_work);
376 }
377 } catch (...) {
378 loop.failed();
379 }
380 loop.check();
381
382 m_stats.liquified_ice_volume = ((double) liquifiedCount) * dz * m_grid->cell_area();
383 }
384
385 void EnthalpyModel::define_model_state_impl(const File &output) const {
386 m_ice_enthalpy.define(output);
387 m_basal_melt_rate.define(output);
388 }
389
390 void EnthalpyModel::write_model_state_impl(const File &output) const {
391 m_ice_enthalpy.write(output);
392 m_basal_melt_rate.write(output);
393 }
394
395 } // end of namespace energy
396 } // end of namespace pism