PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
MohrCoulombYieldStress.cc
Go to the documentation of this file.
1 // Copyright (C) 2004--2023 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 "pism/basalstrength/MohrCoulombYieldStress.hh"
20 #include "pism/basalstrength/MohrCoulombPointwise.hh"
21 
22 #include "pism/util/Grid.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/array/CellType.hh"
30 #include "pism/util/array/Forcing.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 */
90 MohrCoulombYieldStress::MohrCoulombYieldStress(std::shared_ptr<const Grid> grid)
91  : YieldStress(grid),
92  m_till_phi(m_grid, "tillphi") {
93 
94  m_name = "Mohr-Coulomb yield stress model";
95 
97  .long_name("friction angle for till under grounded ice sheet")
98  .units("degrees")
99  .set_time_independent(true);
100 
101  // in this model; need not be time-independent in general
102 
103  {
104  std::string hydrology_tillwat_max = "hydrology.tillwat_max";
105  bool till_is_present = m_config->get_number(hydrology_tillwat_max) > 0.0;
106 
107  if (not till_is_present) {
109  "The Mohr-Coulomb yield stress model cannot be used without till.\n"
110  "Reset %s or choose a different yield stress model.",
111  hydrology_tillwat_max.c_str());
112  }
113  }
114 
115  auto delta_file = m_config->get_string("basal_yield_stress.mohr_coulomb.delta.file");
116 
117  if (not delta_file.empty()) {
118  ForcingOptions opt(*m_grid->ctx(), "basal_yield_stress.mohr_coulomb.delta");
119 
120  unsigned int buffer_size = m_config->get_number("input.forcing.buffer_size");
121 
123 
124  m_delta = std::make_shared<array::Forcing>(m_grid,
125  file,
126  "mohr_coulomb_delta",
127  "", // no standard name
128  buffer_size,
129  opt.periodic, LINEAR);
130  m_delta->metadata()
131  .long_name("minimum effective pressure on till as a fraction of overburden pressure")
132  .units("1");
133  }
134 }
135 
136 void MohrCoulombYieldStress::restart_impl(const File &input_file, int record) {
137  m_basal_yield_stress.read(input_file, record);
138  m_till_phi.read(input_file, record);
139 }
140 
141 
142 //! Initialize the pseudo-plastic till mechanical model.
144  const YieldStressInputs &inputs) {
145 
146  auto tauc_to_phi_file = m_config->get_string("basal_yield_stress.mohr_coulomb.tauc_to_phi.file");
147 
148  if (not tauc_to_phi_file.empty()) {
149  m_basal_yield_stress.regrid(tauc_to_phi_file, io::Default::Nil());
150 
151  m_log->message(2,
152  " Will compute till friction angle (tillphi) as a function"
153  " of the yield stress (tauc)...\n");
154 
156  *inputs.till_water_thickness,
157  inputs.geometry->ice_thickness,
158  inputs.geometry->cell_type,
159  m_till_phi);
160 
161  } else if (m_config->get_flag("basal_yield_stress.mohr_coulomb.topg_to_phi.enabled")) {
162 
163  m_log->message(2, " creating till friction angle map from bed elevation...\n");
164 
165  if (input_file.find_variable(m_till_phi.metadata().get_name())) {
166  // till_phi is present in the input file
167  m_log->message(2,
168  "PISM WARNING: -topg_to_phi computation will override the '%s' field\n"
169  " present in the input file '%s'!\n",
170  m_till_phi.metadata().get_name().c_str(), input_file.filename().c_str());
171  }
172 
174 
175  } else {
176  double till_phi_default = m_config->get_number("basal_yield_stress.mohr_coulomb.till_phi_default");
177  m_till_phi.regrid(input_file, io::Default(till_phi_default));
178  }
179 
180  finish_initialization(inputs);
181 }
182 
184  double till_phi_default = m_config->get_number("basal_yield_stress.mohr_coulomb.till_phi_default");
185  m_till_phi.set(till_phi_default);
186 
187  finish_initialization(inputs);
188 }
189 
190 /*!
191  * Finish initialization after bootstrapping or initializing using constants.
192  */
194  // regrid if requested, regardless of how initialized
195  regrid(name(), m_till_phi);
196 
197  if (m_delta) {
198  ForcingOptions opt(*m_grid->ctx(), "basal_yield_stress.mohr_coulomb.delta");
199 
200  m_delta->init(opt.filename, opt.periodic);
201  }
202 
203  // We use a short time step length because we can get away with it here, but we can
204  // probably do better...
205  this->update(inputs, time().current(), 1.0 /* one second time step */);
206 }
207 
209  (void) t;
210 
211  if (m_delta) {
212  auto dt = m_delta->max_timestep(t);
213 
214  if (dt.finite()) {
215  return MaxTimestep(dt.value(), name());
216  }
217  }
218 
219  return MaxTimestep(name());
220 }
221 
223  m_till_phi.copy_from(input);
224 }
225 
229 }
230 
232  m_basal_yield_stress.write(output);
233  m_till_phi.write(output);
234 }
235 
236 //! Update the till yield stress for use in the pseudo-plastic till basal stress
237 //! model. See also IceBasalResistancePlasticLaw.
238 /*!
239 Updates yield stress @f$ \tau_c @f$ based on modeled till water layer thickness
240 from a Hydrology object. We implement the Mohr-Coulomb criterion allowing
241 a (typically small) till cohesion @f$ c_0 @f$
242 and by expressing the coefficient as the tangent of a till friction angle
243  @f$ \varphi @f$ :
244  @f[ \tau_c = c_0 + (\tan \varphi) N_{till}. @f]
245 See [@ref Paterson] table 8.1 regarding values.
246 
247 The effective pressure on the till is empirically-related
248 to the amount of water in the till. We use this formula derived from
249 [@ref Tulaczyketal2000] and documented in [@ref BuelervanPeltDRAFT]:
250 
251 @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]
252 
253 where @f$ s = W_{till} / W_{till}^{max} @f$, @f$ W_{till}^{max} @f$ =`hydrology_tillwat_max`,
254 @f$ \delta @f$ =`basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden`, @f$ P_o @f$ is the
255 overburden pressure, @f$ N_0 @f$ =`basal_yield_stress.mohr_coulomb.till_reference_effective_pressure` is a
256 reference effective pressure, @f$ e_0 @f$ =`basal_yield_stress.mohr_coulomb.till_reference_void_ratio` is the void ratio
257 at the reference effective pressure, and @f$ C_c @f$ =`basal_yield_stress.mohr_coulomb.till_compressibility_coefficient`
258 is the coefficient of compressibility of the till. Constants @f$ N_0, e_0, C_c @f$ are
259 found by [@ref Tulaczyketal2000] from laboratory experiments on samples of
260 till.
261 
262 If `basal_yield_stress.add_transportable_water` is yes then @f$ s @f$ in the above formula
263 becomes @f$ s = (W + W_{till}) / W_{till}^{max} @f$,
264 that is, the water amount is the sum @f$ W+W_{till} @f$.
265  */
267  double t, double dt) {
268  (void) t;
269  (void) dt;
270 
271  bool slippery_grounding_lines = m_config->get_flag("basal_yield_stress.slippery_grounding_lines"),
272  add_transportable_water = m_config->get_flag("basal_yield_stress.add_transportable_water");
273 
274  const double
275  ice_density = m_config->get_number("constants.ice.density"),
276  standard_gravity = m_config->get_number("constants.standard_gravity");
277 
278  const double high_tauc = m_config->get_number("basal_yield_stress.ice_free_bedrock"),
279  W_till_max = m_config->get_number("hydrology.tillwat_max"),
280  delta = m_config->get_number("basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden"),
281  tlftw = m_config->get_number("basal_yield_stress.mohr_coulomb.till_log_factor_transportable_water");
282 
284 
285  const array::Scalar
286  &W_till = *inputs.till_water_thickness,
287  &W_subglacial = *inputs.subglacial_water_thickness,
288  &ice_thickness = inputs.geometry->ice_thickness;
289 
290  const auto &cell_type = inputs.geometry->cell_type;
291  const auto &bed_topography = inputs.geometry->bed_elevation;
292  const auto &sea_level = inputs.geometry->sea_level_elevation;
293 
294  array::AccessScope list{&W_till, &m_till_phi, &m_basal_yield_stress, &cell_type,
295  &bed_topography, &sea_level, &ice_thickness};
296 
297  if (add_transportable_water) {
298  list.add(W_subglacial);
299  }
300 
301  if (m_delta) {
302  m_delta->update(t, dt);
303  m_delta->average(t, dt);
304  list.add(*m_delta);
305  }
306 
307  for (auto p = m_grid->points(); p; p.next()) {
308  const int i = p.i(), j = p.j();
309 
310  if (cell_type.ice_free(i, j)) {
311  m_basal_yield_stress(i, j) = high_tauc; // large yield stress if ice-free
312  } else { // grounded and there is some ice
313 
314  // user can ask that marine grounding lines get special treatment
315  double water = W_till(i,j); // usual case
316 
317  if (slippery_grounding_lines and
318  bed_topography(i, j) <= sea_level(i, j) and
319  (cell_type.next_to_floating_ice(i, j) or cell_type.next_to_ice_free_ocean(i, j))) {
320  water = W_till_max;
321  } else if (add_transportable_water) {
322  water = W_till(i, j) + tlftw * log(1.0 + W_subglacial(i, j) / tlftw);
323  }
324 
325  double P_overburden = ice_density * standard_gravity * ice_thickness(i, j);
326 
327  m_basal_yield_stress(i, j) = mc.yield_stress(m_delta ? (*m_delta)(i, j) : delta,
328  P_overburden, water, m_till_phi(i, j));
329  }
330  }
331 
333 }
334 
335 //! Computes the till friction angle phi as a piecewise linear function of bed elevation, according to user options.
336 /*!
337 Computes the till friction angle \f$\phi(x,y)\f$ at a location as the following
338 increasing, piecewise-linear function of the bed elevation \f$b(x,y)\f$. Let
339  \f[ M = (\phi_{\text{max}} - \phi_{\text{min}}) / (b_{\text{max}} - b_{\text{min}}) \f]
340 be the slope of the nontrivial part. Then
341  \f[ \phi(x,y) = \begin{cases}
342  \phi_{\text{min}}, & b(x,y) \le b_{\text{min}}, \\
343  \phi_{\text{min}} + (b(x,y) - b_{\text{min}}) \,M,
344  & b_{\text{min}} < b(x,y) < b_{\text{max}}, \\
345  \phi_{\text{max}}, & b_{\text{max}} \le b(x,y), \end{cases} \f]
346 where \f$\phi_{\text{min}}=\f$`phi_min`, \f$\phi_{\text{max}}=\f$`phi_max`,
347 \f$b_{\text{min}}=\f$`topg_min`, \f$b_{\text{max}}=\f$`topg_max`.
348 
349 The default values are vaguely suitable for Antarctica. See src/pism_config.cdl.
350 */
352  array::Scalar &result) {
353 
354  const double
355  phi_min = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.phi_min"),
356  phi_max = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.phi_max"),
357  topg_min = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.topg_min"),
358  topg_max = m_config->get_number("basal_yield_stress.mohr_coulomb.topg_to_phi.topg_max");
359 
360  m_log->message(2,
361  " till friction angle (phi) is piecewise-linear function of bed elev (topg):\n"
362  " / %5.2f for topg < %.f\n"
363  " phi = | %5.2f + (topg - (%.f)) * (%.2f / %.f) for %.f < topg < %.f\n"
364  " \\ %5.2f for %.f < topg\n",
365  phi_min, topg_min,
366  phi_min, topg_min, phi_max - phi_min, topg_max - topg_min, topg_min, topg_max,
367  phi_max, topg_max);
368 
369  if (phi_min >= phi_max) {
371  "invalid -topg_to_phi arguments: phi_min < phi_max is required");
372  }
373 
374  if (topg_min >= topg_max) {
376  "invalid -topg_to_phi arguments: topg_min < topg_max is required");
377  }
378 
379  const double slope = (phi_max - phi_min) / (topg_max - topg_min);
380 
381  array::AccessScope list{&bed_topography, &result};
382 
383  for (auto p = m_grid->points(); p; p.next()) {
384  const int i = p.i(), j = p.j();
385  const double bed = bed_topography(i, j);
386 
387  if (bed <= topg_min) {
388  result(i, j) = phi_min;
389  } else if (bed >= topg_max) {
390  result(i, j) = phi_max;
391  } else {
392  result(i, j) = phi_min + (bed - topg_min) * slope;
393  }
394  }
395 
396  // communicate ghosts so that the tauc computation can be performed locally
397  // (including ghosts of tauc, that is)
398  result.update_ghosts();
399 }
400 
401 /*!
402  * Compute the till friction angle in grounded areas using available basal yield stress,
403  * till water thickness, and overburden pressure.
404  *
405  * This is the inverse of the formula used by `update_impl()`.
406  */
408  const array::Scalar &till_water_thickness,
409  const array::Scalar &ice_thickness,
410  const array::CellType &cell_type,
411  array::Scalar &result) {
412 
414 
415  double
416  ice_density = m_config->get_number("constants.ice.density"),
417  standard_gravity = m_config->get_number("constants.standard_gravity");
418 
419  double
420  delta = m_config->get_number("basal_yield_stress.mohr_coulomb.till_effective_fraction_overburden");
421 
422  const array::Scalar
423  &W_till = till_water_thickness;
424 
425  array::AccessScope list{&cell_type, &basal_yield_stress, &W_till, &ice_thickness, &result};
426 
427  for (auto p = m_grid->points(); p; p.next()) {
428  const int i = p.i(), j = p.j();
429 
430  if (cell_type.ocean(i, j)) {
431  // no change
432  } else if (cell_type.ice_free(i, j)) {
433  // no change
434  } else { // grounded and there is some ice
435  double P_overburden = ice_density * standard_gravity * ice_thickness(i, j);
436 
437  result(i, j) = mc.till_friction_angle(delta, P_overburden, W_till(i, j), basal_yield_stress(i, j));
438  }
439  }
440 
441  result.update_ghosts();
442 }
443 
445  return combine({{"tillphi", Diagnostic::wrap(m_till_phi)}},
447 }
448 
449 } // end of namespace pism
const Time & time() const
Definition: Component.cc:109
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
void regrid(const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:159
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:160
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:361
std::string filename() const
Definition: File.cc:307
High-level PISM I/O class.
Definition: File.hh:56
array::Scalar1 sea_level_elevation
Definition: Geometry.hh:48
array::CellType2 cell_type
Definition: Geometry.hh:55
array::Scalar2 ice_thickness
Definition: Geometry.hh:51
array::Scalar2 bed_elevation
Definition: Geometry.hh:47
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
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
MohrCoulombYieldStress(std::shared_ptr< const Grid > g)
void set_till_friction_angle(const array::Scalar &input)
void write_model_state_impl(const File &output) const
The default (empty implementation).
void define_model_state_impl(const File &output) const
void finish_initialization(const YieldStressInputs &inputs)
MaxTimestep max_timestep_impl(double t) const
std::shared_ptr< array::Forcing > m_delta
void till_friction_angle(const array::Scalar &bed_topography, array::Scalar &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.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
std::string get_name() const
VariableMetadata & set_time_independent(bool flag)
const array::Scalar * till_water_thickness
Definition: YieldStress.hh:34
const array::Scalar * subglacial_water_thickness
Definition: YieldStress.hh:36
const Geometry * geometry
Definition: YieldStress.hh:32
std::string m_name
Definition: YieldStress.hh:76
void update(const YieldStressInputs &inputs, double t, double dt)
Definition: YieldStress.cc:75
std::string name() const
Definition: YieldStress.cc:108
DiagnosticList diagnostics_impl() const
Definition: YieldStress.cc:104
array::Scalar2 m_basal_yield_stress
Definition: YieldStress.hh:74
The PISM basal yield stress model interface (virtual base class)
Definition: YieldStress.hh:43
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
void add(double alpha, const Array2D< T > &x)
Definition: Array2D.hh:65
void read(const std::string &filename, unsigned int time)
Definition: Array.cc:809
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition: Array.cc:540
void write(const std::string &filename) const
Definition: Array.cc:800
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void regrid(const std::string &filename, io::Default default_value)
Definition: Array.cc:814
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
bool ice_free(int i, int j) const
Definition: CellType.hh:54
bool ocean(int i, int j) const
Definition: CellType.hh:34
"Cell type" mask. Adds convenience methods to array::Scalar.
Definition: CellType.hh:30
static Default Nil()
Definition: IO_Flags.hh:97
#define PISM_ERROR_LOCATION
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
T combine(const T &a, const T &b)
std::string filename
Definition: options.hh:33