v489 - [BENCHMARKED] - support for frozen wetlands, -template argument, TOTAL_SWE, convolution with transport

Freezing of wetlands (FrozenLake.cpp; PartitionPrecip.cpp)

Support for RVP template generation from command line (ParseInput.cpp; RavenMain.cpp)

support for convolution + transport (StateVariable.h/.cpp; RavenInclude.h)

Added TOTAL_SWE state variable tracking (idea from G. Jost) (RavenInclude.h; Solvers.cpp)

added exit gracefully if SAT_WILT>FIELD_CAPACITY (HydroUnits.cpp)

handle no name of :OutflowControlStructure (ParseHRUFile.cpp)

proper handling of Lakes and wetlands in SNOBAL_HBV (SnowBalance.cpp)

Author: James Craig

src/CommonFunctions.cpp

@@ -50,7 +50,7 @@ string GetProcessName(process_type p)
   case(SNOWSQUEEZE):        {name="Liquid snow release";      break;}
   case(REFREEZE):           {name="Snow Refreeze";            break;}
   case(SUBLIMATION):        {name="Sublimation";              break;}
-  case(SNOW_BALANCE):       {name="Snow Melt & Refreeze";     break;}
+  case(SNOW_BALANCE):       {name="Snow Balance";             break;}
   case(GLACIER_MELT):       {name="Glacier Melt";             break;}
   case(GLACIER_RELEASE):    {name="Glacier Release";          break;}
   case(GLACIER_INFIL):      {name="Glacier Infiltration";     break;}

src/ControlStructures.cpp

@@ -283,7 +283,7 @@ bool      COutflowRegime::AreConditionsMet(const time_struct& tt) const
 //////////////////////////////////////////////////////////////////
 /// \brief returns outflow from control structure for given stage, with prioritized constraints applied
 /// \param h [in] stage, in [m]
-/// \param Q_start [in] total outflow at start of timestep [m3/s]
+/// \param Q_start [in] control structure outflow at start of timestep [m3/s]
 /// \returns outflow, in [m3/s]
 //
 double          COutflowRegime::GetOutflow(const double &h, const double &h_start, const double &Q_start, const long &target_SBID, const double &drefelev) const
@@ -293,7 +293,7 @@ double          COutflowRegime::GetOutflow(const double &h, const double &h_star
 
   double Q = _pCurve->GetDischarge(h, h_start, Q_start, rivdepth, drefelev);
 
-  double dQdt = (Q-Q_start)/tstep;
+  double dQdt = (Q-Q_start)/tstep; //m3/s/d
 
   //apply constraints here
   for (int j = _nConstraints-1; j >=0; j--) //priority from first to last

src/ControlStructures.h

@@ -78,7 +78,7 @@ class COutflowRegime
   void   AddRegimeConstraint(RegimeConstraint *pCond);
 
   bool      AreConditionsMet(const time_struct &tt) const;
-  double          GetOutflow(const double &h, const double &Qstart, const double &hstart, const long &target_SBID, const double &drefelev) const;
+  double          GetOutflow(const double &h, const double &hstart, const double &Qstart, const long &target_SBID, const double &drefelev) const;
 };
 
 

src/FrozenLake.cpp

@@ -15,9 +15,9 @@
 
 //////////////////////////////////////////////////////////////////
 /// \brief Implementation of the standard constructor
-/// \param cr_type [in] Model of capillary rise selected
-/// \param In_index [in] Soil storage unit index from which water is lost
-/// \param Out_index [in] Soil storage unit index to which water rises
+/// \param lakefreeze_type [in] Model of lake freeze selected
+/// \param pTransMod [in] pointer to transport model
+/// \param pModel [in] - ponter to hydrologic model
 //
 CmvFrozenLake::CmvFrozenLake(lakefreeze_type       type,
                              const CTransportModel *pTransMod,
@@ -63,11 +63,10 @@ void CmvFrozenLake::GetParticipatingParamList(string  *aP , class_type *aPC , in
   {
     nP=1;
     aP [0]="LAKESNOW_BUFFER_HT";     aPC [0]=CLASS_LANDUSE;
-
   }
   else if (_type == LFREEZE_THERMAL)
   {
-    ExitGracefully("LFREEZE_THERMAL",STUB); //need to allow snow buffering
+    
   }
   else
   {
@@ -106,7 +105,8 @@ void   CmvFrozenLake::GetRatesOfChange( const double      *state_vars,
                                           const time_struct &tt,
                                                 double      *rates) const
 {
-  if (pHRU->GetHRUType()!=HRU_LAKE){return;}//Lakes  only (includes reservoirs)
+  if ((pHRU->GetHRUType()!=HRU_LAKE) && 
+      (pHRU->GetHRUType()!=HRU_WETLAND)){return;}//Lakes and wetlands only (includes reservoirs)
 
   double ice_thick=state_vars[iFrom[0]];
 
@@ -130,6 +130,7 @@ void   CmvFrozenLake::GetRatesOfChange( const double      *state_vars,
   else if (_type == LFREEZE_THERMAL)
   {
     //just queries thermal model and updates ice thickness variable
+    // \todo[funct] need to allow snow buffering
 
     ice_thick=state_vars[iFrom[0]];
 
@@ -142,7 +143,6 @@ void   CmvFrozenLake::GetRatesOfChange( const double      *state_vars,
 
     rates[0]=(new_ice_thick-ice_thick)/Options.timestep;
   }
-
 }
 
 //////////////////////////////////////////////////////////////////

src/HydroUnits.cpp

@@ -278,7 +278,11 @@ double CHydroUnit::GetSoilTensionStorageCapacity(const int m) const
   ExitGracefullyIf((m<0) || (m>=_pModel->GetNumSoilLayers()),
                    "CHydroUnit GetSoilTensionStorageCapacity::improper index",BAD_DATA);
 #endif
-  return aThickness[m]*MM_PER_METER*_pSoil[m]->cap_ratio*(_pSoil[m]->field_capacity-_pSoil[m]->sat_wilt);
+  double storcap=aThickness[m]*MM_PER_METER*_pSoil[m]->cap_ratio*(_pSoil[m]->field_capacity-_pSoil[m]->sat_wilt);
+  if (storcap<-REAL_SMALL){
+    ExitGracefully("Wilting point saturation (SAT_WILT) is greater than field capacity (FIELD_CAPACITY) for one or more soil classes, leading to negative tension storage",BAD_DATA);
+  }
+  return storcap;
 }
 
 //////////////////////////////////////////////////////////////////

src/LandUseClass.cpp

@@ -508,6 +508,7 @@ void  CLandUseClass::SetSurfaceProperty(surface_struct &S,
   else if (!name.compare("PET_LIN_COEFF"          )){S.pet_lin_coeff=value;}
   else if (!name.compare("PET_VAP_COEFF"          )){S.pet_vap_coeff=value;}
   else if (!name.compare("RELHUM_CORR"            )){S.relhum_corr=value;}
+  else if (!name.compare("WINDVEL_CORR"           )){S.wind_vel_corr=value;}
   else if (!name.compare("WIND_VEL_CORR"          )){S.wind_vel_corr=value;}
   else if (!name.compare("GR4J_X4"                )){S.GR4J_x4=value;}
   else if (!name.compare("UBC_ICEPT_FACTOR"       )){S.UBC_icept_factor=value;}
@@ -614,6 +615,7 @@ double CLandUseClass::GetSurfaceProperty(const surface_struct &S, string param_n
   else if (!name.compare("PET_LIN_COEFF"          )){return S.pet_lin_coeff;}
   else if (!name.compare("PET_VAP_COEFF"          )){return S.pet_vap_coeff;}
   else if (!name.compare("RELHUM_CORR"            )){return S.relhum_corr;}
+  else if (!name.compare("WINDVEL_CORR"           )){return S.wind_vel_corr;}
   else if (!name.compare("WIND_VEL_CORR"          )){return S.wind_vel_corr;}
   else if (!name.compare("GR4J_X4"                )){return S.GR4J_x4;}
   else if (!name.compare("UBC_ICEPT_FACTOR"       )){return S.UBC_icept_factor;}

src/ParseHRUFile.cpp

@@ -1657,6 +1657,9 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long
       if (pContStruct != NULL) {
         ExitGracefully("ReservoirParse: new control structure started before finishing earlier one with :EndOutflowControlStructure",BAD_DATA_WARN);
       }
+
+      string name="unnamed";
+      if (Len>1){name=s[1];}
       long downID=pModel->GetSubBasinByID(SBID)->GetDownstreamID(); //default target basin
       pContStruct=new CControlStructure(s[1],SBID,downID);//assumes SBID appears first
     }

src/ParseInput.cpp

@@ -108,8 +108,10 @@ bool ParseInputFiles (CModel      *&pModel,
     }
     ExitGracefully("Cannot find or read .rvi file",BAD_DATA);return false;
   }
-  if (!ParseNetCDFRunInfoFile(pModel, Options, runname_overridden,runmode_overridden)){
-    ExitGracefully("Cannot find or read NetCDF runinfo file", BAD_DATA); return false;
+  if (!Options.create_rvp_template) {//otherwise, jump right to parse rvp, where template is created  
+    if (!ParseNetCDFRunInfoFile(pModel, Options, runname_overridden,runmode_overridden)){
+      ExitGracefully("Cannot find or read NetCDF runinfo file", BAD_DATA); return false;
+    }
   }
 
   // Class Property file (.rvp)
@@ -3531,15 +3533,20 @@ bool ParseMainInputFile (CModel     *&pModel,
   } //end while (!end_of_file)
   INPUT.close();
 
+  // Add TOTAL_SWE state variable if any snow is simulated 
+  if (pModel->GetStateVarIndex(SNOW) != -1) {
+    tmpS[0] = TOTAL_SWE; tmpLev[0]=0; tmpN=1;
+    pModel->AddStateVariables(tmpS,tmpLev,tmpN);
+  }
+  
+
   //===============================================================================================
   //Check input quality
   //===============================================================================================
   ExitGracefullyIf(Options.timestep<=0,
                    "ParseMainInputFile::Must have a postitive time step",BAD_DATA);
   ExitGracefullyIf(Options.duration<0,
                    "ParseMainInputFile::Model duration less than zero. Make sure :EndDate is after :StartDate.",BAD_DATA_WARN);
-  ExitGracefullyIf((pModel->GetStateVarIndex(CONVOLUTION,0)!=DOESNT_EXIST) && (pModel->GetTransportModel()->GetNumConstituents()>0),
-                   "ParseMainInputFile: cannot currently perform transport with convolution processes",BAD_DATA);
 
   if((Options.nNetCDFattribs>0) && (Options.output_format!=OUTPUT_NETCDF)){
     WriteAdvisory("ParseMainInputFile: NetCDF attributes were specified but output format is not NetCDF.",Options.noisy);

src/ParsePropertyFile.cpp

@@ -217,7 +217,7 @@ bool ParseClassPropertiesFile(CModel         *&pModel,
     }
   }
 
-  if (Options.create_rvp_template){
+  if (Options.create_rvp_template){ //can proceed even if .rvp doesnt exist.
     CreateRVPTemplate(aPmaster,aPCmaster,nPmaster,Options);
     ExitGracefully("Template RVP File Created.",SIMULATION_DONE);
   }

src/PartitionPrecip.cpp

@@ -56,12 +56,15 @@ void CmvPrecipitation::Initialize()
     if (iFrom[q]==DOESNT_EXIST){iFrom[q]=iAtmos; } //since this routine reacts to presence or absence of individual state variables
   }
 
-  //lake snow handling
+  //frozen lake & wetland snow handling
   if (pModel->StateVarExists(ICE_THICKNESS))
   {
     iFrom[13]=pModel->GetStateVarIndex(SNOW);         iTo[13]=pModel->GetLakeStorageIndex();
     iFrom[14]=pModel->GetStateVarIndex(SNOW_LIQ);     iTo[14]=pModel->GetLakeStorageIndex();
     iFrom[15]=pModel->GetStateVarIndex(PONDED_WATER); iTo[15]=pModel->GetLakeStorageIndex();
+    iFrom[16]=pModel->GetStateVarIndex(SNOW);         iTo[16]=pModel->GetStateVarIndex(DEPRESSION);
+    iFrom[17]=pModel->GetStateVarIndex(SNOW_LIQ);     iTo[17]=pModel->GetStateVarIndex(DEPRESSION);
+    iFrom[18]=pModel->GetStateVarIndex(PONDED_WATER); iTo[18]=pModel->GetStateVarIndex(DEPRESSION);
   }
 }
 //////////////////////////////////////////////////////////////////
@@ -167,9 +170,9 @@ void CmvPrecipitation::GetRatesOfChange(const double             *state_vars,
   int qSnowToLake=13;
   int qSnLiqToLake=14;
   int qPondToLake=15;
-  //int qSnowToDep=16;
-  //int qSnLiqToDep=17;
-  //int qPondToDep=18;
+  int qSnowToDep=16;
+  int qSnLiqToDep=17;
+  int qPondToDep=18;
 
   //get precipitation information from gauges
   total_precip=pHRU->GetForcingFunctions()->precip;//[mm/day]
@@ -237,12 +240,15 @@ void CmvPrecipitation::GetRatesOfChange(const double             *state_vars,
     rainthru = rainfall;
   }
 
-  bool lake_frozen=false;
+  bool lake_frozen   =false;
   bool wetland_frozen=false;
   if (pModel->StateVarExists(ICE_THICKNESS))
   {
     int iIceThick=pModel->GetStateVarIndex(ICE_THICKNESS);
-    if (state_vars[iIceThick]>10){lake_frozen=true;}
+    if (state_vars[iIceThick]>10){ //10mm of ice can support snow?
+      lake_frozen   =true;
+      wetland_frozen=true;
+    }
   }
 
   if ((pHRU->IsLake()) && (!lake_frozen))
@@ -254,7 +260,7 @@ void CmvPrecipitation::GetRatesOfChange(const double             *state_vars,
       rates[qLake]=snowthru+rainthru; // in lake, all water just added
     }
 
-    if (pModel->StateVarExists(ICE_THICKNESS))//all remaining snow and ponded water dumped to lake
+    if (pModel->StateVarExists(ICE_THICKNESS))//all remaining snow and ponded water dumped to lake \todo[funct]: handle multilayer snowpack models
     {
        rates[qSnowToLake ]+=state_vars[iSWE  ]/Options.timestep;
        rates[qSnLiqToLake]+=state_vars[iSnLiq]/Options.timestep;
@@ -264,11 +270,15 @@ void CmvPrecipitation::GetRatesOfChange(const double             *state_vars,
   {
     rates[qDep] =snowthru + rainthru; //all water goes to depression in wetland
 
-    /*if (pModel->StateVarExists(ICE_THICKNESS))// \todo[funct]: all remaining snow and ponded water dumped to wetland/depression
+    if (pModel->StateVarExists(ICE_THICKNESS))//all remaining snow and ponded water dumped to wetland/depression
     {
        rates[qSnowToDep ]+=state_vars[iSWE  ]/Options.timestep;
        rates[qSnLiqToDep]+=state_vars[iSnLiq]/Options.timestep;
-    }*/
+    }
+  }
+  else if (pHRU->GetHRUType() == HRU_WATER) 
+  {
+    rates[qSW] =snowthru + rainthru; //all water goes to surface water
   }
   else //regular land, frozen lakes or frozen wetlands - treat snow
   {
@@ -307,7 +317,7 @@ void CmvPrecipitation::GetRatesOfChange(const double             *state_vars,
       rates[qNewSnow]=snowthru;
     }
 
-    if ( pHRU->IsLake())
+    if (pHRU->IsLake())
     {
       if (pHRU->IsLinkedToReservoir()) {
         rates[qSW]=rainthru; //handles case when both reservoirs and lake storage used
@@ -325,13 +335,13 @@ void CmvPrecipitation::GetRatesOfChange(const double             *state_vars,
     {
       rates[qDep] =rainthru; //ponded water goes to depression in wetland
 
-      /*if (pModel->StateVarExists(ICE_THICKNESS))
+      if (pModel->StateVarExists(ICE_THICKNESS))
       {
-         rates[qPondToDep ]+=state_vars[iPond ]/Options.timestep; // \todo[funct]: handles on-wetland snowmelt whether ice is present or not
-      }*/
+         rates[qPondToDep ]+=state_vars[iPond ]/Options.timestep; //handles on-wetland snowmelt whether ice is present or not
+      }
     }
     else{
-      //everything remaining goes to ponded water variable, usually handled by infiltration & abstraction algorithms
+      //everything remaining goes to ponded water variable, usually handled by infiltration and/or abstraction algorithms
       rates[qPond]=rainthru;
     }
   }