pism

[fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
git clone git://src.adamsgaard.dk/pism # fast
git clone https://src.adamsgaard.dk/pism.git # slow
Log | Files | Refs | README | LICENSE Back to index

basal_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