Index: doc/_images/wflow_routing.png =================================================================== diff -u Binary files differ Index: doc/wflow_routing.rst =================================================================== diff -u -r609e9a8b2b084b2a8b6ec809e3e0479765a51b03 -r34c2b30d41c043c3a36fb6b265acea8484cd201e --- doc/wflow_routing.rst (.../wflow_routing.rst) (revision 609e9a8b2b084b2a8b6ec809e3e0479765a51b03) +++ doc/wflow_routing.rst (.../wflow_routing.rst) (revision 34c2b30d41c043c3a36fb6b265acea8484cd201e) @@ -8,9 +8,42 @@ Introduction ------------ -A +The wflow routing module uses the pcraster kinematic wave to route water over a DEM. By adding a bankflull level +and a floodplainwidth to the configuration the model can also inclue estimate flow over a floodplain. +Method +====== + """ + Updates the kinematic wave reservoir. Should be run after updates to Q + + WL = A * pow(Q,B)/Bw + WL/A * Bw = pow(Q,B) + Q = pow(WL/A * Bw,1/B) + """ + # TODO: This is not correct, need to subtract Channel runoff + self.Qbankfull = pow(self.bankFull/self.AlphaCh * self.Bw,1.0/self.Beta) + self.Qchannel = min(self.SurfaceRunoff,self.Qbankfull) + self.floodcells = boolean(ifthenelse(self.WaterLevelCH > self.bankFull, boolean(1), boolean(0))) + self.Qfloodplain = max(0.0,self.SurfaceRunoff - self.Qbankfull) + + self.WaterLevelCH = self.AlphaCh * pow(self.Qchannel, self.Beta) / (self.Bw) + self.WaterLevelFP = self.AlphaFP * pow(self.Qfloodplain, self.Beta) / (self.Bw + self.Pfp) + #self.WaterLevelFP = max(0.0,self.WaterLevelCH - self.bankFull) + #self.Alpha * pow(self.SurfaceRunoff-self.Qchannel, self.Beta) / (self.Bw + self.Pch) + self.bankFull) + self.Qtot = pow(self.WaterLevelCH/self.AlphaCh * self.Bw,1.0/self.Beta) + pow(self.WaterLevelFP/self.AlphaFP * self.Pfp,1.0/self.Beta) + + # wetted perimeter (m) + self.Pch = self.wetPerimiterCH(self.WaterLevelCH,self.Bw) + self.Pfp = self.wetPerimiterFP(self.WaterLevelFP,self.floodPlainWidth) + # Alpha + self.AlphaFP = self.AlpTerm * pow(self.Pfp, self.AlpPow) + self.AlphaCh = self.AlpTerm * pow(self.Pch, self.AlpPow) + self.Alpha = self.AlpTerm * pow(self.Pch + self.Pfp, self.AlpPow) + self.OldKinWaveVolume = self.KinWaveVolume + self.KinWaveVolume = self.WaterLevelCH * self.Bw * self.DCL + + Dependencies ------------ T @@ -45,6 +78,7 @@ + wflow_routing module documentation ---------------------------------- Index: examples/wflow_rhine_sbm/wflow_routing.ini =================================================================== diff -u -re574688d9f85c8893b24b7e53e676bb2cd4f0763 -r34c2b30d41c043c3a36fb6b265acea8484cd201e --- examples/wflow_rhine_sbm/wflow_routing.ini (.../wflow_routing.ini) (revision e574688d9f85c8893b24b7e53e676bb2cd4f0763) +++ examples/wflow_rhine_sbm/wflow_routing.ini (.../wflow_routing.ini) (revision 34c2b30d41c043c3a36fb6b265acea8484cd201e) @@ -1,10 +1,8 @@ -[API] -IW=0,4 [inputmapstacks] -Inwater = /inmaps/IW +Inwater = /run_default/outmaps/IW # this is optional inflow (positive) or outflow (negative) of the kin-wave reservoir in m^3/sec Inflow = /inmaps/IF @@ -47,5 +45,18 @@ [outputmaps] self.SurfaceRunoff=_run -#self.SubCellFrac=scf - +self.Qfloodplain=_qfp +self.Qchannel=_qch +self.Qbankfull=_qbn +self.WaterLevelFP=_levfp +self.WaterLevelCH=_levch +self.InwaterMM=_IW +self.Qbankfull=_qbnk +self.floodcells=fcel +self.Qtot=QQQ +self.Pch =ch +self.Pfp = fp +self.Alpha = al +self.AlphaCh = alch +self.AlphaFP = alfp +self.Ncombined = nc Index: examples/wflow_rhine_sbm/wflow_sbm.ini =================================================================== diff -u -r0c107ad8319af1405fc28d9a192a0d6aced12bdb -r34c2b30d41c043c3a36fb6b265acea8484cd201e --- examples/wflow_rhine_sbm/wflow_sbm.ini (.../wflow_sbm.ini) (revision 0c107ad8319af1405fc28d9a192a0d6aced12bdb) +++ examples/wflow_rhine_sbm/wflow_sbm.ini (.../wflow_sbm.ini) (revision 34c2b30d41c043c3a36fb6b265acea8484cd201e) @@ -67,8 +67,8 @@ # netcdfoutput requires also outputformat = 1 (default) and additionally the name of the file #netcdfoutput = outmaps.nc -netcdfstaticoutput = outsum.nc -netcdfwritebuffer=100 +#netcdfstaticoutput = outsum.nc +#netcdfwritebuffer=100 [layout] # if set to zero the cell-size is given in lat/long (the default) @@ -80,14 +80,15 @@ [outputmaps] self.zi=zi self.SurfaceRunoff=run +self.WaterLevel=lev #self.SubCellFrac=scf -#self.Inwater=inw +#self.Inwater=IW #self.DistToUpdPt=dist #self.SnowMelt=sno self.FirstZoneDepth=fzd #self.TopoLdd=ldd -self.InwaterMM=iwm +self.InwaterMM=IW #self.watbal=wat self.PotenEvap=PET self.Precipitation=P Index: wflow-py/UnitTests/wflow_sceleton/outstate/TSoil.map =================================================================== diff -u -rec91b46ffc9c4e6e2c442f66cd720176ea9632f9 -r34c2b30d41c043c3a36fb6b265acea8484cd201e Binary files differ Index: wflow-py/wflow/wflow_hbv.py =================================================================== diff -u -r43965a4158a9e88bc6aae38bb1aaccf2e44aba45 -r34c2b30d41c043c3a36fb6b265acea8484cd201e --- wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision 43965a4158a9e88bc6aae38bb1aaccf2e44aba45) +++ wflow-py/wflow/wflow_hbv.py (.../wflow_hbv.py) (revision 34c2b30d41c043c3a36fb6b265acea8484cd201e) @@ -98,7 +98,7 @@ wflow = "wflow_hbv: " -wflowVersion = "$Revision: 904 $ $Date: 2014-01-13 15:39:24 +0100 (Mon, 13 Jan 2014) $" +wflowVersion = "_" """revision of the model""" Index: wflow-py/wflow/wflow_routing.py =================================================================== diff -u -r43965a4158a9e88bc6aae38bb1aaccf2e44aba45 -r34c2b30d41c043c3a36fb6b265acea8484cd201e --- wflow-py/wflow/wflow_routing.py (.../wflow_routing.py) (revision 43965a4158a9e88bc6aae38bb1aaccf2e44aba45) +++ wflow-py/wflow/wflow_routing.py (.../wflow_routing.py) (revision 34c2b30d41c043c3a36fb6b265acea8484cd201e) @@ -93,6 +93,7 @@ """ + def __init__(self, cloneMap, Dir, RunDir, configfile): DynamicModel.__init__(self) @@ -104,6 +105,41 @@ self.configfile = configfile self.SaveDir = os.path.join(self.Dir,self.runId) + + def wetPerimiterFP(self,Waterlevel, floodplainwidth,threshold=0.0,sharpness=0.5): + """ + + :param Waterlevel: + :param bottomwidth: + :param bankfull: + :param floodplainwidth: + :return P: wetted perimiter + """ + a = threshold + b = 1.0 + c = sharpness # not very sharp + + floodplainfact = max(0.001,sCurve(Waterlevel, a=a, c=c,b=b) -0.5) + floodplainperimiter = min(1.0,2.0 * floodplainfact) * floodplainwidth + + return floodplainperimiter + + + def wetPerimiterCH(self,Waterlevel,channelWidth): + """ + + :param Waterlevel: + :param bottomwidth: + :param bankfull: + :param floodplainwidth: + :return P: wetted perimiter + """ + + channelperimiter = 2.0 * Waterlevel + channelWidth + + return channelperimiter + + def _initAPIVars(self): """ Sets vars in the API that are forcing variables to the model @@ -115,15 +151,40 @@ def updateRunOff(self): """ - Updates the kinematic wave reservoir. Should be run after updates to Q - """ - self.WaterLevel = (self.Alpha * pow(self.SurfaceRunoff, self.Beta)) / self.Bw + Updates the kinematic wave reservoir. Should be run after updates to Q + + WL = A * pow(Q,B)/Bw + WL/A * Bw = pow(Q,B) + Q = pow(WL/A * Bw,1/B) + """ + + self.Qbankfull = pow(self.bankFull/self.AlphaCh * self.Bw,1.0/self.Beta) + self.Qchannel = min(self.SurfaceRunoff,self.Qbankfull) + self.floodcells = boolean(ifthenelse(self.WaterLevelCH > self.bankFull, boolean(1), boolean(0))) + self.Qfloodplain = max(0.0,self.SurfaceRunoff - self.Qbankfull) + + self.WaterLevelCH = self.AlphaCh * pow(self.Qchannel, self.Beta) / (self.Bw) + self.WaterLevelFP = ifthenelse(self.River,self.AlphaFP * pow(self.Qfloodplain, self.Beta) / (self.Bw + self.Pfp),0.0) + self.WaterLevel = self.WaterLevelCH + self.WaterLevelFP + + # Dtermien Qtot as a check + self.Qtot = pow(self.WaterLevelCH/self.AlphaCh * self.Bw,1.0/self.Beta) + pow(self.WaterLevelFP/self.AlphaFP * self.Pfp,1.0/self.Beta) # wetted perimeter (m) - P = self.Bw + (2 * self.WaterLevel) + self.Pch = self.wetPerimiterCH(self.WaterLevelCH,self.Bw) + self.Pfp = ifthenelse(self.River,self.wetPerimiterFP(self.WaterLevelFP,self.floodPlainWidth + self.Bw,sharpness=self.floodPlainDist),0.0) # Alpha - self.Alpha = self.AlpTerm * pow(P, self.AlpPow) + + self.WetPComb = self.Pch + self.Pfp + self.Ncombined = (self.Pch/self.WetPComb*self.N**1.5 + self.Pfp/self.WetPComb*self.NFloodPlain**1.5)**(2./3.) + + self.AlpTermFP = pow((self.NFloodPlain / (sqrt(self.SlopeDCL))), self.Beta) + self.AlpTermComb = pow((self.Ncombined / (sqrt(self.SlopeDCL))), self.Beta) + + self.AlphaFP = self.AlpTermFP * pow(self.Pfp, self.AlpPow) + self.AlphaCh = self.AlpTerm * pow(self.Pch, self.AlpPow) + self.Alpha = ifthenelse(self.River,self.AlpTermComb * pow(self.Pch + self.Pfp, self.AlpPow),self.AlphaCh) self.OldKinWaveVolume = self.KinWaveVolume - self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL + self.KinWaveVolume = (self.WaterLevelCH * self.Bw * self.DCL) + (self.WaterLevelFP * self.Pfp * self.DCL) def stateVariables(self): """ @@ -135,7 +196,7 @@ :var self.SurfaceRunoff: Surface runoff in the kin-wave resrvoir [m^3/s] :var self.WaterLevel: Water level in the kin-wave resrvoir [m] """ - states = ['SurfaceRunoff', 'WaterLevel'] + states = ['SurfaceRunoff', 'WaterLevelCH','WaterLevelFP'] return states @@ -215,6 +276,10 @@ wflow_gauges = configget(self.config, "model", "wflow_gauges", "staticmaps/wflow_gauges.map") wflow_inflow = configget(self.config, "model", "wflow_inflow", "staticmaps/wflow_inflow.map") wflow_riverwidth = configget(self.config, "model", "wflow_riverwidth", "staticmaps/wflow_riverwidth.map") + wflow_floodplainwidth = configget(self.config, "model", "wflow_floodplainwidth", "staticmaps/wflow_floodplainwidth.map") + wflow_bankfulldepth = configget(self.config, "model", "wflow_bankfulldepth", "staticmaps/wflow_bankfulldepth.map") + wflow_floodplaindist = configget(self.config, "model", "wflow_bankfulldepth", "staticmaps/wflow_floodplaindist.map") + wflow_landuse = configget(self.config, "model", "wflow_landuse", "staticmaps/wflow_landuse.map") wflow_soil = configget(self.config, "model", "wflow_soil", "staticmaps/wflow_soil.map") @@ -241,26 +306,27 @@ self.OutputLoc = readmap(os.path.join(self.Dir,wflow_gauges)) # location of output gauge(s) self.InflowLoc = pcrut.readmapSave(os.path.join(self.Dir,wflow_inflow), 0.0) # location abstractions/inflows. self.RiverWidth = pcrut.readmapSave(os.path.join(self.Dir,wflow_riverwidth), 0.0) + self.bankFull = pcrut.readmapSave(os.path.join(self.Dir,wflow_bankfulldepth), 16.0) + self.floodPlainWidth = pcrut.readmapSave(os.path.join(self.Dir,wflow_floodplainwidth), 8000.0) + self.floodPlainDist = pcrut.readmapSave(os.path.join(self.Dir,wflow_floodplaindist), 0.5) + self.OutputId = readmap(os.path.join(self.Dir,wflow_subcatch)) # location of subcatchment self.ZeroMap = 0.0 * scalar(subcatch) #map with only zero's - self.IW_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Inwater", - "/inmaps/IW") # timeseries for specific runoff - self.Inflow_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Inflow", - "/inmaps/IF") # timeseries for rainfall "/inmaps/IF" # in/outflow locations (abstractions) - self.Latitude = ycoordinate(boolean(self.Altitude)) self.Longitude = xcoordinate(boolean(self.Altitude)) self.logger.info("Linking parameters to landuse, catchment and soil...") + self.wf_updateparameters() self.Beta = scalar(0.6) # For sheetflow self.N = self.readtblDefault(self.Dir + "/" + self.intbl + "/N.tbl", self.LandUse, subcatch, self.Soil, 0.072) # Manning overland flow self.NRiver = self.readtblDefault(self.Dir + "/" + self.intbl + "/N_River.tbl", self.LandUse, subcatch, self.Soil, 0.036) # Manning river - + self.NFloodPlain = self.readtblDefault(self.Dir + "/" + self.intbl + "/N_FloodPlain.tbl", self.LandUse, subcatch, + self.Soil, self.NRiver * 2.0) # Manning river self.xl, self.yl, self.reallength = pcrut.detRealCellLength(self.ZeroMap, sizeinmetres) self.Slope = slope(self.Altitude) #self.Slope=ifthen(boolean(self.TopoId),max(0.001,self.Slope*celllength()/self.reallength)) @@ -377,23 +443,61 @@ return lst + + def parameters(self): + """ + Define all model parameters here that the framework should handle for the model + See wf_updateparameters and the parameters section of the ini file + If you use this make sure to all wf_updateparameters at the start of the dynamic section + and at the start/end of the initial section + """ + modelparameters = [] + + #Static model parameters e.g. + #modelparameters.append(self.ParamType(name="RunoffGeneratingGWPerc",stack="intbl/RunoffGeneratingGWPerc.tbl",type="static",default=0.1)) + # 3: Input time series ################################################### + self.IW_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Inwater", + "/inmaps/IW") # timeseries for specific runoff + + self.Inflow_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Inflow", + "/inmaps/IF") # timeseries for rainfall "/inmaps/IF" # in/outflow locations (abstractions) + + # Meteo and other forcing + modelparameters.append(self.ParamType(name="InwaterForcing",stack=self.IW_mapstack ,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) + modelparameters.append(self.ParamType(name="Inflow",stack=self.Inflow_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[])) + + + return modelparameters + def resume(self): if self.reinit == 1: self.logger.info("Setting initial conditions to default") - self.WaterLevel = self.ZeroMap + self.WaterLevelCH = self.ZeroMap + self.WaterLevelFP = self.ZeroMap self.SurfaceRunoff = self.ZeroMap + self.WaterLevel = self.WaterLevelCH + self.WaterLevelFP else: self.logger.info("Setting initial conditions from state files") self.wf_resume(os.path.join(self.Dir,"instate")) - P = self.Bw + (2.0 * self.WaterLevel) - self.Alpha = self.AlpTerm * pow(P, self.AlpPow) + self.Pch = self.wetPerimiterCH(self.WaterLevelCH,self.Bw) + self.Pfp = ifthenelse(self.River,self.wetPerimiterFP(self.WaterLevelFP,self.floodPlainWidth + self.Bw,sharpness=self.floodPlainDist),0.0) + self.WetPComb = self.Pch + self.Pfp + self.Ncombined = (self.Pch/self.WetPComb*self.N**1.5 + self.Pfp/self.WetPComb*self.NFloodPlain**1.5)**(2./3.) + + + self.AlpTermFP = pow((self.NFloodPlain / (sqrt(self.SlopeDCL))), self.Beta) + self.AlpTermComb = pow((self.Ncombined / (sqrt(self.SlopeDCL))), self.Beta) + + self.AlphaFP = self.AlpTermFP * pow(self.Pfp, self.AlpPow) + self.AlphaCh = self.AlpTerm * pow(self.Pch, self.AlpPow) + self.Alpha = ifthenelse(self.River,self.AlpTermComb * pow(self.Pch + self.Pfp, self.AlpPow),self.AlphaCh) self.OldSurfaceRunoff = self.SurfaceRunoff self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv # Determine initial kinematic wave volume - self.KinWaveVolume = self.WaterLevel * self.Bw * self.DCL + self.KinWaveVolume = (self.WaterLevelCH * self.Bw * self.DCL) + (self.WaterLevelFP * self.Pfp * self.DCL) self.OldKinWaveVolume = self.KinWaveVolume self.SurfaceRunoffMM = self.SurfaceRunoff * self.QMMConv @@ -421,16 +525,8 @@ self.logger.debug("Step: " + str(int(self.currentStep)) + "/" + str(int(self._d_nrTimeSteps))) self.thestep = self.thestep + 1 - self.InwaterForcing = cover(self.wf_readmap(self.IW_mapstack, 1.0), scalar(0.0)) + self.wf_updateparameters() - if self.thestep > 28: - self.InwaterForcing = cover(0.0) - #self.Inflow=cover(self.wf_readmap(self.Inflow),0) - if (os.path.exists(self.caseName + self.inflowTss)): - self.Inflow = cover(timeinputscalar(self.caseName + self.inflowTss, nominal(self.InflowLoc)), 0) - else: - self.Inflow = cover(self.wf_readmap(self.Inflow_mapstack, 0.0,verbose=False),0) - # The MAx here may lead to watbal error. Howevere, if inwaterMMM becomes < 0, the kinematic wave becomes very slow...... self.InwaterMM = max(0.0,self.InwaterForcing) self.Inwater = self.InwaterMM * self.ToCubic # m3/s Index: wflow-py/wflow/wflow_sbm.py =================================================================== diff -u -r43965a4158a9e88bc6aae38bb1aaccf2e44aba45 -r34c2b30d41c043c3a36fb6b265acea8484cd201e --- wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 43965a4158a9e88bc6aae38bb1aaccf2e44aba45) +++ wflow-py/wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 34c2b30d41c043c3a36fb6b265acea8484cd201e) @@ -87,6 +87,7 @@ -l: loglevel (most be one of DEBUG, WARNING, ERROR) + """ #TODO: add Et reduction in unsat zone based on deficit @@ -403,11 +404,8 @@ self.SSSF = False setglobaloption("unittrue") - self.precipTss = "/intss/P.tss" - self.evapTss = "/intss/PET.tss" - self.tempTss = "/intss/T.tss" - self.inflowTss = "/intss/Inflow.tss" + self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps")