Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 24063)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.cpp	(revision 24064)
@@ -54,4 +54,11 @@
 	int e1 = iomodel->faces[4*index+2];
 	int e2 = iomodel->faces[4*index+3];
+
+	if(e2==-1){
+		this->boundary = true;
+	}
+	else{
+		this->boundary = false;
+	}
 
 	/*Set Element hook (4th column may be -1 for boundary edges)*/
@@ -545,7 +552,6 @@
 		if(Ngrad<1.e-12) Ngrad = 1.e-12;
 
-		/*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
-		IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
-		IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
+		/*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
+		IssmDouble Ks = ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.);
 
 		/*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
@@ -585,5 +591,5 @@
 	/*Initialize Element matrix and return if necessary*/
 	Tria*  tria=(Tria*)element;
-	if(!tria->IsIceInElement()){
+	if(!tria->IsIceInElement() || this->boundary){
 		this->S = 0.;
 		return;
@@ -674,5 +680,8 @@
 
 	/*Solve ODE*/
-	this->S = ODE1(alpha,beta,this->S,dt,0);
+	this->S = ODE1(alpha,beta,this->S,dt,2);
+
+	/*Make sure Area > 0*/
+	if(this->S<0.) this->S = 0.;
 
 	/*Clean up and return*/
Index: /issm/trunk-jpl/src/c/classes/Loads/Channel.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Loads/Channel.h	(revision 24063)
+++ /issm/trunk-jpl/src/c/classes/Loads/Channel.h	(revision 24064)
@@ -20,4 +20,5 @@
 	private: 
 		IssmDouble S;
+		bool       boundary;
 
 	public: