Merge pull request #50 from openstreams/glacier_params

Fix conversion of time-dependent glacier parameters and HBV time dependent parameters

Author: vers-w

tests/test_wflow_hbv2.py

@@ -64,11 +64,11 @@ def testapirunhr(self):
         )
 
         print("Checking  water budget ....")
-        self.assertAlmostEqual(0.0011249125109316083, my_data[:, 2].sum(), places=4)
+        self.assertAlmostEqual(0.001239627579707303, my_data[:, 2].sum(), places=4)
 
         my_data = np.genfromtxt(os.path.join(caseName, runId, "run.csv"), delimiter=",")
         print("Checking  discharge ....")
-        self.assertAlmostEqual(1671.0379650115967, my_data[:, 2].mean(), places=4)      
+        self.assertAlmostEqual(163.88664160221816, my_data[:, 2].mean(), places=4)      
 
 if __name__ == "__main__":
     unittest.main()

wflow/wflow_funcs.py

@@ -532,15 +532,15 @@ def glacierHBV(GlacierFrac,
     :ivar Temperature: Air temperature
     :ivar TT: Temperature threshold for ice melting
     :ivar Cfmax: Ice degree-day factor in mm/(°C/day)
-    :ivar G_SIfrac: Fraction of the snow part turned into ice each timestep
+    :ivar G_SIfrac: Fraction of the snow part turned into ice each daily timestep
     :ivar timestepsecs: Model timestep in seconds
     :ivar basetimestep: Model base timestep (86 400 seconds)
 
     :returns: Snow,Snow2Glacier,GlacierStore,GlacierMelt,
     """
 
     #Fraction of the snow transformed into ice (HBV-light model)
-    Snow2Glacier = G_SIfrac * Snow
+    Snow2Glacier = G_SIfrac * (timestepsecs / basetimestep) * Snow
 
     Snow2Glacier = pcr.ifthenelse(
         GlacierFrac > 0.0, Snow2Glacier, pcr.scalar(0.0)
@@ -554,7 +554,7 @@ def glacierHBV(GlacierFrac,
     GlacierStore = GlacierStore + Snow2Glacier
 
     PotMelt = pcr.ifthenelse(
-        Temperature > TT, Cfmax * (Temperature - TT), pcr.scalar(0.0)
+        Temperature > TT, Cfmax * (timestepsecs / basetimestep) * (Temperature - TT), pcr.scalar(0.0)
     )  # Potential snow melt, based on temperature
 
     GlacierMelt = pcr.ifthenelse(

wflow/wflow_hbv.py

@@ -726,15 +726,15 @@ def initial(self):
             subcatch,
             self.Soil,
             0.02307,
-        )  # Recession constant baseflow   #K4=0.07; BASEFLOW:LINEARRESERVOIR
+        ) * self.timestepsecs / self.basetimestep # Recession constant baseflow   #K4=0.07; BASEFLOW:LINEARRESERVOIR
         if self.SetKquickFlow:
             self.KQuickFlow = self.readtblDefault(
                 self.Dir + "/" + self.intbl + "/KQuickFlow.tbl",
                 self.LandUse,
                 subcatch,
                 self.Soil,
                 0.09880,
-            )  # recession rate at flow HQ     #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
+            ) * self.timestepsecs / self.basetimestep # recession rate at flow HQ     #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
             self.SUZ = self.readtblDefault(
                 self.Dir + "/" + self.intbl + "/SUZ.tbl",
                 self.LandUse,
@@ -748,22 +748,22 @@ def initial(self):
                 subcatch,
                 self.Soil,
                 0.3,
-            )  # K0
+            ) * self.timestepsecs / self.basetimestep # K0
         else:
             self.KHQ = self.readtblDefault(
                 self.Dir + "/" + self.intbl + "/KHQ.tbl",
                 self.LandUse,
                 subcatch,
                 self.Soil,
                 0.09880,
-            )  # recession rate at flow HQ     #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
+            ) * self.timestepsecs / self.basetimestep # recession rate at flow HQ     #KHQ=0.2; OUTFLOWUPPERZONE_NONLINEARRESERVOIR
             self.HQ = self.readtblDefault(
                 self.Dir + "/" + self.intbl + "/HQ.tbl",
                 self.LandUse,
                 subcatch,
                 self.Soil,
                 3.27000,
-            )  # high flow rate HQ for which recession rate of upper reservoir is known   #HQ=3.76;
+            ) * self.timestepsecs / self.basetimestep # high flow rate HQ for which recession rate of upper reservoir is known   #HQ=3.76;
             self.AlphaNL = self.readtblDefault(
                 self.Dir + "/" + self.intbl + "/AlphaNL.tbl",
                 self.LandUse,
@@ -778,7 +778,7 @@ def initial(self):
             subcatch,
             self.Soil,
             0.4000,
-        )  # percolation from Upper to Lowerzone (mm/day)
+        ) * self.timestepsecs / self.basetimestep # percolation from Upper to Lowerzone (mm/day)
         self.CFR = self.readtblDefault(
             self.Dir + "/" + self.intbl + "/CFR.tbl",
             self.LandUse,
@@ -814,7 +814,7 @@ def initial(self):
             subcatch,
             self.Soil,
             2.0,
-        )  # maximum capillary rise from runoff response routine to soil moisture routine
+        ) * self.timestepsecs / self.basetimestep # maximum capillary rise from runoff response routine to soil moisture routine
         self.ICF = self.readtblDefault(
             self.Dir + "/" + self.intbl + "/ICF.tbl",
             self.LandUse,
@@ -871,7 +871,7 @@ def initial(self):
             subcatch,
             self.Soil,
             3.75653,
-        )
+        ) * self.timestepsecs / self.basetimestep
         # WHC= 0.10000        # fraction of Snowvolume that can store water
         self.WHC = self.readtblDefault(
             self.Dir + "/" + self.intbl + "/WHC.tbl",