Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 14694)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 14695)
@@ -5870,6 +5870,5 @@
 	/* Intermediaries */
 	IssmDouble  D_scalar,Jdet;
-	IssmDouble  sediment_compressibility, sediment_porosity, sediment_thickness,
-		sediment_transmitivity, water_compressibility, rho_freshwater, gravity, dt;
+	IssmDouble 	sediment_transmitivity,dt;
 	IssmDouble  sediment_storing;
 	IssmDouble  xyz_list[NUMVERTICES][3];
@@ -5884,18 +5883,9 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
-	sediment_compressibility			= matpar->GetSedimentCompressibility();
-	sediment_porosity							= matpar->GetSedimentPorosity();
-	sediment_thickness						= matpar->GetSedimentThickness();
-	sediment_transmitivity				= matpar->GetSedimentTransmitivity();
-	water_compressibility					= matpar->GetWaterCompressibility();
-	rho_freshwater								= matpar->GetRhoFreshwater();
-	gravity												= matpar->GetG();
+
+	sediment_storing			 = matpar->GetSedimentStoring();
+	sediment_transmitivity = matpar->GetSedimentTransmitivity();
 	
 	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
-
-	/*Compute Storing coefficient from the above */
-
-	sediment_storing = rho_freshwater*gravity*sediment_porosity*sediment_thickness*
-		(water_compressibility+(sediment_compressibility/sediment_porosity));
 
 	/* Start looping on the number of gaussian points: */
@@ -6010,5 +6000,7 @@
 	IssmDouble Jdet;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble dt,gravity,scalar,water_head;
+	IssmDouble dt,scalar,water_head;
+	IssmDouble water_load;
+	IssmDouble sediment_storing;
 	IssmDouble basis[numdof];
 	GaussTria* gauss=NULL;
@@ -6018,11 +6010,16 @@
 
 	/*Retrieve all inputs and parameters*/
+
+	sediment_storing = matpar->GetSedimentStoring();
+
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	gravity = matpar->GetG();
+
+	Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
 	Input* old_wh_input=NULL; 
+	
 	if(reCast<bool,IssmDouble>(dt)){
 		old_wh_input=inputs->GetInput(SedimentHeadEnum); _assert_(old_wh_input);
 	}
-
+	
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
@@ -6034,6 +6031,7 @@
 		GetNodalFunctions(basis, gauss);
 
-		/*Gravity term*/
-		scalar = Jdet*gauss->weight*gravity;
+		/*Loading term*/
+		water_input->GetInputValue(&water_load,gauss);
+		scalar = Jdet*gauss->weight*water_load;
 		if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
 		for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
@@ -6042,5 +6040,5 @@
 		if(reCast<bool,IssmDouble>(dt)){
 			old_wh_input->GetInputValue(&water_head,gauss);
-			scalar = Jdet*gauss->weight*water_head;
+			scalar = Jdet*gauss->weight*water_head*sediment_storing;
 			for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
 		}
@@ -6091,4 +6089,20 @@
 void  Tria::InputUpdateFromSolutionHydrology(IssmDouble* solution){
 
+	int hydrology_model;
+  parameters->FindParam(&hydrology_model,HydrologyEnum);
+	
+  switch(hydrology_model){
+	case HydrologyshreveEnum:
+		return InputUpdateFromSolutionHydrologyShreve(solution);
+	case HydrologydcEnum:
+		return InputUpdateFromSolutionHydrologyDC(solution);
+	default:
+		_error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
+  }
+}
+/*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionHydrologyShreve{{{*/
+void  Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){
+
 	/*Intermediaries*/
 	const int   numdof         = NDOF1 *NUMVERTICES;
@@ -6108,4 +6122,28 @@
 	/*Add input to the element: */
 	this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values));
+
+	/*Free ressources:*/
+	xDelete<int>(doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDC{{{*/
+void  Tria::InputUpdateFromSolutionHydrologyDC(IssmDouble* solution){
+
+	/*Intermediaries*/
+	const int   numdof         = NDOF1 *NUMVERTICES;
+	int        *doflist        = NULL;
+	IssmDouble  values[numdof];
+
+	/*Get dof list: */
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+	/*Use the dof list to index into the solution vector: */
+	for(int i=0;i<numdof;i++){
+		values[i]=solution[doflist[i]];
+		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Add input to the element: */
+	this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
 
 	/*Free ressources:*/
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h	(revision 14694)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h	(revision 14695)
@@ -244,7 +244,9 @@
 		ElementVector* CreatePVectorHydrologyShreve(void);
 		ElementVector* CreatePVectorHydrologyDC(void);
-		void      CreateHydrologyWaterVelocityInput(void);
+		void    CreateHydrologyWaterVelocityInput(void);
 		void	  GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution);
 		void	  InputUpdateFromSolutionHydrology(IssmDouble* solution);
+		void	  InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
+		void	  InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
 		#endif
 		#ifdef _HAVE_BALANCED_
Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp	(revision 14694)
+++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp	(revision 14695)
@@ -315,17 +315,8 @@
 }		 
 /*}}}*/ 
-/*FUNCTION Matpar::GetSedimentCompressibility {{{*/
-IssmDouble Matpar::GetSedimentCompressibility(){
-	return sediment_compressibility;		 
-}		 
-/*}}}*/ 
-/*FUNCTION Matpar::GetSedimentPorosity {{{*/
-IssmDouble Matpar::GetSedimentPorosity(){
-	return sediment_porosity;		 
-}		 
-/*}}}*/ 
-/*FUNCTION Matpar::GetSedimentThickness {{{*/
-IssmDouble Matpar::GetSedimentThickness(){
-	return sediment_thickness;		 
+/*FUNCTION Matpar::GetSedimentStoring {{{*/
+IssmDouble Matpar::GetSedimentStoring(){
+	return this->rho_freshwater* this->g* this->sediment_porosity* this->sediment_thickness*
+    ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity));		 
 }		 
 /*}}}*/ 
@@ -334,9 +325,4 @@
 	return sediment_transmitivity;		 
 }		 
-/*}}}*/ 
-/*FUNCTION Matpar::GetWaterCompressibility {{{*/
-IssmDouble Matpar::GetWaterCompressibility(){
-	return water_compressibility;		 
-}		 
 /*}}}*/			 
 /*FUNCTION Matpar::TMeltingPoint {{{*/
Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h	(revision 14694)
+++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h	(revision 14695)
@@ -106,9 +106,6 @@
 		IssmDouble GetHydrologyCR();
 		IssmDouble GetHydrologyN();
-		IssmDouble GetSedimentCompressibility();
-		IssmDouble GetSedimentPorosity();
-		IssmDouble GetSedimentThickness();
+		IssmDouble GetSedimentStoring();
 		IssmDouble GetSedimentTransmitivity();
-		IssmDouble GetWaterCompressibility();
 		IssmDouble TMeltingPoint(IssmDouble pressure);
 		IssmDouble PureIceEnthalpy(IssmDouble pressure);