computewaterconductivity.c

compute Water Conductivity

double  computeWaterConductivity(double uc, double porosity ) 
{

	double waterDensity = 1.0e3 ;
	double kcc;
	
	if(aKw_Type==WaterCond_Van_Genuchten) //303
	{  //pressure head form
		double afa, n , ks ;
		double se , m ;

		afa=aParameters[0];	    n =aParameters[1] ;
		m = 1.0 - 1.0 / n ;
		ks =aParameters[2] ; 		

		double h = uc ;
		if(h>=0.0)
		{
			se = 1.0 ;
			kcc = ks ;

		}
		else
		
		{
			se = pow(-afa*h,n) ;
			se = pow(1.0+se,-m) ;
		kcc = pow(se,1.0/m);
		kcc = 1.0-pow(kcc,m) ;
		kcc = ks * sqrt(se) *kcc * kcc ;
		}

		
	}
	
	else if(aKw_Type==WaterCond_Brooks_Corey_1964) //308
	{
		double Sres,ks ;

		Sres = aParameters[0] ;
		ks= aParameters[1] ;
		double sw = this->computeWaterSaturation (uc) ;
		
		double se = (sw-Sres) / ( 1.0 - Sres) ;		
		if(se<0.0)
			se = 0.0 ;
		else if (se>1.0)
			se = 1.0 ; 

		kcc = ks * se * se* se ;  // conductivity (m/s)
	}

	else if(aKw_Type==WaterCond_Haverkamp_1977) //310
	{
		double kr,ks ;

		double C = aParameters[3] ;
		double D = aParameters[4];
		ks= aParameters[5] ;

		if(uc<-0.615)
			uc = -0.615 ;

		if(uc>=0.0)
		{
			kcc = ks ;
		}
		else
		{
			kr = C / (C+ pow(-100.0*uc,D) ) ;
			kcc = kr * ks ;
		}
		//	kcc /= waterDensity *9.8 ;
	}

	else if(aKw_Type==WaterCond_Kovacs)  //306
	{
        double A,B,C,D ;
		//Kw = A* pow(10,B*e)*( (Sw-C)/(1.0-C) )^D

		A = aParameters[0] ;
		B= aParameters[1] ;
		C = aParameters[2] ;
		D = aParameters[3] ;

		if (uc<0.0)
			uc = 0.0 ;

		double e6 =  B * porosity / (1.0 - porosity ); 
		double sw = this->computeWaterSaturation (uc) ;
		double tp = A*pow(10.0,e6 );
		kcc = tp * pow((sw-C)/(1.0-C)  ,D) ;
	}

	else if(aKw_Type==Touma_Vauclin_1986)  //312
	{
        double Aw,Bw,theta,D ;
		//Kw = A* pow(10,B*e)*( (Sw-C)/(1.0-C) )^D

		Aw = aParameters[0] ;
		Bw= aParameters[1] ;

		theta = this->computeWaterSaturation (uc) * 0.312;

		if (uc<0.0)
			uc = 0.0 ;

		kcc = Aw * pow(theta  ,Bw) ;
	//	kcc /= waterDensity *9.8 ;
	}

	return	kcc  ;

}