Index: issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp	(revision 21525)
+++ issm/trunk-jpl/src/c/analyses/HydrologySommersAnalysis.cpp	(revision 21526)
@@ -148,4 +148,5 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.friction.law",FrictionLawEnum));
    parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.relaxation",HydrologyRelaxationEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.storage",HydrologyStorageEnum));
 }/*}}}*/
 
@@ -173,4 +174,5 @@
 	ElementMatrix* Ke     = element->NewElementMatrix();
 	IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
+	IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
@@ -179,4 +181,9 @@
 	/*Get conductivity from inputs*/
 	IssmDouble conductivity = GetConductivity(element);
+
+	/*Get englacial storage coefficient*/
+	IssmDouble storage,dt;
+	element->FindParam(&storage,HydrologyStorageEnum);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
 
 	/* Start  looping on the number of gaussian points: */
@@ -187,8 +194,10 @@
 		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
 		element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+		element->NodalFunctions(basis,gauss);
 
 		for(int i=0;i<numnodes;i++){
 			for(int j=0;j<numnodes;j++){
-				Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+				Ke->values[i*numnodes+j] += conductivity*gauss->weight*Jdet*(dbasis[0*numnodes+i]*dbasis[0*numnodes+j] + dbasis[1*numnodes+i]*dbasis[1*numnodes+j])
+				  + gauss->weight*Jdet*storage/dt*basis[i]*basis[j];
 			}
 		}
@@ -244,4 +253,9 @@
 	IssmDouble conductivity = GetConductivity(element);
 
+	/*Get englacial storage coefficient*/
+	IssmDouble storage,dt;
+   element->FindParam(&storage,HydrologyStorageEnum);
+   element->FindParam(&dt,TimesteppingTimeStepEnum);
+
 	/*Build friction element, needed later: */
 	Friction* friction=new Friction(element,2);
@@ -308,4 +322,5 @@
 		  -beta*sqrt(vx*vx+vy*vy)
 		  +ieb
+		  +storage*head_old/dt
 		  )*basis[i];     	
 	}
Index: issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21525)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 21526)
@@ -176,4 +176,5 @@
    HydrologyRelaxationEnum,
 	HydrologyBasalFluxEnum,
+	HydrologyStorageEnum,
 	InversionControlParametersEnum,
 	InversionControlScalingFactorsEnum,
Index: issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21525)
+++ issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 21526)
@@ -182,4 +182,5 @@
 		case HydrologyRelaxationEnum : return "HydrologyRelaxation";
 		case HydrologyBasalFluxEnum : return "HydrologyBasalFlux";
+		case HydrologyStorageEnum : return "HydrologyStorage";
 		case InversionControlParametersEnum : return "InversionControlParameters";
 		case InversionControlScalingFactorsEnum : return "InversionControlScalingFactors";
Index: issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21525)
+++ issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 21526)
@@ -185,4 +185,5 @@
 	      else if (strcmp(name,"HydrologyRelaxation")==0) return HydrologyRelaxationEnum;
 	      else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum;
+	      else if (strcmp(name,"HydrologyStorage")==0) return HydrologyStorageEnum;
 	      else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
 	      else if (strcmp(name,"InversionControlScalingFactors")==0) return InversionControlScalingFactorsEnum;
@@ -259,9 +260,9 @@
 	      else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
 	      else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
-	      else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
          else stage=3;
    }
    if(stage==3){
-	      if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
+	      if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
+	      else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
 	      else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
 	      else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
@@ -382,9 +383,9 @@
 	      else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
 	      else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
-	      else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
+	      if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
+	      else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
 	      else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum;
 	      else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
@@ -505,9 +506,9 @@
 	      else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
 	      else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
-	      else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
          else stage=5;
    }
    if(stage==5){
-	      if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
+	      if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
+	      else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
 	      else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
 	      else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
@@ -628,9 +629,9 @@
 	      else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum;
 	      else if (strcmp(name,"Outputdefinition39")==0) return Outputdefinition39Enum;
-	      else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
          else stage=6;
    }
    if(stage==6){
-	      if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
+	      if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum;
+	      else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum;
 	      else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum;
 	      else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum;
@@ -751,9 +752,9 @@
 	      else if (strcmp(name,"Mpi")==0) return MpiEnum;
 	      else if (strcmp(name,"Mumps")==0) return MumpsEnum;
-	      else if (strcmp(name,"Gsl")==0) return GslEnum;
          else stage=7;
    }
    if(stage==7){
-	      if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
+	      if (strcmp(name,"Gsl")==0) return GslEnum;
+	      else if (strcmp(name,"Cuffey")==0) return CuffeyEnum;
 	      else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
 	      else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum;
@@ -874,9 +875,9 @@
 	      else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
 	      else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
-	      else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
          else stage=8;
    }
    if(stage==8){
-	      if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
+	      if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
+	      else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum;
 	      else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
 	      else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
@@ -997,9 +998,9 @@
 	      else if (strcmp(name,"Materials")==0) return MaterialsEnum;
 	      else if (strcmp(name,"Nodes")==0) return NodesEnum;
-	      else if (strcmp(name,"Contours")==0) return ContoursEnum;
          else stage=9;
    }
    if(stage==9){
-	      if (strcmp(name,"Parameters")==0) return ParametersEnum;
+	      if (strcmp(name,"Contours")==0) return ContoursEnum;
+	      else if (strcmp(name,"Parameters")==0) return ParametersEnum;
 	      else if (strcmp(name,"Vertices")==0) return VerticesEnum;
 	      else if (strcmp(name,"Results")==0) return ResultsEnum;
Index: issm/trunk-jpl/src/m/classes/hydrologysommers.m
===================================================================
--- issm/trunk-jpl/src/m/classes/hydrologysommers.m	(revision 21525)
+++ issm/trunk-jpl/src/m/classes/hydrologysommers.m	(revision 21526)
@@ -16,4 +16,5 @@
 		neumannflux     = NaN;
 		relaxation      = 0;
+		storage         = 0;
 	end
 	methods
@@ -33,4 +34,5 @@
 	      % Set under-relaxation parameter to be 1 (no under-relaxation of nonlinear iteration)	
 			self.relaxation=1;
+			self.storage=0;
 		end % }}}
 		function md = checkconsistency(self,md,solution,analyses) % {{{
@@ -51,4 +53,5 @@
 			md = checkfield(md,'fieldname','hydrology.spchead','size',[md.mesh.numberofvertices 1]);	
          md = checkfield(md,'fieldname','hydrology.relaxation','>=',0);	
+			md = checkfield(md,'fieldname','hydrology.storage','>=',0);
 		end % }}}
 		function disp(self) % {{{
@@ -64,4 +67,5 @@
 			fielddisplay(self,'spchead','water head constraints (NaN means no constraint) (m)');
 			fielddisplay(self,'relaxation','under-relaxation coefficient for nonlinear iteration');
+			fielddisplay(self,'storage','englacial storage coefficient (void ratio)');
 		end % }}}
 		function marshall(self,prefix,md,fid) % {{{
@@ -80,4 +84,5 @@
 			WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','spchead','format','DoubleMat','mattype',1);
          WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','relaxation','format','Double');
+			WriteData(fid,prefix,'object',self,'class','hydrology','fieldname','storage','format','Double');
 		end % }}}
 	end
