Merge branch 'fracture_density_1.0' into pism_pik_1.0

solved conflicts:
	src/icemodel/diagnostics.cc

src/icemodel/diagnostics.cc

@@ -2326,8 +2326,18 @@ void IceModel::init_diagnostics() {
 
     // earth
     {"beddef_load_thickness", d::wrap(m_beddef_load)},
+
   };
 
+  // fractures
+  if (m_config->get_boolean("fracture_density.enabled")) {
+    m_diagnostics["fracture_density"] = d::wrap(m_fracture->density);
+    m_diagnostics["fracture_growth_rate"] = d::wrap(m_fracture->growth_rate);
+    m_diagnostics["fracture_healing_rate"] = d::wrap(m_fracture->healing_rate);
+    m_diagnostics["fracture_flow_enhancement"] = d::wrap(m_fracture->flow_enhancement);
+    m_diagnostics["fracture_toughness"] = d::wrap(m_fracture->toughness);
+  }
+
 #if (PISM_USE_PROJ4==1)
   std::string proj4 = m_grid->get_mapping_info().proj4;
   if (not proj4.empty()) {

src/icemodel/fracture_density.cc

@@ -29,9 +29,13 @@
 #include "pism/util/pism_options.hh"
 #include "IceModel.hh"
 
+#include "pism/rheology/FlowLaw.hh"
+#include "pism/stressbalance/ShallowStressBalance.hh"
+#include "pism/energy/EnergyModel.hh"
+
 namespace pism {
 
-//! \file iMfractures.cc implementing calculation of fracture density with PIK options -fractures.
+//! \file fracture_density.cc implementing calculation of fracture density with PIK options -fractures.
 
 void IceModel::update_fracture_density() {
   const double dx = m_grid->dx(), dy = m_grid->dy(), Mx = m_grid->Mx(), My = m_grid->My();
@@ -63,6 +67,19 @@ void IceModel::update_fracture_density() {
     A_new.copy_from(A);
   }
 
+  double glenexp = m_config->get_double("stress_balance.ssa.Glen_exponent");
+  double softness = m_config->get_double("flow_law.isothermal_Glen.ice_softness"); //const
+  double minH = m_config->get_double("stress_balance.ice_free_thickness_standard");
+
+  const bool borstad_limit = m_config->get_boolean("fracture_density.borstad_limit");
+  //if (borstad_limit){
+    const IceModelVec3 &enthalpy = m_energy_model->enthalpy();
+    list.add(enthalpy);
+    const double *z = &m_grid->z()[0];
+    const rheology::FlowLaw*
+    flow_law = m_stress_balance->shallow()->flow_law();
+  //}
+
   //options
   /////////////////////////////////////////////////////////
   double soft_residual = options::Real("-fracture_softening", "soft_residual", 1.0);
@@ -221,12 +238,60 @@ void IceModel::update_fracture_density() {
       sigmat = KSImax;
     }
 
-    //////////////////////////////////////////////////////////////////////////////
 
+    //////////////////////////////////////////////////////////////////////////////
     //fracture density
-    double fdnew = gamma * (strain_rates(i, j, 0) - 0.0) * (1 - D_new(i, j));
-    if (sigmat > initThreshold) {
-      D_new(i, j) += fdnew * m_dt;
+    double fdnew = 0.0;
+    bool add_fractures;
+
+    //Borstad et al. 2016, constitutive framework for ice weakening
+    if (borstad_limit) {
+
+      double H = m_geometry.ice_thickness(i, j);
+      if (H > minH) {
+
+      //get vertical average hardness
+      unsigned int k = m_grid->kBelowHeight(H);
+      double hardness = averaged_hardness(*flow_law, H, k, &z[0], enthalpy.get_column(i, j));
+      softness = pow(hardness, -glenexp);
+
+      //mean parameters from paper
+      double t0 = 130000.0; //Pa
+      t0 = initThreshold;
+      double kappa = 2.8;
+
+      //effective strain rate
+      double e1=strain_rates(i, j, 0);
+      double e2=strain_rates(i, j, 1);
+      double ee = sqrt(PetscSqr(e1) + PetscSqr(e2) - e1 * e2);
+
+      //threshold for unfractured ice
+      double e0 = pow( (t0/hardness) , glenexp); 
+
+      //threshold for fractured ice (exponential law)
+      double ex = exp((e0-ee)/(e0*(kappa-1)));
+
+      // stress threshold for fractures ice
+      double te = t0 * ex; 
+      sigmat    = te;
+
+      //actual effective stress
+      double ts = hardness * pow(ee,1/glenexp) * (1-D_new(i, j));
+
+      //fracture formation if threshold is hit
+      add_fractures = (ts > te && ee > e0);
+      if (add_fractures) {
+        fdnew = 1.0- ( ex * pow((ee/e0),-1/glenexp) ); //new fracture density
+        D_new(i, j) = fdnew;
+      }
+      }
+
+    } else { //default fracture growth
+      fdnew = gamma * (strain_rates(i, j, 0) - 0.0) * (1 - D_new(i, j));
+      add_fractures = (sigmat > initThreshold);
+      if (add_fractures) {
+        D_new(i, j) += fdnew * m_dt;
+      }
     }
 
     //healing
@@ -261,16 +326,15 @@ void IceModel::update_fracture_density() {
     // write related fracture quantities to nc-file
     // if option -write_fd_fields is set
     if (write_fd && m_geometry.ice_thickness(i, j) > 0.0) {
+
       //fracture toughness
       m_fracture->toughness(i, j) = sigmat;
 
       // fracture growth rate
-      if (sigmat > initThreshold) {
-        m_fracture->growth_rate(i, j) = fdnew;
-        //m_fracture->growth_rate(i,j)=gamma*(vPrinStrain1(i,j)-0.0)*(1-D_new(i,j));
-      } else {
-        m_fracture->growth_rate(i, j) = 0.0;
-      }
+      if (add_fractures)
+          m_fracture->growth_rate(i, j) = fdnew;
+      else
+          m_fracture->growth_rate(i, j) = 0.0;
 
       // fracture healing rate
       if (m_geometry.ice_thickness(i, j) > 0.0) {
@@ -289,15 +353,14 @@ void IceModel::update_fracture_density() {
       A_new(i, j) -= m_dt * uvel * (uvel < 0 ? A(i + 1, j) - A(i, j) : A(i, j) - A(i - 1, j)) / dx;
       A_new(i, j) -= m_dt * vvel * (vvel < 0 ? A(i, j + 1) - A(i, j) : A(i, j) - A(i, j - 1)) / dy;
       A_new(i, j) += m_dt / one_year;
-      if (sigmat > initThreshold) {
+      if (add_fractures) {
         A_new(i, j) = 0.0;
       }
 
       // additional flow enhancement due to fracture softening
-      double phi_exp   = 3.0; //flow_law->exponent();
-      double softening = pow((1.0 - (1.0 - soft_residual) * D_new(i, j)), -phi_exp);
+      double softening = pow((1.0 - (1.0 - soft_residual) * D_new(i, j)), -glenexp);
       if (m_geometry.ice_thickness(i, j) > 0.0) {
-        m_fracture->flow_enhancement(i, j) = 1.0 / pow(softening, 1 / 3.0);
+        m_fracture->flow_enhancement(i, j) = 1.0 / pow(softening, 1 / glenexp);
       } else {
         m_fracture->flow_enhancement(i, j) = 1.0;
       }

src/pism_config.cdl

@@ -570,7 +570,7 @@ netcdf pism_config {
     pism_config:flow_law.isothermal_Glen.ice_softness_units = "Pascal-3 second-1";
 
     pism_config:fracture_density.constant_fd = "no";
-    pism_config:fracture_density.constant_fd_doc = "FIXME";
+    pism_config:fracture_density.constant_fd_doc = "Do not update fracture density field";
     pism_config:fracture_density.constant_fd_option = "constant_fd";
     pism_config:fracture_density.constant_fd_type = "boolean";
 
@@ -626,6 +626,11 @@ netcdf pism_config {
     pism_config:fracture_density.write_fields_option = "write_fd_fields";
     pism_config:fracture_density.write_fields_type = "boolean";
 
+    pism_config:fracture_density.borstad_limit = "no";
+    pism_config:fracture_density.borstad_limit_doc = "Define fracture growth dependent on empirical function from Bortstad et al., 2016";
+    pism_config:fracture_density.borstad_limit_option = "constitutive_stress_limit";
+    pism_config:fracture_density.borstad_limit_type = "boolean";
+
     pism_config:geometry.grounded_cell_fraction  = "false";
     pism_config:geometry.grounded_cell_fraction_doc = "Linear interpolation scheme ('LI' in Gladstone et al. 2010) expanded to two dimensions is used if switched on in order to evaluate the position of the grounding line on a subgrid scale.";
     pism_config:geometry.grounded_cell_fraction_option = "subgl";