Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17614)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17615)
@@ -994,4 +994,16 @@
 /*}}}*/
 void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+	/*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
+	 * and Cuffey p 303 eq 8.18:
+	 *
+	 *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
+	 *
+	 *  or
+	 *
+	 *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
+	 *
+	 *          = 1/sqrt(2) ||eps||_F
+	 *
+	 *  where ||.||_F is the Frobenius norm */
 
 	/*Intermediaries*/
@@ -1008,7 +1020,7 @@
 	}
 	else{
-		/* eps_eff^2 = exx^2 + eyy^2 + 2*exy^2 */
+		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
 		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
 	}
 
@@ -1084,5 +1096,5 @@
 	IssmDouble viscosity;
 	IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
-	IssmDouble epsilon2d[2]; /* epsilon=[exx,exy];           */
+	IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
 	IssmDouble eps_eff;
 
@@ -1093,7 +1105,7 @@
 	}
 	else{
-		/* eps_eff^2 = exx^2 + exy^2 */
+		/* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/
 		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1]);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]);
 	}
 
@@ -1113,12 +1125,12 @@
 
 	 if(dim==2){
-		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/
+		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
 		 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
 		 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
 	 }
 	 else{
-		 /* eps_eff^2 = exx^2*/
+		 /* eps_eff^2 = 1/2 exx^2*/
 		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
-		 eps_eff = sqrt(epsilon1d*epsilon1d);
+		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
 	 }