Index: /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 28232)
+++ /issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp	(revision 28233)
@@ -191,4 +191,6 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.omega",HydrologyOmegaEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.istransition",HydrologyIsTransitionEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.isincludesheetthickness",HydrologyIsIncludeSheetThicknessEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.creep_open_flag",HydrologyCreepOpenFlagEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.englacial_void_ratio",HydrologyEnglacialVoidRatioEnum));
 
@@ -241,5 +243,9 @@
 	/*Get all inputs and parameters*/
 	bool istransition;
+	bool isincludesheetthickness;
+	bool creep_open_flag;
 	element->FindParam(&istransition,HydrologyIsTransitionEnum);
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+	element->FindParam(&creep_open_flag,HydrologyCreepOpenFlagEnum);
 	IssmDouble alpha     = element->FindParam(HydrologySheetAlphaEnum);
 	IssmDouble beta      = element->FindParam(HydrologySheetBetaEnum);
@@ -305,7 +311,13 @@
 
 		/*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
-		phi_0 = rho_water*g*b + rho_ice*g*H;
+		phi_0   = rho_water*g*b + rho_ice*g*H;
+		if(isincludesheetthickness) phi_0 += rho_water*g*h;
 		A=pow(B,-n);
 		v1 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*( - n));
+		if (!creep_open_flag) {
+			if (phi_0-phi<0) {
+				v1 = 0;
+			}
+		}
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
@@ -354,4 +366,8 @@
 
 	/*Retrieve all inputs and parameters*/
+	bool isincludesheetthickness;
+	bool creep_open_flag;
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+	element->FindParam(&creep_open_flag,HydrologyCreepOpenFlagEnum);
 	element->GetVerticesCoordinates(&xyz_list);
 	element->FindParam(&meltflag,HydrologyMeltFlagEnum);
@@ -422,7 +438,13 @@
 
 		/*Compute closing rate*/
-		phi_0 = rho_water*g*b + rho_ice*g*H;
+		phi_0   = rho_water*g*b + rho_ice*g*H;
+		if(isincludesheetthickness) phi_0 += rho_water*g*h;
 		A=pow(B,-n);
 		v2 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
+		if (!creep_open_flag) {
+			if (phi_0-phi<0) {
+				v2 = 0.;
+			}
+		}
 
 		for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
@@ -645,4 +667,8 @@
 
 	/*Retrieve all inputs and parameters*/
+	bool isincludesheetthickness;
+	bool creep_open_flag;
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+	element->FindParam(&creep_open_flag,HydrologyCreepOpenFlagEnum);
 	IssmDouble  dt       = element->FindParam(TimesteppingTimeStepEnum);
 	IssmDouble  l_r      = element->FindParam(HydrologyCavitySpacingEnum);
@@ -687,5 +713,6 @@
 
 		/*Get values for a few potentials*/
-		phi_0 = rho_water*g*b + rho_ice*g*H;
+		phi_0   = rho_water*g*b + rho_ice*g*H;
+		if(isincludesheetthickness) phi_0 += rho_water*g*h_old;
 		N = phi_0 - phi;
 
@@ -699,8 +726,18 @@
 		if(h_old<h_r){
 			alpha = -ub/l_r - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+			if (!creep_open_flag) {
+				if (N<0) {
+					alpha = -ub/l_r;
+				}
+			}
 			beta  = ub*h_r/l_r;
 		}
 		else{
 			alpha = - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+			if (!creep_open_flag) {
+				if (N<0) {
+					alpha = 0.;
+				}
+			}
 			beta  = 0.;
 		}
@@ -732,5 +769,5 @@
 	/*Intermediary*/
 	IssmDouble phi_0, phi_m, p_i;
-	IssmDouble H,b,phi;
+	IssmDouble H,b,phi,h;
 	IssmDouble oceanLS,iceLS;
 
@@ -738,4 +775,7 @@
 
 	/*Get thickness and base on nodes to apply cap on water head*/
+	bool isincludesheetthickness;
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
+	Input *h_input       = element->GetInput(HydrologySheetThicknessEnum);    _assert_(h_input);
    IssmDouble* N = xNew<IssmDouble>(numnodes);
 	IssmDouble  rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -767,4 +807,5 @@
 		oceanLS_input->GetInputValue(&oceanLS,gauss);
 		iceLS_input->GetInputValue(&iceLS,gauss);
+		h_input->GetInputValue(&h,gauss);
 
 		/*Elevation potential*/
@@ -776,4 +817,6 @@
 		/*Compute overburden potential*/
 		phi_0 = phi_m + p_i;
+		phi_0   = rho_water*g*b + rho_ice*g*H;
+		if(isincludesheetthickness) phi_0 += rho_water*g*h;
 
 		/*Calculate effective pressure*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 28232)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 28233)
@@ -380,4 +380,6 @@
 	bool istransition;
 	element->FindParam(&istransition,HydrologyIsTransitionEnum);
+	bool isincludesheetthickness;
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
 	IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
 	IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
@@ -435,4 +437,5 @@
 		/*Get values for a few potentials*/
 		phi_0   = rho_water*g*b + rho_ice*g*H;
+		if(isincludesheetthickness) phi_0 += rho_water*g*h;
 		dphids  = dphi[0]*tx + dphi[1]*ty;
 		dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
@@ -541,4 +544,6 @@
 	bool istransition;
 	element->FindParam(&istransition,HydrologyIsTransitionEnum);
+	bool isincludesheetthickness;
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
 	IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
 	IssmDouble mu_water  = element->FindParam(MaterialsMuWaterEnum);
@@ -592,4 +597,5 @@
 		/*Get values for a few potentials*/
 		phi_0   = rho_water*g*b + rho_ice*g*H;
+		if(isincludesheetthickness) phi_0 += rho_water*g*h;
 		dphids  = dphi[0]*tx + dphi[1]*ty;
 		dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
@@ -691,4 +697,6 @@
 	bool istransition;
 	element->FindParam(&istransition,HydrologyIsTransitionEnum);
+	bool isincludesheetthickness;
+	element->FindParam(&isincludesheetthickness,HydrologyIsIncludeSheetThicknessEnum);
 	IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
 	IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -738,4 +746,5 @@
 	/*Get values for a few potentials*/
 	phi_0   = rho_water*g*b + rho_ice*g*H;
+	if(isincludesheetthickness) phi_0 += rho_water*g*h;
 	dphids  = dphi[0]*tx + dphi[1]*ty;
 	dphimds = rho_water*g*(db[0]*tx + db[1]*ty);