Changeset 23143

Timestamp:

08/21/18 07:02:50 (7 years ago)

Author:

hakesson

Message:

CHG: working on PISM friction law

File:

• issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -676 +676 @@
 }/*}}}*/
 void Friction::GetAlpha2PISM(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
-        /*Here, we want to parameterize the friction as
-         *
-         * ... more comments to come
+        /*Here, we want to parameterize the friction using a pseudoplastic friction law,
+         * computing the basal shear stress as
+         *
+         * alpha2 = tau_c (u_b/(abs(u_b)^(1-q)*u_0^q))
+         *
+         * The yield stress tau_c is a function of the effective pressure N
+         * using a Mohr-Coloumb criterion, so that
+         * tau_c = tan(phi)*N,
+         * where phi is the till friction angle, representing sediment strength
+         *
+         * The effective pressure is given by Eq. (5) in Aschwanden et al. 2016:
+         *
+         * N = delta * P0 * 10^((e_0/Cc)(1-(W/Wmax)))
+         *
+         * W is calculated by a non-conserving hydrology model in HydrologyPismAnalysis.cpp
+         *
+         * see Aschwanden et al. 2016 and Bueler and Brown, 2009 for more details
          */
 
-        /*compute ice pressure P0*/
+        /*compute ice overburden pressure P0*/
         IssmDouble thickness,base,P0;
         element->GetInputValue(&thickness, gauss,ThicknessEnum);
@@ -697 +711 @@
         element->GetInputValue(&W,gauss,WatercolumnEnum);
         element->GetInputValue(&Wmax,gauss,FrictionWatercolumnMaxEnum);
+
+//      /*Check that water column height is within 0 and upper bound, correct if needed*/
+//      // if watercolumn height is higher than the maximum allowed height, set height to upper bound
+//      if(W>Wmax){
+//              W=Wmax;
+//      }
+//      // if watercolumn height is negative (shouldn't happen), set it to 0
+//      if(W<0){
+//              W=0;
+//      }
+//      // if watercolumn height is within 0 and upper bound, nothing to be done
+//      else{
+//              //do nothing
+//      }
+
         N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));
 
+        /*Get till friction angles, defined by user [deg]*/
+        IssmDouble phi,pi;
+        element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);
+
+        /*Convert till friction angle from user-defined deg to rad, which Matlab uses*/
+        phi = phi*pi/180.;
+
         /*Compute yield stress following a Mohr-Colomb criterion*/
-        IssmDouble phi;
-        element->GetInputValue(&phi,gauss,FrictionTillFrictionAngleEnum);
         IssmDouble tau_c = N*tan(phi);