tMohrCoulombYieldStress.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
---
tMohrCoulombYieldStress.cc (18929B)
---
1 // Copyright (C) 2004--2020 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 "MohrCoulombYieldStress.hh"
20 #include "MohrCoulombPointwise.hh"
21
22 #include "pism/util/IceGrid.hh"
23 #include "pism/util/Mask.hh"
24 #include "pism/util/error_handling.hh"
25 #include "pism/util/io/File.hh"
26 #include "pism/util/MaxTimestep.hh"
27 #include "pism/util/pism_utilities.hh"
28 #include "pism/util/Time.hh"
29 #include "pism/util/IceModelVec2CellType.hh"
30 #include "pism/geometry/Geometry.hh"
31 #include "pism/coupler/util/options.hh" // ForcingOptions
32
33 namespace pism {
34
35 //! \file MohrCoulombYieldStress.cc Process model which computes pseudo-plastic yield stress for the subglacial layer.
36
37 /*! \file MohrCoulombYieldStress.cc
38 The output variable of this submodel is `tauc`, the pseudo-plastic yield stress
39 field that is used in the ShallowStressBalance objects. This quantity is
40 computed by the Mohr-Coulomb criterion [\ref SchoofTill], but using an empirical
41 relation between the amount of water in the till and the effective pressure
42 of the overlying glacier resting on the till [\ref Tulaczyketal2000].
43
44 The "dry" strength of the till is a state variable which is private to
45 the submodel, namely `tillphi`. Its initialization is nontrivial: either the
46 `-topg_to_phi` heuristic is used or inverse modeling can be used. (In the
47 latter case `tillphi` can be read-in at the beginning of the run. Currently
48 `tillphi` does not evolve during the run.)
49
50 The effective pressure is derived from the till (pore) water amount (the effective water
51 layer thickness). Then the effective pressure is combined with tillphi to compute an
52 updated `tauc` by the Mohr-Coulomb criterion.
53
54 This submodel is inactive in floating areas.
55 */
56
57
58 /*!
59 The pseudo-plastic till basal resistance model is governed by this power law
60 equation,
61 @f[ \tau_b = - \frac{\tau_c}{|\mathbf{U}|^{1-q} U_{\mathtt{th}}^q} \mathbf{U}, @f]
62 where @f$\tau_b=(\tau_{(b)x},\tau_{(b)y})@f$ is the basal shear stress and
63 @f$U=(u,v)@f$ is the sliding velocity.
64
65 We call the scalar field @f$\tau_c(t,x,y)@f$ the *yield stress* even when
66 the power @f$q@f$ is not zero; when that power is zero the formula describes
67 a plastic material with an actual yield stress. The constant
68 @f$U_{\mathtt{th}}@f$ is the *threshold speed*, and @f$q@f$ is the *pseudo*
69 *plasticity exponent*. The current class computes this yield stress field.
70 See also IceBasalResistancePlasticLaw::drag().
71
72 The strength of the saturated till material, the yield stress, is modeled by a
73 Mohr-Coulomb relation [\ref Paterson, \ref SchoofStream],
74 @f[ \tau_c = c_0 + (\tan \varphi) N_{till}, @f]
75 where @f$N_{till}@f$ is the effective pressure of the glacier on the mineral
76 till.
77
78 The determination of the till friction angle @f$\varphi(x,y)@f$ is important.
79 It is assumed in this default model to be a time-independent factor which
80 describes the strength of the unsaturated "dry" (mineral) till material. Thus
81 it is assumed to change more slowly than the till water pressure, and it follows
82 that it changes more slowly than the yield stress and the basal shear stress.
83
84 Option `-topg_to_phi` causes call to topg_to_phi() at the beginning of the run.
85 This determines the map of @f$\varphi(x,y)@f$. If this option is note given,
86 the current method leaves `tillphi` unchanged, and thus either in its
87 read-in-from-file state or with a default constant value from the config file.
88 */
89 MohrCoulombYieldStress::MohrCoulombYieldStress(IceGrid::ConstPtr grid)
90 : YieldStress(grid),
91 m_till_phi(m_grid, "tillphi", WITHOUT_GHOSTS) {
92
93 m_name = "Mohr-Coulomb yield stress model";
94
95 m_till_phi.set_attrs("model_state",
96 "friction angle for till under grounded ice sheet",
97 "degrees", "degrees", "", 0);
98 m_till_phi.set_time_independent(true);
99 // in this model; need not be time-independent in general
100
101 {
102 std::string hydrology_tillwat_max = "hydrology.tillwat_max";
103 bool till_is_present = m_config->get_number(hydrology_tillwat_max) > 0.0;
104
105 if (not till_is_present) {
106 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
107 "The Mohr-Coulomb yield stress model cannot be used without till.\n"
108 "Reset %s or choose a different yield stress model.",
109 hydrology_tillwat_max.c_str());
110 }
111 }
112
113 auto delta_file = m_config->get_string("basal_yield_stress.mohr_coulomb.delta.file");
114
115 if (not delta_file.empty()) {
116 ForcingOptions opt(*m_grid->ctx(), "basal_yield_stress.mohr_coulomb.delta");
117
118 unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
119 unsigned int evaluations_per_year = m_config->get_number("input.forcing.evaluations_per_year");
120 bool periodic = opt.period > 0;
121
122 File file(m_grid->com, opt.filename, PISM_NETCDF3, PISM_READONLY);
123
124 m_delta = IceModelVec2T::ForcingField(m_grid,
125 file,
126 "mohr_coulomb_delta",
127 "", // no standard name
128 buffer_size,
129 evaluations_per_year,
130 periodic, LINEAR);
131 m_delta->set_attrs("", "minimum effective pressure on till as a fraction of overburden pressure",
132 "1", "1", "", 0);
133 }
134 }
135
136 MohrCoulombYieldStress::~MohrCoulombYieldStress() {
137 // empty
138 }
139
140 void MohrCoulombYieldStress::restart_impl(const File &input_file, int record) {
141 m_basal_yield_stress.read(input_file, record);
142 m_till_phi.read(input_file, record);
143 }
144
145
146 //! Initialize the pseudo-plastic till mechanical model.
147 void MohrCoulombYieldStress::bootstrap_impl(const File &input_file,
148 const YieldStressInputs &inputs) {
149
150 auto tauc_to_phi_file = m_config->get_string("basal_yield_stress.mohr_coulomb.tauc_to_phi.file");
151
152 if (not tauc_to_phi_file.empty()) {
153 m_basal_yield_stress.regrid(tauc_to_phi_file, CRITICAL);
154
155 m_log->message(2,
156 " Will compute till friction angle (tillphi) as a function"
157 " of the yield stress (tauc)...\n");
158
159 till_friction_angle(m_basal_yield_stress,
160 *inputs.till_water_thickness,
161 inputs.geometry->ice_thickness,
162 inputs.geometry->cell_type,
163 m_till_phi);
164
165 } else if (m_config->get_flag("basal_yield_stress.mohr_coulomb.topg_to_phi.enabled")) {
166
167 m_log->message(2, " creating till friction angle map from bed elevation...\n");
168
169 if (input_file.find_variable(m_till_phi.metadata().get_name())) {
170 // till_phi is present in the input file
171 m_log->message(2,
172 "PISM WARNING: -topg_to_phi computation will override the '%s' field\n"
173 " present in the input file '%s'!\n",
174 m_till_phi.metadata().get_name().c_str(), input_file.filename().c_str());
175 }
176
177 till_friction_angle(inputs.geometry->bed_elevation, m_till_phi);
178
179 } else {
180 double till_phi_default = m_config->get_number("basal_yield_stress.mohr_coulomb.till_phi_default");
181 m_till_phi.regrid(input_file, OPTIONAL, till_phi_default);
182 }
183
184 finish_initialization(inputs);
185 }
186
187 void MohrCoulombYieldStress::init_impl(const YieldStressInputs &inputs) {
188 double till_phi_default = m_config->get_number("basal_yield_stress.mohr_coulomb.till_phi_default");
189 m_till_phi.set(till_phi_default);
190
191 finish_initialization(inputs);
192 }
193
194 /*!
195 * Finish initialization after bootstrapping or initializing using constants.
196 */
197 void MohrCoulombYieldStress::finish_initialization(const YieldStressInputs &inputs) {
198 // regrid if requested, regardless of how initialized
199 regrid(name(), m_till_phi);
200
201 if (m_delta) {
202 ForcingOptions opt(*m_grid->ctx(), "basal_yield_stress.mohr_coulomb.delta");
203
204 m_delta->init(opt.filename, opt.period, opt.reference_time);
205 }
206
207 // We use a short time step length because we can get away with it here, but we can
208 // probably do better...
209 this->update(inputs, m_grid->ctx()->time()->current(), 1.0 /* one second time step */);
210 }
211
212 MaxTimestep MohrCoulombYieldStress::max_timestep_impl(double t) const {
213 (void) t;
214
215 if (m_delta) {
216 auto dt = m_delta->max_timestep(t);
217
218 if (dt.finite()) {
219 return MaxTimestep(dt.value(), name());
220 }
221 }
222
223 return MaxTimestep(name());
224 }
225
226 void MohrCoulombYieldStress::set_till_friction_angle(const IceModelVec2S &input) {
227 m_till_phi.copy_from(input);
228 }
229
230 void MohrCoulombYieldStress::define_model_state_impl(const File &output) const {
231 m_basal_yield_stress.define(output);
232 m_till_phi.define(output);
233 }
234
235 void MohrCoulombYieldStress::write_model_state_impl(const File &output) const {
236 m_basal_yield_stress.write(output);
237 m_till_phi.write(output);
238 }
239
240 //! Update the till yield stress for use in the pseudo-plastic till basal stress
241 //! model. See also IceBasalResistancePlasticLaw.
242 /*!
243 Updates yield stress @f$ \tau_c @f$ based on modeled till water layer thickness
244 from a Hydrology object. We implement the Mohr-Coulomb criterion allowing
245 a (typically small) till cohesion @f$ c_0 @f$
246 and by expressing the coefficient as the tangent of a till friction angle
247 @f$ \varphi @f$ :
248 @f[ \tau_c = c_0 + (\tan \varphi) N_{till}. @f]
249 See [@ref Paterson] table 8.1 regarding values.
250
251 The effective pressure on the till is empirically-related
252 to the amount of water in the till. We use this formula derived from
253 [@ref Tulaczyketal2000] and documented in [@ref BuelervanPeltDRAFT]:
254
255 @f[ N_{till} = \min\left\{P_o, N_0 \left(\frac{\delta P_o}{N_0}\right)^s 10^{(e_0/C_c) (1 - s)}\right\} @f]
256
257 where @f$ s = W_{till} / W_{till}^{max} @f$, @f$ W_{till}^{max} @f$ =`hydrology_tillwat_max`,
258 @f$ \delta @f$ =`basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden`, @f$ P_o @f$ is the
259 overburden pressure, @f$ N_0 @f$ =`basal_yield_stress.mohr_coulomb.till_reference_effective_pressure` is a
260 reference effective pressure, @f$ e_0 @f$ =`basal_yield_stress.mohr_coulomb.till_reference_void_ratio` is the void ratio
261 at the reference effective pressure, and @f$ C_c @f$ =`basal_yield_stress.mohr_coulomb.till_compressibility_coefficient`
262 is the coefficient of compressibility of the till. Constants @f$ N_0, e_0, C_c @f$ are
263 found by [@ref Tulaczyketal2000] from laboratory experiments on samples of
264 till.
265
266 If `basal_yield_stress.add_transportable_water` is yes then @f$ s @f$ in the above formula
267 becomes @f$ s = (W + W_{till}) / W_{till}^{max} @f$,
268 that is, the water amount is the sum @f$ W+W_{till} @f$.
269 */
270 void MohrCoulombYieldStress::update_impl(const YieldStressInputs &inputs,
271 double t, double dt) {
272 (void) t;
273 (void) dt;
274
275 bool slippery_grounding_lines = m_config->get_flag("basal_yield_stress.slippery_grounding_lines"),
276 add_transportable_water = m_config->get_flag("basal_yield_stress.add_transportable_water");
277
278 const double
279 ice_density = m_config->get_number("constants.ice.density"),
280 standard_gravity = m_config->get_number("constants.standard_gravity");
281
282 const double high_tauc = m_config->get_number("basal_yield_stress.ice_free_bedrock"),
283 W_till_max = m_config->get_number("hydrology.tillwat_max"),
284 delta = m_config->get_number("basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden"),
285 tlftw = m_config->get_number("basal_yield_stress.mohr_coulomb.till_log_factor_transportable_water");
286
287 MohrCoulombPointwise mc(m_config);
288
289 const IceModelVec2S
290 &W_till = *inputs.till_water_thickness,
291 &W_subglacial = *inputs.subglacial_water_thickness,
292 &ice_thickness = inputs.geometry->ice_thickness;
293
294 const IceModelVec2CellType &cell_type = inputs.geometry->cell_type;
295 const IceModelVec2S &bed_topography = inputs.geometry->bed_elevation;
296 const IceModelVec2S &sea_level = inputs.geometry->sea_level_elevation;
297
298 IceModelVec::AccessList list{&W_till, &m_till_phi, &m_basal_yield_stress, &cell_type,
299 &bed_topography, &sea_level, &ice_thickness};
300
301 if (add_transportable_water) {
302 list.add(W_subglacial);
303 }
304
305 if (m_delta) {
306 m_delta->update(t, dt);
307 m_delta->average(t, dt);
308 list.add(*m_delta);
309 }
310
311 for (Points p(*m_grid); p; p.next()) {
312 const int i = p.i(), j = p.j();
313
314 if (cell_type.ice_free(i, j)) {
315 m_basal_yield_stress(i, j) = high_tauc; // large yield stress if ice-free
316 } else { // grounded and there is some ice
317
318 // user can ask that marine grounding lines get special treatment
319 double water = W_till(i,j); // usual case
320
321 if (slippery_grounding_lines and
322 bed_topography(i, j) <= sea_level(i, j) and
323 (cell_type.next_to_floating_ice(i, j) or cell_type.next_to_ice_free_ocean(i, j))) {
324 water = W_till_max;
325 } else if (add_transportable_water) {
326 water = W_till(i, j) + tlftw * log(1.0 + W_subglacial(i, j) / tlftw);
327 }
328
329 double P_overburden = ice_density * standard_gravity * ice_thickness(i, j);
330
331 m_basal_yield_stress(i, j) = mc.yield_stress(m_delta ? (*m_delta)(i, j) : delta,
332 P_overburden, water, m_till_phi(i, j));
333 }
334 }
335
336 m_basal_yield_stress.update_ghosts();
337 }
338
339 //! Computes the till friction angle phi as a piecewise linear function of bed elevation, according to user options.
340 /*!
341 Computes the till friction angle \f$\phi(x,y)\f$ at a location as the following
342 increasing, piecewise-linear function of the bed elevation \f$b(x,y)\f$. Let
343 \f[ M = (\phi_{\text{max}} - \phi_{\text{min}}) / (b_{\text{max}} - b_{\text{min}}) \f]
344 be the slope of the nontrivial part. Then
345 \f[ \phi(x,y) = \begin{cases}
346 \phi_{\text{min}}, & b(x,y) \le b_{\text{min}}, \\
347 \phi_{\text{min}} + (b(x,y) - b_{\text{min}}) \,M,
348 & b_{\text{min}} < b(x,y) < b_{\text{max}}, \\
349 \phi_{\text{max}}, & b_{\text{max}} \le b(x,y), \end{cases} \f]
350 where \f$\phi_{\text{min}}=\f$`phi_min`, \f$\phi_{\text{max}}=\f$`phi_max`,
351 \f$b_{\text{min}}=\f$`topg_min`, \f$b_{\text{max}}=\f$`topg_max`.
352
353 The default values are vaguely suitable for Antarctica. See src/pism_config.cdl.
354 */
355 void MohrCoulombYieldStress::till_friction_angle(const IceModelVec2S &bed_topography,
356 IceModelVec2S &result) {
357
358 const double
359 phi_min = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.phi_min"),
360 phi_max = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.phi_max"),
361 topg_min = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.topg_min"),
362 topg_max = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.topg_max");
363
364 m_log->message(2,
365 " till friction angle (phi) is piecewise-linear function of bed elev (topg):\n"
366 " / %5.2f for topg < %.f\n"
367 " phi = | %5.2f + (topg - (%.f)) * (%.2f / %.f) for %.f < topg < %.f\n"
368 " \\ %5.2f for %.f < topg\n",
369 phi_min, topg_min,
370 phi_min, topg_min, phi_max - phi_min, topg_max - topg_min, topg_min, topg_max,
371 phi_max, topg_max);
372
373 if (phi_min >= phi_max) {
374 throw RuntimeError(PISM_ERROR_LOCATION,
375 "invalid -topg_to_phi arguments: phi_min < phi_max is required");
376 }
377
378 if (topg_min >= topg_max) {
379 throw RuntimeError(PISM_ERROR_LOCATION,
380 "invalid -topg_to_phi arguments: topg_min < topg_max is required");
381 }
382
383 const double slope = (phi_max - phi_min) / (topg_max - topg_min);
384
385 IceModelVec::AccessList list{&bed_topography, &result};
386
387 for (Points p(*m_grid); p; p.next()) {
388 const int i = p.i(), j = p.j();
389 const double bed = bed_topography(i, j);
390
391 if (bed <= topg_min) {
392 result(i, j) = phi_min;
393 } else if (bed >= topg_max) {
394 result(i, j) = phi_max;
395 } else {
396 result(i, j) = phi_min + (bed - topg_min) * slope;
397 }
398 }
399
400 // communicate ghosts so that the tauc computation can be performed locally
401 // (including ghosts of tauc, that is)
402 result.update_ghosts();
403 }
404
405 /*!
406 * Compute the till friction angle in grounded areas using available basal yield stress,
407 * till water thickness, and overburden pressure.
408 *
409 * This is the inverse of the formula used by `update_impl()`.
410 */
411 void MohrCoulombYieldStress::till_friction_angle(const IceModelVec2S &basal_yield_stress,
412 const IceModelVec2S &till_water_thickness,
413 const IceModelVec2S &ice_thickness,
414 const IceModelVec2CellType &cell_type,
415 IceModelVec2S &result) {
416
417 MohrCoulombPointwise mc(m_config);
418
419 double
420 ice_density = m_config->get_number("constants.ice.density"),
421 standard_gravity = m_config->get_number("constants.standard_gravity");
422
423 double
424 delta = m_config->get_number("basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden");
425
426 const IceModelVec2S
427 &W_till = till_water_thickness;
428
429 IceModelVec::AccessList list{&cell_type, &basal_yield_stress, &W_till, &ice_thickness, &result};
430
431 for (Points p(*m_grid); p; p.next()) {
432 const int i = p.i(), j = p.j();
433
434 if (cell_type.ocean(i, j)) {
435 // no change
436 } else if (cell_type.ice_free(i, j)) {
437 // no change
438 } else { // grounded and there is some ice
439 double P_overburden = ice_density * standard_gravity * ice_thickness(i, j);
440
441 result(i, j) = mc.till_friction_angle(delta, P_overburden, W_till(i, j), basal_yield_stress(i, j));
442 }
443 }
444
445 result.update_ghosts();
446 }
447
448 DiagnosticList MohrCoulombYieldStress::diagnostics_impl() const {
449 return combine({{"tillphi", Diagnostic::wrap(m_till_phi)}},
450 YieldStress::diagnostics_impl());
451 }
452
453 } // end of namespace pism