PISM, A Parallel Ice Sheet Model  stable v2.0.6 committed by Constantine Khrulev on 2023-01-23 15:14:38 -0900
GoldsbyKohlstedt.cc
Go to the documentation of this file.
1 /* Copyright (C) 2015, 2016, 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  */
19 
20 #include <cmath>
21 #include <stdexcept>
22 #include <gsl/gsl_math.h> // M_PI
23 
24 #include "GoldsbyKohlstedt.hh"
25 
26 namespace pism {
27 namespace rheology {
28 
29 // Goldsby-Kohlstedt (forward) ice flow law
30 
31 GoldsbyKohlstedt::GoldsbyKohlstedt(const std::string &prefix,
32  const Config &config, EnthalpyConverter::Ptr ec)
33  : FlowLaw(prefix, config, ec) {
34  m_name = "Goldsby-Kohlstedt / Paterson-Budd (hybrid)";
35 
36  m_V_act_vol = -13.e-6; // m^3/mol
37  m_d_grain_size = 1.0e-3; // m (see p. ?? of G&K paper)
38  //--- dislocation creep ---
39  m_disl_crit_temp = 258.0; // Kelvin
40  //disl_A_cold = 4.0e5; // MPa^{-4.0} s^{-1}
41  //disl_A_warm = 6.0e28; // MPa^{-4.0} s^{-1}
42  m_disl_A_cold = 4.0e-19; // Pa^{-4.0} s^{-1}
43  m_disl_A_warm = 6.0e4; // Pa^{-4.0} s^{-1} (GK)
44  m_disl_n = 4.0; // stress exponent
45  m_disl_Q_cold = 60.e3; // J/mol Activation energy
46  m_disl_Q_warm = 180.e3; // J/mol Activation energy (GK)
47  //--- grain boundary sliding ---
48  m_gbs_crit_temp = 255.0; // Kelvin
49  //gbs_A_cold = 3.9e-3; // MPa^{-1.8} m^{1.4} s^{-1}
50  //gbs_A_warm = 3.e26; // MPa^{-1.8} m^{1.4} s^{-1}
51  m_gbs_A_cold = 6.1811e-14; // Pa^{-1.8} m^{1.4} s^{-1}
52  m_gbs_A_warm = 4.7547e15; // Pa^{-1.8} m^{1.4} s^{-1}
53  m_gbs_n = 1.8; // stress exponent
54  m_gbs_Q_cold = 49.e3; // J/mol Activation energy
55  m_gbs_Q_warm = 192.e3; // J/mol Activation energy
56  m_p_grain_sz_exp = 1.4; // from Peltier
57  //--- easy slip (basal) ---
58  //basal_A = 5.5e7; // MPa^{-2.4} s^{-1}
59  m_basal_A = 2.1896e-7; // Pa^{-2.4} s^{-1}
60  m_basal_n = 2.4; // stress exponent
61  m_basal_Q = 60.e3; // J/mol Activation energy
62  //--- diffusional flow ---
63  m_diff_crit_temp = 258.0; // when to use enhancement factor
64  m_diff_V_m = 1.97e-5; // Molar volume (m^3/mol)
65  m_diff_D_0v = 9.10e-4; // Preexponential volume diffusion (m^2 second-1)
66  m_diff_Q_v = 59.4e3; // activation energy, vol. diff. (J/mol)
67  m_diff_D_0b = 5.8e-4; // preexponential grain boundary coeff.
68  m_diff_Q_b = 49.e3; // activation energy, g.b. (J/mol)
69  m_diff_delta = 9.04e-10; // grain boundary width (m)
70 }
71 
72 double GoldsbyKohlstedt::flow_impl(double stress, double E,
73  double pressure, double grainsize) const {
74  double temp = m_EC->temperature(E, pressure);
75  return flow_from_temp(stress, temp, pressure, grainsize);
76 }
77 
78 double GoldsbyKohlstedt::hardness_impl(double enthalpy, double pressure) const {
79 
80  // We use the Paterson-Budd relation for the hardness parameter. It would be nice if we didn't
81  // have to, but we currently need ice hardness to compute the strain heating. See
82  // SIAFD::compute_volumetric_strain_heating().
83  double
84  T_pa = m_EC->pressure_adjusted_temperature(enthalpy, pressure),
85  A = softness_paterson_budd(T_pa);
86 
87  return pow(A, m_hardness_power);
88 }
89 
90 double GoldsbyKohlstedt::softness_impl(double , double) const {
91  throw std::runtime_error("double GoldsbyKohlstedt::softness is not implemented");
92 
93 #ifndef __GNUC__
94  return 0;
95 #endif
96 }
97 
98 /*!
99  This is the (forward) Goldsby-Kohlstedt flow law. See:
100  D. L. Goldsby & D. L. Kohlstedt (2001), "Superplastic deformation
101  of ice: experimental observations", J. Geophys. Res. 106(M6), 11017-11030.
102 */
103 double GoldsbyKohlstedt::flow_from_temp(double stress, double temp,
104  double pressure, double gs) const {
105  double eps_diff, eps_disl, eps_basal, eps_gbs, diff_D_b;
106 
107  if (fabs(stress) < 1e-10) {
108  return 0;
109  }
110  const double T = temp + (m_beta_CC_grad / (m_rho * m_standard_gravity)) * pressure;
111  const double pV = pressure * m_V_act_vol;
112  const double RT = m_ideal_gas_constant * T;
113  // Diffusional Flow
114  const double diff_D_v = m_diff_D_0v * exp(-m_diff_Q_v/RT);
115  diff_D_b = m_diff_D_0b * exp(-m_diff_Q_b/RT);
116  if (T > m_diff_crit_temp) {
117  diff_D_b *= 1000; // Coble creep scaling
118  }
119  eps_diff = 42 * m_diff_V_m *
120  (diff_D_v + M_PI * m_diff_delta * diff_D_b / gs) / (RT*(gs*gs));
121  // Dislocation Creep
122  if (T > m_disl_crit_temp) {
123  eps_disl = m_disl_A_warm * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_warm + pV)/RT);
124  } else {
125  eps_disl = m_disl_A_cold * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_cold + pV)/RT);
126  }
127  // Basal Slip
128  eps_basal = m_basal_A * pow(stress, m_basal_n-1) * exp(-(m_basal_Q + pV)/RT);
129  // Grain Boundary Sliding
130  if (T > m_gbs_crit_temp) {
131  eps_gbs = m_gbs_A_warm * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
132  exp(-(m_gbs_Q_warm + pV)/RT);
133  } else {
134  eps_gbs = m_gbs_A_cold * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
135  exp(-(m_gbs_Q_cold + pV)/RT);
136  }
137 
138  return eps_diff + eps_disl + (eps_basal * eps_gbs) / (eps_basal + eps_gbs);
139 }
140 
141 
142 /*****************
143 THE NEXT PROCEDURE REPEATS CODE; INTENDED ONLY FOR DEBUGGING
144 *****************/
145 GKparts GoldsbyKohlstedt::flowParts(double stress, double temp, double pressure) const {
146  double gs, eps_diff, eps_disl, eps_basal, eps_gbs, diff_D_b;
147  GKparts p;
148 
149  if (fabs(stress) < 1e-10) {
150  p.eps_total=0.0;
151  p.eps_diff=0.0; p.eps_disl=0.0; p.eps_gbs=0.0; p.eps_basal=0.0;
152  return p;
153  }
154  const double T = temp + (m_beta_CC_grad / (m_rho * m_standard_gravity)) * pressure;
155  const double pV = pressure * m_V_act_vol;
156  const double RT = m_ideal_gas_constant * T;
157  // Diffusional Flow
158  const double diff_D_v = m_diff_D_0v * exp(-m_diff_Q_v/RT);
159  diff_D_b = m_diff_D_0b * exp(-m_diff_Q_b/RT);
160  if (T > m_diff_crit_temp) {
161  diff_D_b *= 1000; // Coble creep scaling
162  }
163  gs = m_d_grain_size;
164  eps_diff = 14 * m_diff_V_m *
165  (diff_D_v + M_PI * m_diff_delta * diff_D_b / gs) / (RT*(gs*gs));
166  // Dislocation Creep
167  if (T > m_disl_crit_temp) {
168  eps_disl = m_disl_A_warm * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_warm + pV)/RT);
169  } else {
170  eps_disl = m_disl_A_cold * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_cold + pV)/RT);
171  }
172  // Basal Slip
173  eps_basal = m_basal_A * pow(stress, m_basal_n-1) * exp(-(m_basal_Q + pV)/RT);
174  // Grain Boundary Sliding
175  if (T > m_gbs_crit_temp) {
176  eps_gbs = m_gbs_A_warm * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
177  exp(-(m_gbs_Q_warm + pV)/RT);
178  } else {
179  eps_gbs = m_gbs_A_cold * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
180  exp(-(m_gbs_Q_cold + pV)/RT);
181  }
182 
183  p.eps_diff=eps_diff;
184  p.eps_disl=eps_disl;
185  p.eps_basal=eps_basal;
186  p.eps_gbs=eps_gbs;
187  p.eps_total=eps_diff + eps_disl + (eps_basal * eps_gbs) / (eps_basal + eps_gbs);
188  return p;
189 }
190 /*****************/
191 
193  const Config &config,
195  : GoldsbyKohlstedt(prefix, config, ec) {
196  m_name = "Goldsby-Kohlstedt / Paterson-Budd (hybrid, simplified)";
197 
198  m_d_grain_size_stripped = 3.0e-3; // m; = 3mm (see Peltier et al 2000 paper)
199 }
200 
201 
202 double GoldsbyKohlstedtStripped::flow_from_temp(double stress, double temp, double pressure, double) const {
203  // note value of gs is ignored
204  // note pressure only effects the temperature; the "P V" term is dropped
205  // note no diffusional flow
206  double eps_disl, eps_basal, eps_gbs;
207 
208  if (fabs(stress) < 1e-10) {
209  return 0;
210  }
211  const double T = temp + (m_beta_CC_grad / (m_rho * m_standard_gravity)) * pressure;
212  const double RT = m_ideal_gas_constant * T;
213  // NO Diffusional Flow
214  // Dislocation Creep
215  if (T > m_disl_crit_temp) {
216  eps_disl = m_disl_A_warm * pow(stress, m_disl_n-1) * exp(-m_disl_Q_warm/RT);
217  } else {
218  eps_disl = m_disl_A_cold * pow(stress, m_disl_n-1) * exp(-m_disl_Q_cold/RT);
219  }
220  // Basal Slip
221  eps_basal = m_basal_A * pow(stress, m_basal_n-1) * exp(-m_basal_Q/RT);
222  // Grain Boundary Sliding
223  if (T > m_gbs_crit_temp) {
224  eps_gbs = m_gbs_A_warm *
225  (pow(stress, m_gbs_n-1) / pow(m_d_grain_size_stripped, m_p_grain_sz_exp)) *
226  exp(-m_gbs_Q_warm/RT);
227  } else {
228  eps_gbs = m_gbs_A_cold *
229  (pow(stress, m_gbs_n-1) / pow(m_d_grain_size_stripped, m_p_grain_sz_exp)) *
230  exp(-m_gbs_Q_cold/RT);
231  }
232 
233  return eps_disl + (eps_basal * eps_gbs) / (eps_basal + eps_gbs);
234 }
235 
236 
237 } // end of namespace rheology
238 } // end of namespace pism
A class for storing and accessing PISM configuration flags and parameters.
std::shared_ptr< EnthalpyConverter > Ptr
double softness_paterson_budd(double T_pa) const
Return the softness parameter A(T) for a given temperature T.
Definition: FlowLaw.cc:81
double m_standard_gravity
acceleration due to gravity
Definition: FlowLaw.hh:149
std::string m_name
Definition: FlowLaw.hh:116
double m_rho
ice density
Definition: FlowLaw.hh:119
double m_hardness_power
; used to compute hardness
Definition: FlowLaw.hh:135
double m_beta_CC_grad
Clausius-Clapeyron gradient.
Definition: FlowLaw.hh:121
double m_ideal_gas_constant
ideal gas constant
Definition: FlowLaw.hh:151
EnthalpyConverter::Ptr m_EC
Definition: FlowLaw.hh:125
virtual double flow_from_temp(double stress, double temp, double pressure, double gs) const
GoldsbyKohlstedtStripped(const std::string &prefix, const Config &config, EnthalpyConverter::Ptr EC)
virtual double flow_impl(double stress, double E, double pressure, double grainsize) const
double hardness_impl(double E, double p) const
double softness_impl(double E, double p) const __attribute__((noreturn))
GKparts flowParts(double stress, double temp, double pressure) const
virtual double flow_from_temp(double stress, double temp, double pressure, double gs) const
GoldsbyKohlstedt(const std::string &prefix, const Config &config, EnthalpyConverter::Ptr EC)
A hybrid of Goldsby-Kohlstedt (2001) ice (constitutive form) and Paterson-Budd (1982)-Glen (viscosity...