Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 14659)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 14660)
@@ -5758,11 +5758,11 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixHydrologyShreve{{{*/
+/*FUNCTION Tria::CreateKMatrixHydrologyShreve(void){{{*/
 ElementMatrix* Tria::CreateKMatrixHydrologyShreve(void){
-
+	
 	/*Constants*/
 	const int  numdof=NDOF1*NUMVERTICES;
-
-	/*Intermediaries */
+	
+/*Intermediaries */
 	IssmDouble diffusivity;
 	IssmDouble Jdettria,DL_scalar,dt,h;
@@ -5779,15 +5779,15 @@
 	IssmDouble DLprime[2][2]                   ={0.0};
 	GaussTria *gauss=NULL;
-
-	/*Skip if water or ice shelf element*/
+	
+/*Skip if water or ice shelf element*/
 	if(IsOnWater() | IsFloating()) return NULL;
-
-	/*Initialize Element matrix*/
+	
+/*Initialize Element matrix*/
 	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
-	/*Create water velocity vx and vy from current inputs*/
+	
+/*Create water velocity vx and vy from current inputs*/
 	CreateHydrologyWaterVelocityInput();
-
-	/*Retrieve all inputs and parameters*/
+	
+/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
@@ -5796,48 +5796,48 @@
 	Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
 	h=sqrt(2*this->GetArea());
-
-	/* Start  looping on the number of gaussian points: */
+	
+/* Start  looping on the number of gaussian points: */
 	gauss=new GaussTria(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
+		
 		gauss->GaussPoint(ig);
-
+		
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
 		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
+		
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
 		vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
 		vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
-
+		
 		DL_scalar=gauss->weight*Jdettria;
-
+		
 		TripleMultiply( &L[0],1,numdof,1,
-					&DL_scalar,1,1,0,
-					&L[0],1,numdof,0,
-					&Ke->values[0],1);
-
+										&DL_scalar,1,1,0,
+										&L[0],1,numdof,0,
+										&Ke->values[0],1);
+		
 		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
 		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
-
+		
 		dvxdx=dvx[0];
 		dvydy=dvy[1];
 		DL_scalar=dt*gauss->weight*Jdettria;
-
+		
 		DL[0][0]=DL_scalar*dvxdx;
 		DL[1][1]=DL_scalar*dvydy;
 		DLprime[0][0]=DL_scalar*vx;
 		DLprime[1][1]=DL_scalar*vy;
-
+		
 		TripleMultiply( &B[0][0],2,numdof,1,
-					&DL[0][0],2,2,0,
-					&B[0][0],2,numdof,0,
-					&Ke->values[0],1);
-
+										&DL[0][0],2,2,0,
+										&B[0][0],2,numdof,0,
+										&Ke->values[0],1);
+		
 		TripleMultiply( &B[0][0],2,numdof,1,
-					&DLprime[0][0],2,2,0,
-					&Bprime[0][0],2,numdof,0,
-					&Ke->values[0],1);
-
+										&DLprime[0][0],2,2,0,
+										&Bprime[0][0],2,numdof,0,
+										&Ke->values[0],1);
+		
 		/*Artificial diffusivity*/
 		vel=sqrt(vx*vx+vy*vy);
@@ -5850,12 +5850,12 @@
 		KDL[0][1]=DL_scalar*K[0][1];
 		KDL[1][1]=DL_scalar*K[1][1];
-
+		
 		TripleMultiply( &Bprime[0][0],2,numdof,1,
-					&KDL[0][0],2,2,0,
-					&Bprime[0][0],2,numdof,0,
-					&Ke->values[0],1);
-	}
-
-	/*Clean up and return*/
+										&KDL[0][0],2,2,0,
+										&Bprime[0][0],2,numdof,0,
+										&Ke->values[0],1);
+	}
+	
+/*Clean up and return*/
 	delete gauss;
 	return Ke;
@@ -5870,10 +5870,12 @@
 	/* Intermediaries */
 	IssmDouble  D_scalar,Jdet;
-	IssmDouble  sediment_porosity,dt;
+	IssmDouble  sediment_compressibility, sediment_porosity, sediment_thickness,
+		sediment_transmitivity, water_compressibility, rho_freshwater, gravity, dt;
+	IssmDouble  sediment_storing;
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  B[2][numdof];
-	IssmDouble L[NUMVERTICES];
+	IssmDouble  L[NUMVERTICES];
 	IssmDouble  D[2][2];
-	GaussTria  *gauss = NULL;
+	GaussTria   *gauss = NULL;
 
 	/*Initialize Element matrix*/
@@ -5882,6 +5884,18 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
-	sediment_porosity = matpar->GetSedimentPorosity();
+	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();
+	
 	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+
+	/*Compute Storing coefficient fro 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: */
@@ -5900,7 +5914,7 @@
 		GetBHydro(&B[0][0],&xyz_list[0][0],gauss); 
 		TripleMultiply(&B[0][0],2,numdof,1,
-					&D[0][0],2,2,0,
-					&B[0][0],2,numdof,0,
-					&Ke->values[0],1);
+									 &D[0][0],2,2,0,
+									 &B[0][0],2,numdof,0,
+									 &Ke->values[0],1);
 
 		/*Transient*/
@@ -5908,9 +5922,9 @@
 			GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
 			D_scalar=gauss->weight*Jdet;
-
+			
 			TripleMultiply(&L[0],numdof,1,0,
-						&D_scalar,1,1,0,
-						&L[0],1,numdof,0,
-						&Ke->values[0],1);
+										 &D_scalar,1,1,0,
+										 &L[0],1,numdof,0,
+										 &Ke->values[0],1);
 		}
 	}
Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h	(revision 14659)
+++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h	(revision 14660)
@@ -111,5 +111,4 @@
 		IssmDouble GetSedimentTransmitivity();
 		IssmDouble GetWaterCompressibility();
-		IssmDouble GetWaterDensity();
 		IssmDouble TMeltingPoint(IssmDouble pressure);
 		IssmDouble PureIceEnthalpy(IssmDouble pressure);