tbasal_resistance.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
---
tbasal_resistance.cc (8248B)
---
1 // Copyright (C) 2004-2017, 2019 Jed Brown, Ed Bueler, and Constantine Khroulev
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 <cmath> // pow, sqrt
20
21 #include "basal_resistance.hh"
22
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/Logger.hh"
25
26 namespace pism {
27
28 static double square(double x) {
29 return x * x;
30 }
31
32 /* Purely plastic */
33
34 IceBasalResistancePlasticLaw::IceBasalResistancePlasticLaw(const Config &config) {
35 m_plastic_regularize = config.get_number("basal_resistance.plastic.regularization", "m second-1");
36 }
37
38 IceBasalResistancePlasticLaw::~IceBasalResistancePlasticLaw() {
39 // empty
40 }
41
42 void IceBasalResistancePlasticLaw::print_info(const Logger &log,
43 int threshold,
44 units::System::Ptr system) const {
45 log.message(threshold, "Using purely plastic till with eps = %10.5e m year-1.\n",
46 units::convert(system, m_plastic_regularize, "m second-1", "m year-1"));
47 }
48
49
50 //! Compute the drag coefficient for the basal shear stress.
51 double IceBasalResistancePlasticLaw::drag(double tauc, double vx, double vy) const {
52 const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
53
54 return tauc / sqrt(magreg2);
55 }
56
57 //! Compute the drag coefficient and its derivative with respect to \f$ \alpha = \frac 1 2 (u_x^2 + u_y^2) \f$
58 /**
59 * @f{align*}{
60 * \beta &= \frac{\tau_{c}}{|\mathbf{u}|} = \tau_{c}\cdot (|\mathbf{u}|^{2})^{-1/2}\\
61 * \diff{\beta}{\frac12 |\mathbf{u}|^{2}} &= -\frac{1}{|\mathbf{u}|^{2}}\cdot \beta
62 * @f}
63 */
64 void IceBasalResistancePlasticLaw::drag_with_derivative(double tauc, double vx, double vy,
65 double *beta, double *dbeta) const {
66 const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
67
68 *beta = tauc / sqrt(magreg2);
69
70 if (dbeta) {
71 *dbeta = -1 * (*beta) / magreg2;
72 }
73 }
74
75 /* Pseudo-plastic */
76
77 IceBasalResistancePseudoPlasticLaw::IceBasalResistancePseudoPlasticLaw(const Config &config)
78 : IceBasalResistancePlasticLaw(config) {
79 m_pseudo_q = config.get_number("basal_resistance.pseudo_plastic.q");
80 m_pseudo_u_threshold = config.get_number("basal_resistance.pseudo_plastic.u_threshold", "m second-1");
81 m_sliding_scale_factor_reduces_tauc = config.get_number("basal_resistance.pseudo_plastic.sliding_scale_factor");
82 }
83
84 IceBasalResistancePseudoPlasticLaw::~IceBasalResistancePseudoPlasticLaw() {
85 // empty
86 }
87
88 void IceBasalResistancePseudoPlasticLaw::print_info(const Logger &log,
89 int threshold,
90 units::System::Ptr system) const {
91
92 if (m_pseudo_q == 1.0) {
93 log.message(threshold,
94 "Using linearly viscous till with u_threshold = %.2f m year-1.\n",
95 units::convert(system, m_pseudo_u_threshold, "m second-1", "m year-1"));
96 } else {
97 log.message(threshold,
98 "Using pseudo-plastic till with eps = %10.5e m year-1, q = %.4f,"
99 " and u_threshold = %.2f m year-1.\n",
100 units::convert(system, m_plastic_regularize, "m second-1", "m year-1"),
101 m_pseudo_q,
102 units::convert(system, m_pseudo_u_threshold, "m second-1", "m year-1"));
103 }
104 }
105
106 //! Compute the drag coefficient for the basal shear stress.
107 /*!
108
109 The basal shear stress term @f$ \tau_b @f$ in the SSA stress balance
110 for ice is minus the return value here times (vx,vy). Thus this
111 method computes the basal shear stress as
112
113 @f[ \tau_b = - \frac{\tau_c}{|\mathbf{U}|^{1-q} U_{\mathtt{th}}^q} \mathbf{U} @f]
114
115 where @f$ \tau_b=(\tau_{(b)x},\tau_{(b)y}) @f$ , @f$ U=(u,v) @f$ ,
116 @f$ q= @f$ `pseudo_q`, and @f$ U_{\mathtt{th}}= @f$
117 `pseudo_u_threshold`. Typical values for the constants are
118 @f$ q=0.25 @f$ and @f$ U_{\mathtt{th}} = 100 @f$ m year-1.
119
120 The linearly-viscous till case pseudo_q = 1.0 is allowed, in which
121 case @f$ \beta = \tau_c/U_{\mathtt{th}} @f$ . The purely-plastic till
122 case pseudo_q = 0.0 is also allowed; note that there is still a
123 regularization with data member plastic_regularize.
124
125 One can scale tauc if desired:
126
127 A scale factor of @f$ A @f$ is intended to increase basal sliding rate
128 by @f$ A @f$ . It would have exactly this effect \e if the driving
129 stress were \e hypothetically completely held by the basal
130 resistance. Thus this scale factor is used to reduce (if
131 `-sliding_scale_factor_reduces_tauc` @f$ A @f$ with @f$ A > 1 @f$) or increase (if @f$ A <
132 1 @f$) the value of the (pseudo-) yield stress `tauc`. The concept
133 behind this is described at
134 [the SeaRISE wiki](http://websrv.cs.umt.edu/isis/index.php/Category_1:_Whole_Ice_Sheet#Initial_Experiment_-_E1_-_Increased_Basal_Lubrication).
135
136 Specifically, the concept behind this mechanism is to suppose
137 equality of driving and basal shear stresses,
138
139 @f[ \rho g H \nabla h = \frac{\tau_c}{|\mathbf{U}|^{1-q} U_{\mathtt{th}}^q} \mathbf{U}. @f]
140
141 (*For emphasis:* The membrane stress held by the ice itself is
142 missing from this incomplete stress balance.) Thus the pseudo yield
143 stress @f$ \tau_c @f$ would be related to the sliding speed
144 @f$ |\mathbf{U}| @f$ by
145
146 @f[ |\mathbf{U}| = \frac{C}{\tau_c^{1/q}} \f]
147
148 for some (geometry-dependent) constant @f$ C @f$ . Multiplying
149 @f$ |\mathbf{U}| @f$ by @f$ A @f$ in this equation corresponds to
150 dividing @f$ \tau_c @f$ by @f$ A^q @f$ .
151
152 Note that the mechanism has no effect whatsoever if @f$ q=0 @f$ , which
153 is the purely plastic case. In that case there is \e no direct
154 relation between the yield stress and the sliding velocity, and the
155 difference between the driving stress and the yield stress is
156 entirely held by the membrane stresses. (There is also no singular
157 mathematical operation as @f$ A^q = A^0 = 1 @f$ .)
158 */
159 double IceBasalResistancePseudoPlasticLaw::drag(double tauc, double vx, double vy) const {
160 const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
161
162 if (m_sliding_scale_factor_reduces_tauc > 0.0) {
163 double Aq = pow(m_sliding_scale_factor_reduces_tauc, m_pseudo_q);
164 return (tauc / Aq) * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
165 } else {
166 return tauc * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
167 }
168 }
169
170
171 //! Compute the drag coefficient and its derivative with respect to @f$ \alpha = \frac 1 2 (u_x^2 + u_y^2) @f$
172 /**
173 * @f{align*}{
174 * \beta &= \frac{\tau_{c}}{u_{\text{threshold}}}\cdot (|u|^{2})^{\frac{q-1}{2}} \\
175 * \diff{\beta}{\frac12 |\mathbf{u}|^{2}} &= \frac{\tau_{c}}{u_{\text{threshold}}}\cdot \frac{q-1}{2}\cdot (|\mathbf{u}|^{2})^{\frac{q-1}{2} - 1}\cdot 2 \\
176 * &= \frac{q-1}{|\mathbf{u}|^{2}}\cdot \beta(\mathbf{u}) \\
177 * @f}
178 */
179 void IceBasalResistancePseudoPlasticLaw::drag_with_derivative(double tauc, double vx, double vy,
180 double *beta, double *dbeta) const
181 {
182 const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
183
184 if (m_sliding_scale_factor_reduces_tauc > 0.0) {
185 double Aq = pow(m_sliding_scale_factor_reduces_tauc, m_pseudo_q);
186 *beta = (tauc / Aq) * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
187 } else {
188 *beta = tauc * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
189 }
190
191 if (dbeta) {
192 *dbeta = (m_pseudo_q - 1) * (*beta) / magreg2;
193 }
194
195 }
196
197 } // end of namespace pism