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