tGoldsbyKohlstedt.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
---
tGoldsbyKohlstedt.cc (9546B)
---
1 /* Copyright (C) 2015, 2016 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::averaged_hardness_impl(double, int,
79 const double *,
80 const double *) const {
81
82 throw std::runtime_error("double GoldsbyKohlstedt::averaged_hardness is not implemented");
83
84 #ifndef __GNUC__
85 return 0;
86 #endif
87 }
88
89 double GoldsbyKohlstedt::hardness_impl(double enthalpy, double pressure) const {
90
91 // We use the Paterson-Budd relation for the hardness parameter. It would be nice if we didn't
92 // have to, but we currently need ice hardness to compute the strain heating. See
93 // SIAFD::compute_volumetric_strain_heating().
94 double
95 T_pa = m_EC->pressure_adjusted_temperature(enthalpy, pressure),
96 A = softness_paterson_budd(T_pa);
97
98 return pow(A, m_hardness_power);
99 }
100
101 double GoldsbyKohlstedt::softness_impl(double , double) const {
102 throw std::runtime_error("double GoldsbyKohlstedt::softness is not implemented");
103
104 #ifndef __GNUC__
105 return 0;
106 #endif
107 }
108
109 /*!
110 This is the (forward) Goldsby-Kohlstedt flow law. See:
111 D. L. Goldsby & D. L. Kohlstedt (2001), "Superplastic deformation
112 of ice: experimental observations", J. Geophys. Res. 106(M6), 11017-11030.
113 */
114 double GoldsbyKohlstedt::flow_from_temp(double stress, double temp,
115 double pressure, double gs) const {
116 double eps_diff, eps_disl, eps_basal, eps_gbs, diff_D_b;
117
118 if (fabs(stress) < 1e-10) {
119 return 0;
120 }
121 const double T = temp + (m_beta_CC_grad / (m_rho * m_standard_gravity)) * pressure;
122 const double pV = pressure * m_V_act_vol;
123 const double RT = m_ideal_gas_constant * T;
124 // Diffusional Flow
125 const double diff_D_v = m_diff_D_0v * exp(-m_diff_Q_v/RT);
126 diff_D_b = m_diff_D_0b * exp(-m_diff_Q_b/RT);
127 if (T > m_diff_crit_temp) {
128 diff_D_b *= 1000; // Coble creep scaling
129 }
130 eps_diff = 42 * m_diff_V_m *
131 (diff_D_v + M_PI * m_diff_delta * diff_D_b / gs) / (RT*(gs*gs));
132 // Dislocation Creep
133 if (T > m_disl_crit_temp) {
134 eps_disl = m_disl_A_warm * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_warm + pV)/RT);
135 } else {
136 eps_disl = m_disl_A_cold * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_cold + pV)/RT);
137 }
138 // Basal Slip
139 eps_basal = m_basal_A * pow(stress, m_basal_n-1) * exp(-(m_basal_Q + pV)/RT);
140 // Grain Boundary Sliding
141 if (T > m_gbs_crit_temp) {
142 eps_gbs = m_gbs_A_warm * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
143 exp(-(m_gbs_Q_warm + pV)/RT);
144 } else {
145 eps_gbs = m_gbs_A_cold * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
146 exp(-(m_gbs_Q_cold + pV)/RT);
147 }
148
149 return eps_diff + eps_disl + (eps_basal * eps_gbs) / (eps_basal + eps_gbs);
150 }
151
152
153 /*****************
154 THE NEXT PROCEDURE REPEATS CODE; INTENDED ONLY FOR DEBUGGING
155 *****************/
156 GKparts GoldsbyKohlstedt::flowParts(double stress, double temp, double pressure) const {
157 double gs, eps_diff, eps_disl, eps_basal, eps_gbs, diff_D_b;
158 GKparts p;
159
160 if (fabs(stress) < 1e-10) {
161 p.eps_total=0.0;
162 p.eps_diff=0.0; p.eps_disl=0.0; p.eps_gbs=0.0; p.eps_basal=0.0;
163 return p;
164 }
165 const double T = temp + (m_beta_CC_grad / (m_rho * m_standard_gravity)) * pressure;
166 const double pV = pressure * m_V_act_vol;
167 const double RT = m_ideal_gas_constant * T;
168 // Diffusional Flow
169 const double diff_D_v = m_diff_D_0v * exp(-m_diff_Q_v/RT);
170 diff_D_b = m_diff_D_0b * exp(-m_diff_Q_b/RT);
171 if (T > m_diff_crit_temp) {
172 diff_D_b *= 1000; // Coble creep scaling
173 }
174 gs = m_d_grain_size;
175 eps_diff = 14 * m_diff_V_m *
176 (diff_D_v + M_PI * m_diff_delta * diff_D_b / gs) / (RT*(gs*gs));
177 // Dislocation Creep
178 if (T > m_disl_crit_temp) {
179 eps_disl = m_disl_A_warm * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_warm + pV)/RT);
180 } else {
181 eps_disl = m_disl_A_cold * pow(stress, m_disl_n-1) * exp(-(m_disl_Q_cold + pV)/RT);
182 }
183 // Basal Slip
184 eps_basal = m_basal_A * pow(stress, m_basal_n-1) * exp(-(m_basal_Q + pV)/RT);
185 // Grain Boundary Sliding
186 if (T > m_gbs_crit_temp) {
187 eps_gbs = m_gbs_A_warm * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
188 exp(-(m_gbs_Q_warm + pV)/RT);
189 } else {
190 eps_gbs = m_gbs_A_cold * (pow(stress, m_gbs_n-1) / pow(gs, m_p_grain_sz_exp)) *
191 exp(-(m_gbs_Q_cold + pV)/RT);
192 }
193
194 p.eps_diff=eps_diff;
195 p.eps_disl=eps_disl;
196 p.eps_basal=eps_basal;
197 p.eps_gbs=eps_gbs;
198 p.eps_total=eps_diff + eps_disl + (eps_basal * eps_gbs) / (eps_basal + eps_gbs);
199 return p;
200 }
201 /*****************/
202
203 GoldsbyKohlstedtStripped::GoldsbyKohlstedtStripped(const std::string &prefix,
204 const Config &config,
205 EnthalpyConverter::Ptr ec)
206 : GoldsbyKohlstedt(prefix, config, ec) {
207 m_name = "Goldsby-Kohlstedt / Paterson-Budd (hybrid, simplified)";
208
209 m_d_grain_size_stripped = 3.0e-3; // m; = 3mm (see Peltier et al 2000 paper)
210 }
211
212
213 double GoldsbyKohlstedtStripped::flow_from_temp(double stress, double temp, double pressure, double) const {
214 // note value of gs is ignored
215 // note pressure only effects the temperature; the "P V" term is dropped
216 // note no diffusional flow
217 double eps_disl, eps_basal, eps_gbs;
218
219 if (fabs(stress) < 1e-10) {
220 return 0;
221 }
222 const double T = temp + (m_beta_CC_grad / (m_rho * m_standard_gravity)) * pressure;
223 const double RT = m_ideal_gas_constant * T;
224 // NO Diffusional Flow
225 // Dislocation Creep
226 if (T > m_disl_crit_temp) {
227 eps_disl = m_disl_A_warm * pow(stress, m_disl_n-1) * exp(-m_disl_Q_warm/RT);
228 } else {
229 eps_disl = m_disl_A_cold * pow(stress, m_disl_n-1) * exp(-m_disl_Q_cold/RT);
230 }
231 // Basal Slip
232 eps_basal = m_basal_A * pow(stress, m_basal_n-1) * exp(-m_basal_Q/RT);
233 // Grain Boundary Sliding
234 if (T > m_gbs_crit_temp) {
235 eps_gbs = m_gbs_A_warm *
236 (pow(stress, m_gbs_n-1) / pow(m_d_grain_size_stripped, m_p_grain_sz_exp)) *
237 exp(-m_gbs_Q_warm/RT);
238 } else {
239 eps_gbs = m_gbs_A_cold *
240 (pow(stress, m_gbs_n-1) / pow(m_d_grain_size_stripped, m_p_grain_sz_exp)) *
241 exp(-m_gbs_Q_cold/RT);
242 }
243
244 return eps_disl + (eps_basal * eps_gbs) / (eps_basal + eps_gbs);
245 }
246
247
248 } // end of namespace rheology
249 } // end of namespace pism