Initial changes to CEnthalpyModel::RouteMassinReservoir, Reservoir.cpp and Reservoir.h to introduce two lake layers to lake enthalpy model.

Author: mschnorb

src/EnergyTransport.cpp

@@ -270,17 +270,21 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub
   double A_new = pRes->GetSurfaceArea();
   double A_old = pRes->GetOldSurfaceArea();
   double A_avg = 0.5 * (A_new + A_old);
+  double V_h_new = pRes->GetHypolimnionStorage       ();
+  double V_h_old = pRes->GetOldHypolimnionStorage    ();
 
   CHydroUnit*   pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex());
 
   double Acorr=1.0;
 
   double SW(0), LW(0), LW_in(0), temp_air(0), AET(0);
+
   double hstar(0),ksed(0),Vsed=0.001;
   double T_new =ConvertVolumetricEnthalpyToTemperature(_aMres[p]      / V_new);
   double T_old =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p] / V_old);
   double Ts_new=ConvertVolumetricEnthalpyToTemperature(_aMsed[p]      / Vsed );
   double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p] / Vsed );
+ 
 
   if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input
 
@@ -300,7 +304,7 @@ double CEnthalpyModel::GetEnergyLossesFromLake(const int p, double &Q_sens, doub
   }
 
   Q_sens  =hstar* A_avg * (temp_air -0.5*(T_new+T_old));
-  Q_cond  =ksed * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old));
+  Q_cond  =kdiff * A_avg * (0.5*(Ts_new+Ts_old)- 0.5*(T_new+T_old));
   Q_sw_in =(SW      )*A_avg;
   Q_lw_in =(LW_in   )*A_avg;
   Q_lw_out=(LW      )*A_avg;
@@ -340,15 +344,27 @@ void   CEnthalpyModel::RouteMassInReservoir   (const int          p,          //
   double tstep=Options.timestep;
 
   if((V_old<=0.0) || (V_new<=0.0)) { Res_mass=ResSedMass=0.0; return;} //handles dried out reservoir/lake
-
+  
+  double V_h_new = pRes->GetHypolimnionStorage       ();
+  double V_h_old = pRes->GetOldHypolimnionStorage    ();
+  double V_e_new = V_new-V_h_new;
+  double V_e_old = V_old-V_h_old;
+  double Q_dn_new = 0;
+  double Q_dn_old = 0;
+  double Q_up_new = 0;
+  double Q_up_old = 0;
+  double A_h_new = pRes->GetMixingArea();
+  double A_h_old = pRes->GetOldMixingArea();
+  double A_h_avg = 0.5 * (A_h_new + A_h_old);
+   
   CHydroUnit*   pHRU=_pModel->GetHydroUnit(pRes->GetHRUIndex());
 
   double Acorr=1.0;
 
-  double T_old    =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_old);
+  double T_old    =ConvertVolumetricEnthalpyToTemperature(_aMres_last[p]/V_e_old);
 
   double SW(0), LW(0), LW_in(0), temp_air(0), AET(0);
-  double hstar(0), ksed(0), Vsed=0.001;
+  double hstar(0);
   if(pHRU!=NULL) { //otherwise only simulate advective mixing+ rain input
     Acorr    =pHRU->GetArea()*M2_PER_KM2/A_avg; //handles the fact that GetAET() returns mm/d normalized by HRU area, not actual area
 
@@ -359,14 +375,38 @@ void   CEnthalpyModel::RouteMassInReservoir   (const int          p,          //
     LW       =-STEFAN_BOLTZ*EMISS_WATER*pow(T_old+ZERO_CELSIUS,4); //[MJ/m2/d] //TMP DEBUG -time-lagged - should include in the N-R formulation
 
     AET      =pRes->GetAET()*Acorr/ MM_PER_METER ;             //[m/d]
-
-    Vsed     =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2;
+//  Vsed     =pRes->GetLakebedThickness() * pHRU->GetArea()*M2_PER_KM2;
     hstar    =pRes->GetLakeConvectionCoeff(); //[MJ/m2/d/K]
-    ksed     =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K]
+//  ksed     =pRes->GetLakebedConductivity()/0.5/ pRes->GetLakebedThickness();//[MJ/m2/d/K]
   }
 
-  double T_sed_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/Vsed);
-
+  double Ts_old=ConvertVolumetricEnthalpyToTemperature(_aMsed_last[p]/V_h_old);
+  
+  
+  double dens_e,dens_h,vdiff,kdiff;
+  
+  double a0= 999.842594;
+  double a1= 6.793952e-2;
+  double a2=-9.095290e-3;
+  double a3= 1.001685e-4;
+  double a4=-1.120083e-6;
+  double a5= 6.536332e-9;
+  
+  dens_e = a0+a1*T_old + a2*pow(T_old,2) + a3*pow(T_old,3) +a4*pow(T_old,4)+a5*pow(T_old,5);
+  dens_h = a0+a1*Ts_old + a2*pow(Ts_old,2) + a3*pow(Ts_old,3) +a4*pow(Ts_old,4)+a5*pow(Ts_old,5);
+  
+  double N2;
+  
+ if(dens_e<dens_h){
+	  N2 = -1*GRAVITY/dens_e*(dens_e-dens_h)/1.0;
+  }
+double m = -0.65;
+double b = -3.1;
+  
+  vdiff = pow(10,m*log10(N2)+b)/100/100;
+  kdiff = vdiff*HCP_WATER*SEC_PER_DAY;
+  
+  
   // N-R solution of Crank-nicolson problem as set of two non-linear algebraic equations
   // -----------------------------------------------------------------------------------------
   //dE   /dt=Qh_in-Qh_out+(Rnet*A+P*hrain*A)-ET*rho*LH*A+k*(Tair-T(E))*A+ksed*(Tsed-T(E))*A     //[MJ/d]
@@ -393,11 +433,11 @@ void   CEnthalpyModel::RouteMassInReservoir   (const int          p,          //
 
   B[0]+=(1.0-0.5*tstep/V_old*Q_old)*_aMres_last[p];
   B[0]-=0.5*tstep*hstar*A_old*T_old;
-  B[0]-=0.5*tstep*ksed *A_avg*T_old;
-  B[0]+=0.5*tstep*ksed *A_avg*T_sed_old;
+  B[0]-=0.5*tstep*kdiff *A_avg*T_old;
+  B[0]+=0.5*tstep*kdiff *A_avg*Ts_old;
 
-  B[1] =0.5*tstep*ksed *A_avg*T_old;
-  B[1]+=_aMsed_last[p]-0.5*tstep*ksed *A_avg*T_sed_old;
+  B[1] =0.5*tstep*kdiff *A_avg*T_old;
+  B[1]+=_aMsed_last[p]-0.5*tstep*kdiff *A_avg*Ts_old;
 
   int iter=0;
   double change=ALMOST_INF;
@@ -409,28 +449,28 @@ void   CEnthalpyModel::RouteMassInReservoir   (const int          p,          //
 
   while (change > tolerance)
   {
-    T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_new);
-    Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/Vsed);
+    T_guess = ConvertVolumetricEnthalpyToTemperature(E_guess/V_e_new);
+    Ts_guess= ConvertVolumetricEnthalpyToTemperature(Es_guess/V_h_new);
 
     A[0][0] = 1.0 + 0.5 * tstep  * Q_new/V_new;
     A[1][1] = 1.0;
     A[1][0] = 0.0;
     A[0][1] = 0.0;
     if (E_guess!=0.0){
       A[0][0]+= 0.5*tstep*hstar*A_new*T_guess/E_guess;
-      A[0][0]+= 0.5*tstep*ksed *A_avg*T_guess/E_guess;  //[MJ/m2/d/K]*[m2]/[MJ/m3/K] = [m3/d]
-      A[1][0]-= 0.5*tstep*ksed *A_avg*T_guess/E_guess;
+      A[0][0]+= 0.5*tstep*kdiff *A_avg*T_guess/E_guess;  //[MJ/m2/d/K]*[m2]/[MJ/m3/K] = [m3/d]
+      A[1][0]-= 0.5*tstep*kdiff *A_avg*T_guess/E_guess;
     }
     if (Es_guess!=0.0){
-      A[0][1]-= 0.5*tstep*ksed *A_avg*Ts_guess/Es_guess;
-      A[1][1]+= 0.5*tstep*ksed *A_avg*Ts_guess/Es_guess;
+      A[0][1]-= 0.5*tstep*kdiff *A_avg*Ts_guess/Es_guess;
+      A[1][1]+= 0.5*tstep*kdiff *A_avg*Ts_guess/Es_guess;
     }
 
     //J_ij=A_ij+E*dA_ij/dE = Jacobian
-    J[0][0]=A[0][0]+0.5*tstep*(hstar*A_new+ksed*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_new)/ V_new);
-    J[0][1]=A[0][1]-0.5*tstep*(            ksed*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/Vsed )/ Vsed );
-    J[1][0]=A[1][0]-0.5*tstep*(            ksed*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_new)/ V_new);
-    J[1][1]=A[1][1]+0.5*tstep*(            ksed*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/Vsed )/ Vsed );
+    J[0][0]=A[0][0]+0.5*tstep*(hstar*A_new+kdiff*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_e_new)/ V_e_new);
+    J[0][1]=A[0][1]-0.5*tstep*(            kdiff*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/V_h_new )/ V_h_new );
+    J[1][0]=A[1][0]-0.5*tstep*(            kdiff*A_avg) * (TemperatureEnthalpyDerivative(E_guess /V_e_new)/ V_e_new);
+    J[1][1]=A[1][1]+0.5*tstep*(            kdiff*A_avg) * (TemperatureEnthalpyDerivative(Es_guess/V_h_new )/ V_h_new );
 
     R[0]   =-A[0][0]*E_guess-A[0][1]*Es_guess+B[0];
     R[1]   =-A[1][0]*E_guess-A[1][1]*Es_guess+B[1];
@@ -440,7 +480,7 @@ void   CEnthalpyModel::RouteMassInReservoir   (const int          p,          //
     dE  = (J[0][1]*R[1]-J[1][1]*R[0])/den;
     dEs = (J[1][0]*R[0]-J[0][0]*R[1])/den;
 
-    change =sqrt(dE*dE+dEs*dEs)/ V_new / HCP_WATER; //convert to approx temp difference (for >0C water)
+    change =sqrt(dE*dE+dEs*dEs)/ V_e_new / HCP_WATER; //convert to approx temp difference (for >0C water)
 
     E_guess  +=dE;
     Es_guess +=dEs;

src/Makefile

@@ -25,12 +25,12 @@ LDFLAGS  :=
 CXXFLAGS += -std=c++11 -fPIC
 
 # OPTION 1) include netcdf - uncomment following two commands (assumes netCDF path = /usr/local): 
-#CXXFLAGS += -Dnetcdf     
-#LDLIBS   += -L/usr/local -lnetcdf  
+CXXFLAGS += -Dnetcdf     
+LDLIBS   += -L/usr/local -lnetcdf  
 
 # OPTION 2) include lp_solve for water management optimization- uncomment following two commands (assumes liblpsolve55.a in ../lib/lp_solve_unix/ folder):
-CXXFLAGS += -D_LPSOLVE_ 
-LDLIBS   += -L../lib/lp_solve_unix -llpsolve55
+#CXXFLAGS += -D_LPSOLVE_ 
+#LDLIBS   += -L../lib/lp_solve_unix -llpsolve55
 
 # OPTION 3) if you use a OSX/BSD system, uncomment the LDFLAGS line below
 # this is to allow for use a 1Gb stack, see http://linuxtoosx.blogspot.ca/2010/10/stack-overflow-increasing-stack-limit.html

src/Reservoir.cpp

@@ -18,7 +18,7 @@ void CReservoir::BaseConstructor(const string Name,const long SubID)
   _SBID=SubID;
 
   _pDownSB=NULL;
-
+ // _mixing_depth=1;
   _lakebed_thick=1.0;
   _lakebed_cond =0.0;
   _lake_convcoeff=2.0;
@@ -403,6 +403,14 @@ long  CReservoir::GetSubbasinID          () const { return _SBID; }
 //
 double  CReservoir::GetStorage           () const { return GetVolume(_stage); }
 
+/// \returns reservoir hypoliminion storage [m3]
+//
+double  CReservoir::GetHypolimnionStorage           () const { return GetVolume(_stage-_mixing_depth); }
+
+/// \returns reservoir old hypoliminion storage [m3]
+//
+double  CReservoir::GetOldHypolimnionStorage           () const { return GetVolume(_stage_last-_mixing_depth); }
+
 //////////////////////////////////////////////////////////////////
 /// \returns current stage [m]
 //
@@ -449,6 +457,16 @@ double  CReservoir::GetSurfaceArea       () const { return GetArea(_stage); }
 //
 double  CReservoir::GetOldSurfaceArea    () const { return GetArea(_stage_last); }
 
+//////////////////////////////////////////////////////////////////
+/// \returns current mixing area [m2]
+//
+double  CReservoir::GetMixingArea       () const { return GetArea(_stage-_mixing_depth); }
+
+//////////////////////////////////////////////////////////////////
+/// \returns current old mixing area [m2]
+//
+double  CReservoir::GetOldMixingArea       () const { return GetArea(_stage_last-_mixing_depth); }
+
 //////////////////////////////////////////////////////////////////
 /// \returns lakebed thickness [m]
 //

src/Reservoir.h

@@ -55,6 +55,7 @@ class CReservoir
   string       _name;                ///< reservoir name
   long         _SBID;                ///< subbasin ID
   double       _max_capacity;        ///< maximum reservoir capacity [m3]
+  double       _mixing_depth;
 
   double       _lakebed_thick;       ///< lakebed thickness [m]
   double       _lakebed_cond;        ///< lakebed thermal conductivity [MJ/m/K/d]
@@ -178,6 +179,8 @@ class CReservoir
 
   double            GetStorage               () const; //[m3]
   double            GetOldStorage            () const; //[m3]
+  double            GetHypolimnionStorage    () const; //[m3]
+  double            GetOldHypolimnionStorage () const; //[m3]
   double            GetOutflowRate           () const; //[m3/s]
   double            GetOldOutflowRate        () const; //[m3/s]
   double            GetIntegratedOutflow     (const double &tstep) const; //[m3]
@@ -191,6 +194,8 @@ class CReservoir
   double            GetOldStage              () const; //[m]
   double            GetSurfaceArea           () const; //[m2]
   double            GetOldSurfaceArea        () const; //[m2]
+  double            GetMixingArea            () const; //[m2]
+  double            GetOldMixingArea            () const; //[m2]
   double            GetLakebedThickness      () const; //[m]
   double            GetLakebedConductivity   () const; //[MJ/m/K/d]
   double            GetLakeConvectionCoeff   () const; //[MJ/m2/d/K]