Index: wflow-py/wflow/wflow_hbv.py =================================================================== diff -u -r4f8d31cc84a0ba96e9ef23399d92cbe5517aec2d -re81cbd87d29f5a1243914ad813d58883312d2938 --- wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision 4f8d31cc84a0ba96e9ef23399d92cbe5517aec2d) +++ wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision e81cbd87d29f5a1243914ad813d58883312d2938) @@ -1091,319 +1091,212 @@ self.wf_updateparameters() # read forcing an dynamic parameters self.Precipitation = max(0.0, self.Precipitation) * self.Pcorr - # self.Precipitation=cover(self.wf_readmap(self.P_mapstack,0.0),0.0) * self.Pcorr - # self.PotEvaporation=cover(self.wf_readmap(self.PET_mapstack,0.0),0.0) - # self.Inflow=cover(self.wf_readmap(self.Inflow_mapstack,0.0,verbose=False),0.0) - # These ar ALWAYS 0 at present!!! - # self.Inflow=pcrut.readmapSave(self.Inflow_mapstack,0.0) - if self.ExternalQbase: - self.Seepage = cover(self.wf_readmap(self.Seepage_mapstack, 0.0), 0.0) - else: - self.Seepage = cover(0.0) - self.Temperature = cover(self.wf_readmap(self.TEMP_mapstack, 10.0), 10.0) - self.Temperature = self.Temperature + self.TempCor + #Water onto the canopy + Interception=min(self.Precipitation,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep + self.InterceptionStorage=self.InterceptionStorage+Interception #: Current interception storage + self.Precipitation=self.Precipitation-Interception - # Multiply input parameters with a factor (for calibration etc) -p option in command line (no also in ini) - self.wf_multparameters() + self.PotEvaporation=exp(-self.EPF*self.Precipitation)*self.ECORR * self.PotEvaporation # correction for potential evaporation on wet days + self.PotEvaporation=self.CEVPF*self.PotEvaporation # Correct per landuse - RainFrac = ifthenelse( - 1.0 * self.TTI == 0.0, - ifthenelse(self.Temperature <= self.TT, scalar(0.0), scalar(1.0)), - min( - (self.Temperature - (self.TT - self.TTI / 2.0)) / self.TTI, scalar(1.0) - ), - ) - RainFrac = max( - RainFrac, scalar(0.0) - ) # fraction of precipitation which falls as rain - SnowFrac = 1.0 - RainFrac # fraction of self.Precipitation which falls as snow + self.IntEvap=min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage + self.InterceptionStorage=self.InterceptionStorage-self.IntEvap - self.Precipitation = ( - self.SFCF * SnowFrac * self.Precipitation - + self.RFCF * RainFrac * self.Precipitation - ) # different correction for rainfall and snowfall + # I nthe origal HBV code + RestEvap = max(0.0,self.PotEvaporation-self.IntEvap) - # Water onto the canopy - Interception = min( - self.Precipitation, self.ICF - self.InterceptionStorage - ) #: Interception in mm/timestep - self.InterceptionStorage = ( - self.InterceptionStorage + Interception - ) #: Current interception storage - self.Precipitation = self.Precipitation - Interception + if hasattr(self, 'ReserVoirComplexLocs'): + self.ReserVoirPotEvap = self.PotEvaporation + self.ReserVoirPrecip = self.Precipitation - self.PotEvaporation = ( - exp(-self.EPF * self.Precipitation) * self.ECORR * self.PotEvaporation - ) # correction for potential evaporation on wet days - self.PotEvaporation = self.CEVPF * self.PotEvaporation # Correct per landuse + self.PotEvaporation = self.filter_P_PET * self.PotEvaporation + self.Precipitation = self.filter_P_PET * self.Precipitation - self.IntEvap = min( - self.InterceptionStorage, self.PotEvaporation - ) #: Evaporation from interception storage - self.InterceptionStorage = self.InterceptionStorage - self.IntEvap - # I nthe origal HBV code - RestEvap = max(0.0, self.PotEvaporation - self.IntEvap) + SnowFall=SnowFrac*self.Precipitation #: snowfall depth + RainFall=RainFrac*self.Precipitation #: rainfall depth + PotSnowMelt=ifthenelse(self.Temperature > self.TT,self.Cfmax*(self.Temperature-self.TT),scalar(0.0)) #Potential snow melt, based on temperature + PotRefreezing=ifthenelse(self.Temperature < self.TT, self.Cfmax*self.CFR*(self.TT-self.Temperature),0.0) #Potential refreezing, based on temperature - if hasattr(self, "ReserVoirComplexLocs"): - self.ReserVoirPotEvap = self.PotEvaporation - self.ReserVoirPrecip = self.Precipitation + Refreezing=ifthenelse(self.Temperature < self.TT,min(PotRefreezing,self.FreeWater),0.0) #actual refreezing + self.SnowMelt=min(PotSnowMelt,self.DrySnow) #actual snow melt + self.DrySnow=self.DrySnow+SnowFall+Refreezing-self.SnowMelt #dry snow content + self.FreeWater=self.FreeWater-Refreezing #free water content in snow + MaxFreeWater=self.DrySnow*self.WHC + self.FreeWater=self.FreeWater+self.SnowMelt+RainFall + InSoil = max(self.FreeWater-MaxFreeWater,0.0) #abundant water in snow pack which goes into soil + self.FreeWater=self.FreeWater-InSoil + RainAndSnowmelt = RainFall + self.SnowMelt - self.PotEvaporation = self.filter_P_PET * self.PotEvaporation - self.Precipitation = self.filter_P_PET * self.Precipitation + self.SnowCover = ifthenelse(self.DrySnow >0, scalar(1), scalar(0)) + self.NrCell= areatotal(self.SnowCover,self.TopoId) - SnowFall = SnowFrac * self.Precipitation #: snowfall depth - RainFall = RainFrac * self.Precipitation #: rainfall depth - PotSnowMelt = ifthenelse( - self.Temperature > self.TT, - self.Cfmax * (self.Temperature - self.TT), - scalar(0.0), - ) # Potential snow melt, based on temperature - PotRefreezing = ifthenelse( - self.Temperature < self.TT, - self.Cfmax * self.CFR * (self.TT - self.Temperature), - 0.0, - ) # Potential refreezing, based on temperature + #first part of precipitation is intercepted + #Interception=min(InSoil,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep + #self.InterceptionStorage=self.InterceptionStorage+Interception #: Current interception storage + #NetInSoil=InSoil-Interception + NetInSoil=InSoil - Refreezing = ifthenelse( - self.Temperature < self.TT, min(PotRefreezing, self.FreeWater), 0.0 - ) # actual refreezing - self.SnowMelt = min(PotSnowMelt, self.DrySnow) # actual snow melt - self.DrySnow = ( - self.DrySnow + SnowFall + Refreezing - self.SnowMelt - ) # dry snow content - self.FreeWater = self.FreeWater - Refreezing # free water content in snow - MaxFreeWater = self.DrySnow * self.WHC - self.FreeWater = self.FreeWater + self.SnowMelt + RainFall - InSoil = max( - self.FreeWater - MaxFreeWater, 0.0 - ) # abundant water in snow pack which goes into soil - self.FreeWater = self.FreeWater - InSoil - RainAndSnowmelt = RainFall + self.SnowMelt + self.SoilMoisture=self.SoilMoisture+NetInSoil + DirectRunoff=max(self.SoilMoisture-self.FieldCapacity,0.0) #if soil is filled to capacity: abundant water runs of directly + self.SoilMoisture=self.SoilMoisture-DirectRunoff + NetInSoil=NetInSoil-DirectRunoff #net water which infiltrates into soil - self.SnowCover = ifthenelse(self.DrySnow > 0, scalar(1), scalar(0)) - self.NrCell = areatotal(self.SnowCover, self.TopoId) + MaxSnowPack = 10000.0 + if self.MassWasting: + # Masswasting of snow + # 5.67 = tan 80 graden + SnowFluxFrac = min(0.5,self.Slope/5.67) * min(1.0,self.DrySnow/MaxSnowPack) + MaxFlux = SnowFluxFrac * self.DrySnow + self.DrySnow = accucapacitystate(self.TopoLdd,self.DrySnow, MaxFlux) + self.FreeWater = accucapacitystate(self.TopoLdd,self.FreeWater,SnowFluxFrac * self.FreeWater ) + else: + SnowFluxFrac = self.ZeroMap + MaxFlux= self.ZeroMap - # first part of precipitation is intercepted - # Interception=min(InSoil,self.ICF-self.InterceptionStorage)#: Interception in mm/timestep - # self.InterceptionStorage=self.InterceptionStorage+Interception #: Current interception storage - # NetInSoil=InSoil-Interception - NetInSoil = InSoil - self.SoilMoisture = self.SoilMoisture + NetInSoil - DirectRunoff = max( - self.SoilMoisture - self.FieldCapacity, 0.0 - ) # if soil is filled to capacity: abundant water runs of directly - self.SoilMoisture = self.SoilMoisture - DirectRunoff - NetInSoil = NetInSoil - DirectRunoff # net water which infiltrates into soil + #IntEvap=min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage + #self.InterceptionStorage=self.InterceptionStorage-IntEvap - MaxSnowPack = 10000.0 - if self.MassWasting: - # Masswasting of snow - # 5.67 = tan 80 graden - SnowFluxFrac = min(0.5, self.Slope / 5.67) * min( - 1.0, self.DrySnow / MaxSnowPack - ) - MaxFlux = SnowFluxFrac * self.DrySnow - self.DrySnow = accucapacitystate(self.TopoLdd, self.DrySnow, MaxFlux) - self.FreeWater = accucapacitystate( - self.TopoLdd, self.FreeWater, SnowFluxFrac * self.FreeWater - ) - else: - SnowFluxFrac = self.ZeroMap - MaxFlux = self.ZeroMap + # I nthe origal HBV code + #RestEvap = max(0.0,self.PotEvaporation-IntEvap) - # IntEvap=min(self.InterceptionStorage,self.PotEvaporation) #: Evaporation from interception storage - # self.InterceptionStorage=self.InterceptionStorage-IntEvap + self.SoilEvap=ifthenelse(self.SoilMoisture > self.Treshold,min(self.SoilMoisture,RestEvap),\ + min(self.SoilMoisture,min(RestEvap,self.PotEvaporation*(self.SoilMoisture/self.Treshold)))) + #: soil evapotranspiration + self.SoilMoisture=self.SoilMoisture-self.SoilEvap #evaporation from soil moisture storage - # I nthe origal HBV code - # RestEvap = max(0.0,self.PotEvaporation-IntEvap) - self.SoilEvap = ifthenelse( - self.SoilMoisture > self.Treshold, - min(self.SoilMoisture, RestEvap), - min( - self.SoilMoisture, - min( - RestEvap, self.PotEvaporation * (self.SoilMoisture / self.Treshold) - ), - ), - ) - #: soil evapotranspiration - self.SoilMoisture = ( - self.SoilMoisture - self.SoilEvap - ) # evaporation from soil moisture storage + self.ActEvap=self.IntEvap+self.SoilEvap #: Sum of evaporation components (IntEvap+SoilEvap) + self.HBVSeepage=((min(self.SoilMoisture/self.FieldCapacity,1))**self.BetaSeepage)*NetInSoil #runoff water from soil + self.SoilMoisture=self.SoilMoisture-self.HBVSeepage - self.ActEvap = ( - self.IntEvap + self.SoilEvap - ) #: Sum of evaporation components (IntEvap+SoilEvap) - self.HBVSeepage = ( - (min(self.SoilMoisture / self.FieldCapacity, 1)) ** self.BetaSeepage - ) * NetInSoil # runoff water from soil - self.SoilMoisture = self.SoilMoisture - self.HBVSeepage + Backtosoil=min(self.FieldCapacity-self.SoilMoisture,DirectRunoff) #correction for extremely wet periods: soil is filled to capacity + self.DirectRunoff=DirectRunoff-Backtosoil + self.SoilMoisture=self.SoilMoisture+Backtosoil + self.InUpperZone=self.DirectRunoff+self.HBVSeepage # total water available for runoff - Backtosoil = min( - self.FieldCapacity - self.SoilMoisture, DirectRunoff - ) # correction for extremely wet periods: soil is filled to capacity - self.DirectRunoff = DirectRunoff - Backtosoil - self.SoilMoisture = self.SoilMoisture + Backtosoil - self.InUpperZone = ( - self.DirectRunoff + self.HBVSeepage - ) # total water available for runoff + # Steps is always 1 at the moment + # calculations for Upper zone + self.UpperZoneStorage=self.UpperZoneStorage+self.InUpperZone #incoming water from soil + self.Percolation=min(self.PERC,self.UpperZoneStorage-self.InUpperZone/2) #Percolation + self.UpperZoneStorage=self.UpperZoneStorage-self.Percolation + self.CapFlux=self.Cflux*(((self.FieldCapacity-self.SoilMoisture)/self.FieldCapacity)) #: Capillary flux flowing back to soil + self.CapFlux=min(self.UpperZoneStorage,self.CapFlux) + self.CapFlux=min(self.FieldCapacity-self.SoilMoisture,self.CapFlux) + self.UpperZoneStorage=self.UpperZoneStorage-self.CapFlux + self.SoilMoisture=self.SoilMoisture+self.CapFlux - # Steps is always 1 at the moment - # calculations for Upper zone - self.UpperZoneStorage = ( - self.UpperZoneStorage + self.InUpperZone - ) # incoming water from soil - self.Percolation = min( - self.PERC, self.UpperZoneStorage - self.InUpperZone / 2 - ) # Percolation - self.UpperZoneStorage = self.UpperZoneStorage - self.Percolation - self.CapFlux = self.Cflux * ( - ((self.FieldCapacity - self.SoilMoisture) / self.FieldCapacity) - ) #: Capillary flux flowing back to soil - self.CapFlux = min(self.UpperZoneStorage, self.CapFlux) - self.CapFlux = min(self.FieldCapacity - self.SoilMoisture, self.CapFlux) - self.UpperZoneStorage = self.UpperZoneStorage - self.CapFlux - self.SoilMoisture = self.SoilMoisture + self.CapFlux + if not self.SetKquickFlow: + self.QuickFlow=min(ifthenelse(self.Percolation 0: + self.ReservoirVolume, self.Outflow, self.ResPercFull,\ + self.DemandRelease = simplereservoir(self.ReservoirVolume, self.SurfaceRunoff,\ + self.ResMaxVolume, self.ResTargetFullFrac, + self.ResMaxRelease, self.ResDemand, + self.ResTargetMinFrac, self.ReserVoirSimpleLocs, + timestepsecs=self.timestepsecs) + self.OutflowDwn = upstream(self.TopoLddOrg,cover(self.Outflow,scalar(0.0))) + self.Inflow = self.OutflowDwn + cover(self.Inflow,self.ZeroMap) + #else: + # self.Inflow= cover(self.Inflow,self.ZeroMap) - if self.nrresSimple > 0: - self.ReservoirVolume, self.Outflow, self.ResPercFull, self.DemandRelease = simplereservoir( - self.ReservoirVolume, - self.SurfaceRunoff, - self.ResMaxVolume, - self.ResTargetFullFrac, - self.ResMaxRelease, - self.ResDemand, - self.ResTargetMinFrac, - self.ReserVoirSimpleLocs, - timestepsecs=self.timestepsecs, - ) - self.OutflowDwn = upstream( - self.TopoLddOrg, cover(self.Outflow, scalar(0.0)) - ) - self.Inflow = self.OutflowDwn + cover(self.Inflow, self.ZeroMap) - # else: - # self.Inflow= cover(self.Inflow,self.ZeroMap) + elif self.nrresComplex > 0: + self.ReservoirWaterLevel, self.Outflow, self.ReservoirPrecipitation, self.ReservoirEvaporation,\ + self.ReservoirVolume = complexreservoir(self.ReservoirWaterLevel, self.ReserVoirComplexLocs, self.LinkedReservoirLocs, self.ResArea,\ + self.ResThreshold, self.ResStorFunc, self.ResOutflowFunc, self.sh, self.hq, self.Res_b, + self.Res_e, self.SurfaceRunoff,self.ReserVoirPrecip, self.ReserVoirPotEvap, + self.ReservoirComplexAreas, self.wf_supplyJulianDOY(), timestepsecs=self.timestepsecs) + self.OutflowDwn = upstream(self.TopoLddOrg,cover(self.Outflow,scalar(0.0))) + self.Inflow = self.OutflowDwn + cover(self.Inflow,self.ZeroMap) + else: + self.Inflow= cover(self.Inflow,self.ZeroMap) - elif self.nrresComplex > 0: - self.ReservoirWaterLevel, self.Outflow, self.ReservoirPrecipitation, self.ReservoirEvaporation, self.ReservoirVolume = complexreservoir( - self.ReservoirWaterLevel, - self.ReserVoirComplexLocs, - self.LinkedReservoirLocs, - self.ResArea, - self.ResThreshold, - self.ResStorFunc, - self.ResOutflowFunc, - self.sh, - self.hq, - self.Res_b, - self.Res_e, - self.SurfaceRunoff, - self.ReserVoirPrecip, - self.ReserVoirPotEvap, - self.ReservoirComplexAreas, - self.wf_supplyJulianDOY(), - timestepsecs=self.timestepsecs, - ) - self.OutflowDwn = upstream( - self.TopoLddOrg, cover(self.Outflow, scalar(0.0)) - ) - self.Inflow = self.OutflowDwn + cover(self.Inflow, self.ZeroMap) - else: - self.Inflow = cover(self.Inflow, self.ZeroMap) + self.QuickFlowCubic = (self.QuickFlow + self.RealQuickFlow) * self.ToCubic + self.BaseFlowCubic = self.BaseFlow * self.ToCubic - self.QuickFlowCubic = (self.QuickFlow + self.RealQuickFlow) * self.ToCubic - self.BaseFlowCubic = self.BaseFlow * self.ToCubic + self.SurfaceWaterSupply = ifthenelse (self.Inflow < 0.0 , max(-1.0 * self.Inwater,self.SurfaceRunoff), self.ZeroMap) + self.Inwater = self.Inwater + ifthenelse(self.SurfaceWaterSupply> 0, -1.0 * self.SurfaceWaterSupply,self.Inflow) - self.SurfaceWaterSupply = ifthenelse( - self.Inflow < 0.0, - max(-1.0 * self.Inwater, self.SurfaceRunoff), - self.ZeroMap, - ) - self.Inwater = self.Inwater + ifthenelse( - self.SurfaceWaterSupply > 0, -1.0 * self.SurfaceWaterSupply, self.Inflow - ) - ########################################################################## - # Runoff calculation via Kinematic wave ################################## - ########################################################################## - # per distance along stream - q = self.Inwater / self.DCL + self.ForecQ_qmec / self.DCL - self.OldSurfaceRunoff = self.SurfaceRunoff + ########################################################################## + # Runoff calculation via Kinematic wave ################################## + ########################################################################## + # per distance along stream + q=self.Inwater/self.DCL + self.ForecQ_qmec/self.DCL + self.OldSurfaceRunoff=self.SurfaceRunoff - self.SurfaceRunoff = kinematic( - self.TopoLdd, - self.SurfaceRunoff, - q, - self.Alpha, - self.Beta, - self.Tslice, - self.timestepsecs, - self.DCL, - ) # m3/s - self.SurfaceRunoffMM = ( - self.SurfaceRunoff * self.QMMConv - ) # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) + self.SurfaceRunoff = kinematic(self.TopoLdd, self.SurfaceRunoff,q,self.Alpha, self.Beta,self.Tslice,self.timestepsecs,self.DCL) # m3/s + self.SurfaceRunoffMM=self.SurfaceRunoff*self.QMMConv # SurfaceRunoffMM (mm) from SurfaceRunoff (m3/s) + + self.updateRunOff() + InflowKinWaveCell=upstream(self.TopoLdd,self.SurfaceRunoff) + self.MassBalKinWave = (self.KinWaveVolume - self.OldKinWaveVolume)/self.timestepsecs + InflowKinWaveCell + self.Inwater - self.SurfaceRunoff + Runoff=self.SurfaceRunoff + + # Updating + # -------- + # Assume a tss file with as many columns as outpulocs. Start updating for each non-missing value and start with the + # first column (nr 1). Assumes that outputloc and columns match! + + if self.updating: + QM = timeinputscalar(self.updateFile, self.UpdateMap) * self.QMMConv + + # Now update the state. Just add to the Ustore + # self.UStoreDepth = result + # No determine multiplication ratio for each gauge influence area. + # For missing gauges 1.0 is assumed (no change). + # UpDiff = areamaximum(QM, self.UpdateMap) - areamaximum(self.SurfaceRunoffMM, self.UpdateMap) + UpRatio = areamaximum(QM, self.UpdateMap)/areamaximum(self.SurfaceRunoffMM, self.UpdateMap) + + UpRatio = cover(areaaverage(UpRatio,self.TopoId),1.0) + # Now split between Soil and Kyn wave + self.UpRatioKyn = min(self.MaxUpdMult,max(self.MinUpdMult,(UpRatio - 1.0) * self.UpFrac + 1.0)) + UpRatioSoil = min(self.MaxUpdMult,max(self.MinUpdMult,(UpRatio - 1.0) * (1.0 - self.UpFrac) + 1.0)) + + # update/nudge self.UStoreDepth for the whole upstream area, + # not sure how much this helps or worsens things + UpdSoil = True + if UpdSoil: + toadd = (self.UpperZoneStorage * UpRatioSoil) - self.UpperZoneStorage + self.UpperZoneStorage = self.UpperZoneStorage + toadd + + # Update the kinematic wave reservoir up to a maximum upstream distance + # TODO: add (much smaller) downstream updating also? + MM = (1.0 - self.UpRatioKyn)/self.UpdMaxDist + self.UpRatioKyn = MM * self.DistToUpdPt + self.UpRatioKyn + self.updateRunOff() InflowKinWaveCell = upstream(self.TopoLdd, self.SurfaceRunoff) self.MassBalKinWave = (