Index: wflow/sphy/advanced_routing.py =================================================================== diff -u -r5bc2645dde878c05902af9d9e373036e8d9908b2 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/sphy/advanced_routing.py (.../advanced_routing.py) (revision 5bc2645dde878c05902af9d9e373036e8d9908b2) +++ wflow/sphy/advanced_routing.py (.../advanced_routing.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -1,11 +1,13 @@ -#-Function to rout the specific runoff +# -Function to rout the specific runoff def ROUT(self, pcr, rvolume, qold, qout, sres): # Calculate the discharge Q (m3/d) Q = pcr.accufractionflux(self.FlowDir, rvolume, self.QFRAC) # Re-calculate Q, based on qold en kx, and assign Qout for cells being lake/reservoir - Q = pcr.ifthenelse(self.QFRAC==0, qout, (1-self.kx) * Q + (qold*24*3600) * self.kx) + Q = pcr.ifthenelse( + self.QFRAC == 0, qout, (1 - self.kx) * Q + (qold * 24 * 3600) * self.kx + ) # Only calculate inflow for lake/reservoir cells - Qin = pcr.ifthenelse(self.QFRAC==0, pcr.upstream(self.FlowDir, Q), 0) + Qin = pcr.ifthenelse(self.QFRAC == 0, pcr.upstream(self.FlowDir, Q), 0) sres = sres - qout + Qin - Q = Q / (24 * 3600) #-only convert Q to m3/s + Q = Q / (24 * 3600) # -only convert Q to m3/s return sres, Q, Qin Index: wflow/sphy/dynamic_veg.py =================================================================== diff -u -r5bc2645dde878c05902af9d9e373036e8d9908b2 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/sphy/dynamic_veg.py (.../dynamic_veg.py) (revision 5bc2645dde878c05902af9d9e373036e8d9908b2) +++ wflow/sphy/dynamic_veg.py (.../dynamic_veg.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -1,18 +1,25 @@ -#-Function that returns crop factor (Kc) and maximum storage (Smax) -def Veg_function(pcr, ndvi, fpar_max, fpar_min, lai_max, ndvi_min, ndvi_max, kc_min, kc_max): - SR = (1 + ndvi)/(1 - ndvi) - SR_max = (1 + ndvi_max)/(1 - ndvi_max) - SR_min = (1 + ndvi_min)/(1 - ndvi_min) - FPAR = pcr.min((((SR - SR_min) * (fpar_max - fpar_min))/ (SR_max - SR_min)) + 0.001, 0.95) - LAI = lai_max * pcr.log10(1-FPAR)/pcr.log10(1-fpar_max) - Smax = 0.935 + 0.498*LAI - 0.00575*(LAI**2) - Kc = kc_min + (kc_max - kc_min) * pcr.max(pcr.min((ndvi - ndvi_min)/(ndvi_max - ndvi_min), 1), 0) - return Kc, Smax +# -Function that returns crop factor (Kc) and maximum storage (Smax) +def Veg_function( + pcr, ndvi, fpar_max, fpar_min, lai_max, ndvi_min, ndvi_max, kc_min, kc_max +): + SR = (1 + ndvi) / (1 - ndvi) + SR_max = (1 + ndvi_max) / (1 - ndvi_max) + SR_min = (1 + ndvi_min) / (1 - ndvi_min) + FPAR = pcr.min( + (((SR - SR_min) * (fpar_max - fpar_min)) / (SR_max - SR_min)) + 0.001, 0.95 + ) + LAI = lai_max * pcr.log10(1 - FPAR) / pcr.log10(1 - fpar_max) + Smax = 0.935 + 0.498 * LAI - 0.00575 * (LAI ** 2) + Kc = kc_min + (kc_max - kc_min) * pcr.max( + pcr.min((ndvi - ndvi_min) / (ndvi_max - ndvi_min), 1), 0 + ) + return Kc, Smax -#-Function that returns the interception, precipitation throughfall, and remaining storage + +# -Function that returns the interception, precipitation throughfall, and remaining storage def Inter_function(pcr, S, Smax, Etr): PreT = pcr.max(0, S - Smax) S = S - PreT Int = pcr.min(1.5 * Etr, S) S = S - Int - return Int, PreT, S + return Int, PreT, S Index: wflow/sphy/lakes.py =================================================================== diff -u -r5bc2645dde878c05902af9d9e373036e8d9908b2 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/sphy/lakes.py (.../lakes.py) (revision 5bc2645dde878c05902af9d9e373036e8d9908b2) +++ wflow/sphy/lakes.py (.../lakes.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -1,52 +1,134 @@ -#-Function that updates the lake storage and lake level given a measured lake level. If no lake +# -Function that updates the lake storage and lake level given a measured lake level. If no lake # level is measured, then the actual storage is not updated with a measured level. The function # returns the updated storage and lake level def UpdateLakeHStore(self, pcr, pcrm): - #-buffer actual storage + # -buffer actual storage OldStorage = self.StorRES - #-Check if measured lake levels area available + # -Check if measured lake levels area available try: LakeLevel = pcr.readmap(pcrm.generateNameT(self.LLevel, self.counter)) Level = True except: Level = False if Level: - #-update lake storage according to measured water level - self.StorRES = pcr.ifthenelse(self.UpdateLakeLevel, pcr.ifthenelse(pcr.defined(LakeLevel), pcr.ifthenelse(self.LakeSH_Func==1,\ - self.LakeSH_exp_a * pcr.exp(self.LakeSH_exp_b * LakeLevel), pcr.ifthenelse(self.LakeSH_Func==2, self.LakeSH_pol_a1 \ - * LakeLevel + self.LakeSH_pol_b, pcr.ifthenelse(self.LakeSH_Func==3, (self.LakeSH_pol_a2 * LakeLevel**2) + \ - self.LakeSH_pol_a1 * LakeLevel + self.LakeSH_pol_b, (self.LakeSH_pol_a3 * LakeLevel**3) + (self.LakeSH_pol_a2 \ - * LakeLevel**2) + (self.LakeSH_pol_a1 * LakeLevel + self.LakeSH_pol_b)))), self.StorRES), self.StorRES) + # -update lake storage according to measured water level + self.StorRES = pcr.ifthenelse( + self.UpdateLakeLevel, + pcr.ifthenelse( + pcr.defined(LakeLevel), + pcr.ifthenelse( + self.LakeSH_Func == 1, + self.LakeSH_exp_a * pcr.exp(self.LakeSH_exp_b * LakeLevel), + pcr.ifthenelse( + self.LakeSH_Func == 2, + self.LakeSH_pol_a1 * LakeLevel + self.LakeSH_pol_b, + pcr.ifthenelse( + self.LakeSH_Func == 3, + (self.LakeSH_pol_a2 * LakeLevel ** 2) + + self.LakeSH_pol_a1 * LakeLevel + + self.LakeSH_pol_b, + (self.LakeSH_pol_a3 * LakeLevel ** 3) + + (self.LakeSH_pol_a2 * LakeLevel ** 2) + + (self.LakeSH_pol_a1 * LakeLevel + self.LakeSH_pol_b), + ), + ), + ), + self.StorRES, + ), + self.StorRES, + ) # prevent storage becoming negative for whatever reason self.StorRES = pcr.max(self.StorRES, 0) - #-Update the lake level based on the storage for lakes where no levels are measured - LakeLevel = pcr.ifthenelse(self.UpdateLakeLevel, pcr.ifthenelse(pcr.defined(LakeLevel), LakeLevel, \ - pcr.ifthenelse(self.LakeHS_Func==1, self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), pcr.ifthenelse(self.LakeHS_Func==2, self.LakeHS_pol_a1 * \ - self.StorRES + self.LakeHS_pol_b, pcr.ifthenelse(self.LakeHS_Func==3, (self.LakeHS_pol_a2 * self.StorRES**2) + \ - self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, (self.LakeHS_pol_a3 * self.StorRES**3) + (self.LakeHS_pol_a2 *\ - self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b)))), pcr.ifthenelse(self.LakeHS_Func==1, \ - self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), pcr.ifthenelse(self.LakeHS_Func==2, self.LakeHS_pol_a1 * \ - self.StorRES + self.LakeHS_pol_b, pcr.ifthenelse(self.LakeHS_Func==3, (self.LakeHS_pol_a2 * self.StorRES**2) + \ - self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, (self.LakeHS_pol_a3 * self.StorRES**3) + (self.LakeHS_pol_a2 *\ - self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b)))) + # -Update the lake level based on the storage for lakes where no levels are measured + LakeLevel = pcr.ifthenelse( + self.UpdateLakeLevel, + pcr.ifthenelse( + pcr.defined(LakeLevel), + LakeLevel, + pcr.ifthenelse( + self.LakeHS_Func == 1, + self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), + pcr.ifthenelse( + self.LakeHS_Func == 2, + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, + pcr.ifthenelse( + self.LakeHS_Func == 3, + (self.LakeHS_pol_a2 * self.StorRES ** 2) + + self.LakeHS_pol_a1 * self.StorRES + + self.LakeHS_pol_b, + (self.LakeHS_pol_a3 * self.StorRES ** 3) + + (self.LakeHS_pol_a2 * self.StorRES ** 2) + + self.LakeHS_pol_a1 * self.StorRES + + self.LakeHS_pol_b, + ), + ), + ), + ), + pcr.ifthenelse( + self.LakeHS_Func == 1, + self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), + pcr.ifthenelse( + self.LakeHS_Func == 2, + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, + pcr.ifthenelse( + self.LakeHS_Func == 3, + (self.LakeHS_pol_a2 * self.StorRES ** 2) + + self.LakeHS_pol_a1 * self.StorRES + + self.LakeHS_pol_b, + (self.LakeHS_pol_a3 * self.StorRES ** 3) + + (self.LakeHS_pol_a2 * self.StorRES ** 2) + + self.LakeHS_pol_a1 * self.StorRES + + self.LakeHS_pol_b, + ), + ), + ), + ) else: # if no lake level map is available, then calculate the h based on storages - LakeLevel = pcr.ifthenelse(self.LakeHS_Func==1, self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), \ - pcr.ifthenelse(self.LakeHS_Func==2, self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, pcr.ifthenelse(\ - self.LakeHS_Func==3, (self.LakeHS_pol_a2 * self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b,\ - (self.LakeHS_pol_a3 * self.StorRES**3) + (self.LakeHS_pol_a2 * self.StorRES**2) + self.LakeHS_pol_a1 * self.StorRES +\ - self.LakeHS_pol_b))) + LakeLevel = pcr.ifthenelse( + self.LakeHS_Func == 1, + self.LakeHS_exp_a * pcr.exp(self.LakeHS_exp_b * self.StorRES), + pcr.ifthenelse( + self.LakeHS_Func == 2, + self.LakeHS_pol_a1 * self.StorRES + self.LakeHS_pol_b, + pcr.ifthenelse( + self.LakeHS_Func == 3, + (self.LakeHS_pol_a2 * self.StorRES ** 2) + + self.LakeHS_pol_a1 * self.StorRES + + self.LakeHS_pol_b, + (self.LakeHS_pol_a3 * self.StorRES ** 3) + + (self.LakeHS_pol_a2 * self.StorRES ** 2) + + self.LakeHS_pol_a1 * self.StorRES + + self.LakeHS_pol_b, + ), + ), + ) self.StorRES = pcr.ifthenelse(self.LakeID != 0, self.StorRES, OldStorage) return LakeLevel, self.StorRES -#-function that calculates the fraction of lake storage that is available for routing, and the lake outflow + +# -function that calculates the fraction of lake storage that is available for routing, and the lake outflow def QLake(self, pcr, LakeLevel): - Qout = pcr.ifthenelse(self.LakeQH_Func==1, self.LakeQH_exp_a * pcr.exp(self.LakeQH_exp_b * LakeLevel), pcr.ifthenelse(\ - self.LakeQH_Func==2, self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b, pcr.ifthenelse(self.LakeQH_Func==3, \ - (self.LakeQH_pol_a2 * LakeLevel**2) + self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b, (self.LakeQH_pol_a3 * \ - LakeLevel**3) + (self.LakeQH_pol_a2 * LakeLevel**2) + self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b))) + Qout = pcr.ifthenelse( + self.LakeQH_Func == 1, + self.LakeQH_exp_a * pcr.exp(self.LakeQH_exp_b * LakeLevel), + pcr.ifthenelse( + self.LakeQH_Func == 2, + self.LakeQH_pol_a1 * LakeLevel + self.LakeQH_pol_b, + pcr.ifthenelse( + self.LakeQH_Func == 3, + (self.LakeQH_pol_a2 * LakeLevel ** 2) + + self.LakeQH_pol_a1 * LakeLevel + + self.LakeQH_pol_b, + (self.LakeQH_pol_a3 * LakeLevel ** 3) + + (self.LakeQH_pol_a2 * LakeLevel ** 2) + + self.LakeQH_pol_a1 * LakeLevel + + self.LakeQH_pol_b, + ), + ), + ) Qout = pcr.max(0, Qout) - Qout = Qout * 3600 * 24 #-convert to m3/d - Qout = pcr.cover(Qout, 0) #-for non-lake cells, Qout is zero + Qout = Qout * 3600 * 24 # -convert to m3/d + Qout = pcr.cover(Qout, 0) # -for non-lake cells, Qout is zero return Qout Index: wflow/sphy/reservoirs.py =================================================================== diff -u -r5bc2645dde878c05902af9d9e373036e8d9908b2 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/sphy/reservoirs.py (.../reservoirs.py) (revision 5bc2645dde878c05902af9d9e373036e8d9908b2) +++ wflow/sphy/reservoirs.py (.../reservoirs.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -1,29 +1,59 @@ -#-Advanced reservoir +# -Advanced reservoir def QAdv(self, pcr): DayNo = self.timecalc.julian(self)[0] - #-determine if it is flood or dry season - S1 = pcr.ifthenelse(self.ResFlStart < self.ResFlEnd, pcr.ifthenelse(DayNo>=self.ResFlStart, pcr.ifthenelse(DayNo<=self.ResFlEnd, pcr.boolean(1), pcr.boolean(0)), pcr.boolean(0)),\ - pcr.ifthenelse(DayNo>=self.ResFlEnd, pcr.ifthenelse(DayNo>=self.ResFlStart, pcr.boolean(1), pcr.boolean(0)), pcr.ifthenelse(DayNo<=self.ResFlEnd, \ - pcr.ifthenelse(DayNo<=self.ResFlStart, pcr.boolean(1), pcr.boolean(0)), pcr.boolean(0)))) + # -determine if it is flood or dry season + S1 = pcr.ifthenelse( + self.ResFlStart < self.ResFlEnd, + pcr.ifthenelse( + DayNo >= self.ResFlStart, + pcr.ifthenelse(DayNo <= self.ResFlEnd, pcr.boolean(1), pcr.boolean(0)), + pcr.boolean(0), + ), + pcr.ifthenelse( + DayNo >= self.ResFlEnd, + pcr.ifthenelse(DayNo >= self.ResFlStart, pcr.boolean(1), pcr.boolean(0)), + pcr.ifthenelse( + DayNo <= self.ResFlEnd, + pcr.ifthenelse( + DayNo <= self.ResFlStart, pcr.boolean(1), pcr.boolean(0) + ), + pcr.boolean(0), + ), + ), + ) S_avail = pcr.max(self.StorRES - self.ResPVOL, 0) - Q = pcr.max(pcr.ifthenelse(S1, self.ResMaxFl * S_avail / (self.ResEVOL - self.ResPVOL),\ - self.ResDemFl * S_avail / (self.ResEVOL - self.ResPVOL)), self.StorRES - self.ResEVOL) + Q = pcr.max( + pcr.ifthenelse( + S1, + self.ResMaxFl * S_avail / (self.ResEVOL - self.ResPVOL), + self.ResDemFl * S_avail / (self.ResEVOL - self.ResPVOL), + ), + self.StorRES - self.ResEVOL, + ) return Q -#-Simple reservoir + +# -Simple reservoir def QSimple(self, pcr): - Q = pcr.max(self.ResKr * self.StorRES * (self.StorRES / self.ResSmax)**1.5, self.StorRES - self.ResSmax) + Q = pcr.max( + self.ResKr * self.StorRES * (self.StorRES / self.ResSmax) ** 1.5, + self.StorRES - self.ResSmax, + ) return Q - -#-Calculates reservoir outflow and the fraction to release, depending on the type of reservoir (simple or advanced) + + +# -Calculates reservoir outflow and the fraction to release, depending on the type of reservoir (simple or advanced) def QRes(self, pcr): if self.ResSimple and self.ResAdvanced: - Qout = pcr.ifthenelse(self.ResFunc==1, QSimple(self, pcr), pcr.ifthenelse(self.ResFunc==2,\ - QAdv(self, pcr), 0)) + Qout = pcr.ifthenelse( + self.ResFunc == 1, + QSimple(self, pcr), + pcr.ifthenelse(self.ResFunc == 2, QAdv(self, pcr), 0), + ) elif self.ResSimple: - Qout = pcr.ifthenelse(self.ResFunc==1, QSimple(self, pcr), 0) + Qout = pcr.ifthenelse(self.ResFunc == 1, QSimple(self, pcr), 0) else: - Qout = pcr.ifthenelse(self.ResFunc==2, QAdv(self, pcr), 0) - - return Qout \ No newline at end of file + Qout = pcr.ifthenelse(self.ResFunc == 2, QAdv(self, pcr), 0) + + return Qout Index: wflow/wf_DynamicFramework.py =================================================================== diff -u -r4989830126f5201eac4f0d6bb5d4a748fe0c13a3 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision 4989830126f5201eac4f0d6bb5d4a748fe0c13a3) +++ wflow/wf_DynamicFramework.py (.../wf_DynamicFramework.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -636,15 +636,20 @@ setattr( self._userModel(), var, - getattr(self._userModel(), var) * float(self.modelparameters_changes_timestep[cmdd].split('*')[1]) #self.modelparameters_changes_timestep[cmdd], + getattr(self._userModel(), var) + * float( + self.modelparameters_changes_timestep[cmdd].split("*")[1] + ), # self.modelparameters_changes_timestep[cmdd], ) self.logger.warning( "Variable change (apply_timestep) applied to " - + str(var) + " with factor" + self.modelparameters_changes_timestep[cmdd].split('*')[1] + + str(var) + + " with factor" + + self.modelparameters_changes_timestep[cmdd].split("*")[1] ) if self._userModel()._inInitial(): -# import pdb; pdb.set_trace() + # import pdb; pdb.set_trace() for cmdd in self.modelparameters_changes_once: var = cmdd.replace("self._userModel().", "").strip() if not hasattr(self._userModel(), var): @@ -654,14 +659,18 @@ ) else: setattr( - self._userModel(), var, getattr(self._userModel(), var) * float(self.modelparameters_changes_once[cmdd].split('*')[1]) + self._userModel(), + var, + getattr(self._userModel(), var) + * float(self.modelparameters_changes_once[cmdd].split("*")[1]), ) self.logger.warning( "Variable change (apply_once) applied to " - + str(var) + " with factor" + self.modelparameters_changes_once[cmdd].split('*')[1] + + str(var) + + " with factor" + + self.modelparameters_changes_once[cmdd].split("*")[1] ) - def wf_updateparameters(self): """ Update the model Parameters (can be used in static and dynamic part of the model) Index: wflow/wflow_adapt.py =================================================================== diff -u -rc427e9b948a7fbad40a1c828b68355cb277dc5e0 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wflow_adapt.py (.../wflow_adapt.py) (revision c427e9b948a7fbad40a1c828b68355cb277dc5e0) +++ wflow/wflow_adapt.py (.../wflow_adapt.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -150,7 +150,7 @@ 'http://fews.wldelft.nl/schemas/version1.0/pi-schemas/pi_diag.xsd" version="1.2">\n' ) for aline in lines: - translator = aline.maketrans('', '',"><&\"'") + translator = aline.maketrans("", "", "><&\"'") alineesc = aline.translate(translator) fo.write(alineesc) fo.write("\n") Index: wflow/wflow_delwaq.py =================================================================== diff -u -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wflow_delwaq.py (.../wflow_delwaq.py) (revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4) +++ wflow/wflow_delwaq.py (.../wflow_delwaq.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -341,7 +341,12 @@ lowerck = np.absolute(np_ptid) == np.absolute(np_flowto) # mak epointer matrix and add to zero zolumns orgpointer = np.hstack( - (np_ptid, np_flowto, np.zeros((len(np_flowto), 1)), np.zeros((len(np_flowto), 1))) + ( + np_ptid, + np_flowto, + np.zeros((len(np_flowto), 1)), + np.zeros((len(np_flowto), 1)), + ) ) pointer = orgpointer.copy() # Pointer labels: @@ -377,7 +382,9 @@ if bouns == 1: extraboun = np.hstack((bounid, cells, zzerocol, zzerocol)) else: - extraboun = np.vstack((extraboun, np.hstack((bounid, cells, zzerocol, zzerocol)))) + extraboun = np.vstack( + (extraboun, np.hstack((bounid, cells, zzerocol, zzerocol))) + ) pointer_labels = np.hstack((pointer_labels, zzerocol[:, 0] + bouns)) bouns = bouns + 1 start = start + len(cells) @@ -1474,7 +1481,9 @@ # volume for each timestep and number of segments - logger.info("Writing volumes.dat. Nr of points: " + str(np.size(volume_block))) + logger.info( + "Writing volumes.dat. Nr of points: " + str(np.size(volume_block)) + ) dw_WriteSegmentOrExchangeData( i, dwdir + "/includes_flow/volume.dat", volume_block, 1, WriteAscii ) Index: wflow/wflow_routing.py =================================================================== diff -u -r96678c522f159fa455fd7b9eec5983035e274f31 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wflow_routing.py (.../wflow_routing.py) (revision 96678c522f159fa455fd7b9eec5983035e274f31) +++ wflow/wflow_routing.py (.../wflow_routing.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -407,7 +407,7 @@ self.xl, self.yl, self.reallength = pcrut.detRealCellLength( self.ZeroMap, sizeinmetres ) - + # Check if we have reservoirs if hasattr(self, "ReserVoirSimpleLocs"): tt = pcr.pcr2numpy(self.ReserVoirSimpleLocs, 0.0) @@ -937,7 +937,9 @@ self.InflowKinWaveCell = pcr.upstream( self.TopoLdd, self.OldSurfaceRunoff ) - deltasup = float(pcr.mapmaximum(pcr.abs(oldsup - self.SurfaceWaterSupply))) + deltasup = float( + pcr.mapmaximum(pcr.abs(oldsup - self.SurfaceWaterSupply)) + ) if deltasup < self.breakoff or self.nrit >= self.maxitsupply: break Index: wflow/wflow_sbm.py =================================================================== diff -u -ra6b939e08b980c412cc05d2dc0bf707e7f49d2c4 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wflow_sbm.py (.../wflow_sbm.py) (revision a6b939e08b980c412cc05d2dc0bf707e7f49d2c4) +++ wflow/wflow_sbm.py (.../wflow_sbm.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -2430,7 +2430,7 @@ # Calculate simple recharge to groundwater self.Recharge = self.Transfer - self.CapFlux - self.ActEvapSat - + # Now add capflux to the layers one by one (from bottom to top) for n in np.arange(len(self.UStoreLayerThickness) - 1, -1, -1): @@ -2913,7 +2913,9 @@ self.InflowKinWaveCell = pcr.upstream( self.TopoLdd, self.OldSurfaceRunoff ) - deltasup = float(pcr.mapmaximum(pcr.abs(oldsup - self.SurfaceWaterSupply))) + deltasup = float( + pcr.mapmaximum(pcr.abs(oldsup - self.SurfaceWaterSupply)) + ) if deltasup < self.breakoff or self.nrit >= self.maxitsupply: break Index: wflow/wflow_sediment.py =================================================================== diff -u -r2e7a4df821fb7087d463b3cd4219e7552d191921 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wflow_sediment.py (.../wflow_sediment.py) (revision 2e7a4df821fb7087d463b3cd4219e7552d191921) +++ wflow/wflow_sediment.py (.../wflow_sediment.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -37,39 +37,41 @@ from wflow.wf_DynamicFramework import * from wflow.wflow_adapt import * -#import scipy +# import scipy def usage(*args): sys.stdout = sys.stderr - for msg in args: print (msg) - print (__doc__) + for msg in args: + print(msg) + print(__doc__) sys.exit(0) -class WflowModel(DynamicModel): - """ + +class WflowModel(DynamicModel): + """ The user defined model class. This is your work! """ - - def __init__(self, cloneMap,Dir,RunDir,configfile): - """ + + def __init__(self, cloneMap, Dir, RunDir, configfile): + """ *Required* The init function **must** contain what is shown below. Other functionality may be added by you if needed. """ - DynamicModel.__init__(self) - setclone(Dir + "/staticmaps/" + cloneMap) - self.runId=RunDir - self.caseName=Dir - self.Dir = Dir - self.configfile = configfile - self.SaveDir = os.path.join(self.Dir, self.runId) - - def smoothriv(self, river, subcatch, streamorder, param, wlgt): - """ + DynamicModel.__init__(self) + setclone(Dir + "/staticmaps/" + cloneMap) + self.runId = RunDir + self.caseName = Dir + self.Dir = Dir + self.configfile = configfile + self.SaveDir = os.path.join(self.Dir, self.runId) + + def smoothriv(self, river, subcatch, streamorder, param, wlgt): + """ Average a parameter (eg. slope, width...) over river cells only and by streamorder :param river: River map @@ -80,51 +82,84 @@ :return: """ - #Cover river map with zeros - #Rivercov = ifthenelse(river == 1, 1.0, scalar(0.0)) - #Rivercov = ifthen(subcatch > 0, cover(boolean(river), 0)) - Rivercov = cover(boolean(river),0) - Rivercov = scalar(Rivercov) - - order3 = ifthenelse(streamorder == 3, param, scalar(0.0)) - riv3 = ifthenelse(streamorder == 3, Rivercov, scalar(0.0)) - order4 = ifthenelse(streamorder == 4, param, scalar(0.0)) - riv4 = ifthenelse(streamorder == 4, Rivercov, scalar(0.0)) - order5 = ifthenelse(streamorder == 5, param, scalar(0.0)) - riv5 = ifthenelse(streamorder == 5, Rivercov, scalar(0.0)) - order6 = ifthenelse(streamorder == 6, param, scalar(0.0)) - riv6 = ifthenelse(streamorder == 6, Rivercov, scalar(0.0)) - order7 = ifthenelse(streamorder == 7, param, scalar(0.0)) - riv7 = ifthenelse(streamorder == 7, Rivercov, scalar(0.0)) - order8 = ifthenelse(streamorder == 8, param, scalar(0.0)) - riv8 = ifthenelse(streamorder == 8, Rivercov, scalar(0.0)) - order9 = ifthenelse(streamorder == 9, param, scalar(0.0)) - riv9 = ifthenelse(streamorder == 9, Rivercov, scalar(0.0)) - - unitmap = riv9 * 0.0 + 1.0 - - #self.W = ifthenelse(self.streamorder == 1, windowaverage(order1, celllength()*2) * windowtotal(riv1, celllength()*2), self.W) - param = ifthenelse(streamorder == 9, windowaverage(order9, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv9, celllength()*wlgt), param) - param = ifthenelse(streamorder == 8, windowaverage(order8, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv8, celllength()*wlgt), param) - param = ifthenelse(streamorder == 7, windowaverage(order7, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv7, celllength()*wlgt), param) - param = ifthenelse(streamorder == 6, windowaverage(order6, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv6, celllength()*wlgt), param) - param = ifthenelse(streamorder == 5, windowaverage(order5, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv5, celllength()*wlgt), param) - param = ifthenelse(streamorder == 4, windowaverage(order4, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv4, celllength()*wlgt), param) - param = ifthenelse(streamorder == 3, windowaverage(order3, celllength()*wlgt) * windowtotal(unitmap, celllength()*wlgt) \ - / windowtotal(riv3, celllength()*wlgt), param) + # Cover river map with zeros + # Rivercov = ifthenelse(river == 1, 1.0, scalar(0.0)) + # Rivercov = ifthen(subcatch > 0, cover(boolean(river), 0)) + Rivercov = cover(boolean(river), 0) + Rivercov = scalar(Rivercov) + order3 = ifthenelse(streamorder == 3, param, scalar(0.0)) + riv3 = ifthenelse(streamorder == 3, Rivercov, scalar(0.0)) + order4 = ifthenelse(streamorder == 4, param, scalar(0.0)) + riv4 = ifthenelse(streamorder == 4, Rivercov, scalar(0.0)) + order5 = ifthenelse(streamorder == 5, param, scalar(0.0)) + riv5 = ifthenelse(streamorder == 5, Rivercov, scalar(0.0)) + order6 = ifthenelse(streamorder == 6, param, scalar(0.0)) + riv6 = ifthenelse(streamorder == 6, Rivercov, scalar(0.0)) + order7 = ifthenelse(streamorder == 7, param, scalar(0.0)) + riv7 = ifthenelse(streamorder == 7, Rivercov, scalar(0.0)) + order8 = ifthenelse(streamorder == 8, param, scalar(0.0)) + riv8 = ifthenelse(streamorder == 8, Rivercov, scalar(0.0)) + order9 = ifthenelse(streamorder == 9, param, scalar(0.0)) + riv9 = ifthenelse(streamorder == 9, Rivercov, scalar(0.0)) - return param - + unitmap = riv9 * 0.0 + 1.0 - def parameters(self): - """ + # self.W = ifthenelse(self.streamorder == 1, windowaverage(order1, celllength()*2) * windowtotal(riv1, celllength()*2), self.W) + param = ifthenelse( + streamorder == 9, + windowaverage(order9, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv9, celllength() * wlgt), + param, + ) + param = ifthenelse( + streamorder == 8, + windowaverage(order8, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv8, celllength() * wlgt), + param, + ) + param = ifthenelse( + streamorder == 7, + windowaverage(order7, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv7, celllength() * wlgt), + param, + ) + param = ifthenelse( + streamorder == 6, + windowaverage(order6, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv6, celllength() * wlgt), + param, + ) + param = ifthenelse( + streamorder == 5, + windowaverage(order5, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv5, celllength() * wlgt), + param, + ) + param = ifthenelse( + streamorder == 4, + windowaverage(order4, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv4, celllength() * wlgt), + param, + ) + param = ifthenelse( + streamorder == 3, + windowaverage(order3, celllength() * wlgt) + * windowtotal(unitmap, celllength() * wlgt) + / windowtotal(riv3, celllength() * wlgt), + param, + ) + + return param + + 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 @@ -140,25 +175,68 @@ """ - modelparameters = [] - - # Input time series from ini file - self.P_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Precipitation", "/inmaps/P") #precipitation - self.Int_mapstack = self.Dir + configget(self.config, "inputmapstacks", "Interception", "/inmaps/int") #rainfall interception - self.SR_mapstack = self.Dir + configget(self.config, "inputmapstacks", "SurfaceRunoff", "/inmaps/run") #surface runoff - self.WL_mapstack = self.Dir + configget(self.config, "inputmapstacks", "WaterLevel", "/inmaps/levKin") #water level + modelparameters = [] - # Meteo and other forcing from wflow_sbm - modelparameters.append(self.ParamType(name="Precipitation",stack=self.P_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="Interception",stack=self.Int_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[])) - modelparameters.append(self.ParamType(name="SurfaceRunoff",stack=self.SR_mapstack,type="timeseries",default=0.0,verbose=True,lookupmaps=[])) - modelparameters.append(self.ParamType(name="WaterLevel",stack=self.WL_mapstack,type="timeseries",default=0.0,verbose=False,lookupmaps=[])) - + # Input time series from ini file + self.P_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Precipitation", "/inmaps/P" + ) # precipitation + self.Int_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "Interception", "/inmaps/int" + ) # rainfall interception + self.SR_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "SurfaceRunoff", "/inmaps/run" + ) # surface runoff + self.WL_mapstack = self.Dir + configget( + self.config, "inputmapstacks", "WaterLevel", "/inmaps/levKin" + ) # water level - return modelparameters + # Meteo and other forcing from wflow_sbm + modelparameters.append( + self.ParamType( + name="Precipitation", + stack=self.P_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="Interception", + stack=self.Int_mapstack, + type="timeseries", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="SurfaceRunoff", + stack=self.SR_mapstack, + type="timeseries", + default=0.0, + verbose=True, + lookupmaps=[], + ) + ) + modelparameters.append( + self.ParamType( + name="WaterLevel", + stack=self.WL_mapstack, + type="timeseries", + default=0.0, + verbose=False, + lookupmaps=[], + ) + ) - def stateVariables(self): - """ + return modelparameters + + def stateVariables(self): + """ *Required* Returns a list of state variables that are essential to the model. @@ -172,23 +250,38 @@ :var OutSedLoad: Total sediment load transported out of river cells [ton/cell/timestep] :var RivStoreSed: Deposited sediment storage in the river [ton/cell/timestep] """ - - if self.RunRiverModel == 0: - states = [] - else: - states = ['SedLoad', 'ClayLoad', 'SiltLoad', - 'SandLoad', 'SaggLoad', 'LaggLoad', - 'GravLoad', 'OutSedLoad', 'OutClayLoad', - 'OutSiltLoad', 'OutSandLoad', 'OutSaggLoad', - 'OutLaggLoad', 'OutGravLoad', 'RivStoreSed', - 'RivStoreClay', 'RivStoreSilt', 'RivStoreSand', - 'RivStoreSagg', 'RivStoreLagg', 'RivStoreGrav'] - - return states - - - def supplyCurrentTime(self): - """ + + if self.RunRiverModel == 0: + states = [] + else: + states = [ + "SedLoad", + "ClayLoad", + "SiltLoad", + "SandLoad", + "SaggLoad", + "LaggLoad", + "GravLoad", + "OutSedLoad", + "OutClayLoad", + "OutSiltLoad", + "OutSandLoad", + "OutSaggLoad", + "OutLaggLoad", + "OutGravLoad", + "RivStoreSed", + "RivStoreClay", + "RivStoreSilt", + "RivStoreSand", + "RivStoreSagg", + "RivStoreLagg", + "RivStoreGrav", + ] + + return states + + def supplyCurrentTime(self): + """ *Optional* Supplies the current time in seconds after the start of the run @@ -200,11 +293,13 @@ - time in seconds since the start of the model run """ - - return self.currentTimeStep() * int(configget(self.config,'model','timestepsecs','86400')) - - def suspend(self): - """ + + return self.currentTimeStep() * int( + configget(self.config, "model", "timestepsecs", "86400") + ) + + def suspend(self): + """ *Required* Suspends the model to disk. All variables needed to restart the model @@ -213,17 +308,16 @@ This function is required. """ - - self.logger.info("Saving initial conditions...") - #: It is advised to use the wf_suspend() function - #: here which will suspend the variables that are given by stateVariables - #: function. - self.wf_suspend(self.SaveDir + "/outstate/") - - def initial(self): - - """ + self.logger.info("Saving initial conditions...") + #: It is advised to use the wf_suspend() function + #: here which will suspend the variables that are given by stateVariables + #: function. + self.wf_suspend(self.SaveDir + "/outstate/") + + def initial(self): + + """ *Required* Initial part of the model, executed only once. It reads all static model @@ -253,412 +347,714 @@ :var CovRiver: Factor representing bank erodibility reduction due to its vegetation [-] """ - #: pcraster option to calculate with units or cells. Not really an issue - #: in this model but always good to keep in mind. - setglobaloption("unittrue") + #: pcraster option to calculate with units or cells. Not really an issue + #: in this model but always good to keep in mind. + setglobaloption("unittrue") - self.timestepsecs = int(configget(self.config,'model','timestepsecs','86400')) - sizeinmetres = int(configget(self.config, "layout", "sizeinmetres", "0")) - self.basetimestep=86400 - - # Reads all parameter from disk - self.wf_updateparameters() - - ''' Read static model parameters/maps from ini file ''' - - self.intbl = configget(self.config, "model", "intbl", "intbl") - self.RunRiverModel = int(configget(self.config, "model", "runrivermodel", "1")) - self.slopecorr = int(configget(self.config, "model", "slopecorr", "0")) - self.UsleKMethod = int(configget(self.config, "model", "uslekmethod", "2")) - self.UsleCMethod = int(configget(self.config, "model", "uslecmethod", "2")) - self.RainErodMethod = int(configget(self.config, "model", "rainerodmethod", "1")) - self.LandTransportMethod = int(configget(self.config, "model", "landtransportmethod", "1")) - self.RivTransportMethod = int(configget(self.config, "model", "rivtransportmethod", "1")) - if self.RunRiverModel == 1: - self.LandTransportMethod = 1 - - - #Static maps to use - wflow_dem = configget(self.config, "model", "wflow_dem", "staticmaps/wflow_dem.map") - self.Altitude = self.wf_readmap(os.path.join(self.Dir,wflow_dem), 0.0, fail=True) - wflow_landuse = configget(self.config, "model", "wflow_landuse", "staticmaps/wflow_landuse.map") - self.LandUse = self.wf_readmap(os.path.join(self.Dir,wflow_landuse), 0.0, fail=True) - wflow_soil = configget(self.config, "model", "wflow_soil", "staticmaps/wflow_soil.map") - self.Soil = self.wf_readmap(os.path.join(self.Dir,wflow_soil), 0.0, fail=False) - wflow_subcatch = configget(self.config, "model", "wflow_subcatch", "staticmaps/wflow_subcatch.map") - self.TopoId = self.wf_readmap(os.path.join(self.Dir,wflow_subcatch), 0.0, fail=True) - wflow_Hype = configget(self.config, "model", "wflow_Hype", "staticmaps/SUBID-HYPE-Rhine.map") - self.HypeId = self.wf_readmap(os.path.join(self.Dir,wflow_Hype), 0.0, fail=False) - wflow_ldd = configget(self.config, "model", "wflow_ldd", "staticmaps/wflow_ldd.map") - self.TopoLdd = self.wf_readmap(os.path.join(self.Dir,wflow_ldd), 0.0, fail=True) - wflow_river = configget(self.config, "model", "wflow_river", "staticmaps/wflow_river.map") - self.River = self.wf_readmap(os.path.join(self.Dir,wflow_river), 0.0, fail=True) - wflow_riverwidth = configget(self.config, "model", "wflow_riverwidth", "staticmaps/RiverWidth.map") - self.RiverWidth = self.wf_readmap(os.path.join(self.Dir,wflow_riverwidth), 0.0, fail=True) - self.RiverWidth = ifthenelse(self.River == 1, self.RiverWidth, scalar(0.0)) - wflow_dcl = configget(self.config, "model", "wflow_dcl", "staticmaps/DCL.map") - self.DCL = self.wf_readmap(os.path.join(self.Dir,wflow_dcl), 0.0, fail=True) # Drain/River length - wflow_streamorder = configget(self.config, "model", "wflow_streamorder", "staticmaps/wflow_streamorder.map") - self.streamorder = self.wf_readmap(os.path.join(self.Dir,wflow_streamorder), 0.0, fail=True) - - #Soil - wflow_clay = configget(self.config, "model", "wflow_clay", "staticmaps/percent_clay.map") - self.PercentClay = self.wf_readmap(os.path.join(self.Dir,wflow_clay), 0.1, fail=True) - wflow_silt = configget(self.config, "model", "wflow_silt", "staticmaps/percent_silt.map") - self.PercentSilt = self.wf_readmap(os.path.join(self.Dir,wflow_silt), 0.1, fail=True) - wflow_oc = configget(self.config, "model", "wflow_oc", "staticmaps/percent_oc.map") - self.PercentOC = self.wf_readmap(os.path.join(self.Dir,wflow_oc), 0.1, fail=False) -# wflow_bulk = configget(self.config, "model", "wflow_bulk", "staticmaps/bulk_density.map") -# self.BulkDensity = self.wf_readmap(os.path.join(self.Dir,wflow_bulk), 0.1, fail=False) - - - ''' Read tbl parameters and link them to landuse, soil, subcatch...) ''' - - # Soil impervious area - self.PathFrac = self.readtblDefault(self.Dir + "/" + self.intbl + "/PathFrac.tbl", self.LandUse, self.TopoId, self.Soil, 0.01) - - #Soil model parameters - self.ErosOv = self.readtblDefault(self.Dir + "/" + self.intbl + "/eros_ov.tbl", self.LandUse, self.TopoId, self.Soil, 0.90) - - if self.RunRiverModel == 1: - #River model parameters - self.D50River = self.readtblFlexDefault(self.Dir + "/" + self.intbl + "/D50_River.tbl", 0.050, wflow_streamorder) - self.D50River = ifthenelse(self.River == 1, self.D50River, scalar(0.0)) - self.CovRiver = self.readtblDefault(self.Dir + "/" + self.intbl + "/cov_River.tbl", self.LandUse, self.TopoId, self.Soil, 1.0) - - ''' Determine global variables ''' - # Map with zeros - self.ZeroMap = 0.0 * self.Altitude - - # Subcatchment map - self.subcatch = ordinal(self.TopoId) - self.subcatch = ifthen(self.subcatch > 0, self.subcatch) - - # HYPE subcatchment map - self.HYPEcatch = ordinal(self.HypeId) - self.HYPEcatch = ifthen(self.HYPEcatch > 0, self.HYPEcatch) - - # Determine real slope, cell length and cell area - self.xl, self.yl, self.reallength = pcrut.detRealCellLength(self.ZeroMap, sizeinmetres) - self.Slope = slope(self.Altitude) - self.Slope = max(0.00001, self.Slope * celllength() / self.reallength) - - self.cellareaKm = (self.reallength/1000.)**2 - self.UpArea = accuflux(self.TopoLdd, self.cellareaKm) - - #Sine of the slope - self.sinSlope = sin(atan(self.Slope)) - - #If necessary reduce dem and slope for river cells - if self.slopecorr == 1: - self.AltitudeMin = self.wf_readmap(os.path.join(self.Dir,"staticmaps/wflow_demmin.map"), 0.0, fail=True) - #In 90m SRTM correct sea cells (-32768) with wflow_dem - self.AltitudeMin = ifthenelse(self.AltitudeMin == -32768, scalar(0.0), self.AltitudeMin) - #Take this minimum value only for river cells (rest is original DEM) - self.RiverDem = ifthenelse(self.River == 1, self.AltitudeMin, scalar(0.0)) - self.RiverDem = ifthen(self.subcatch > 0, cover(self.RiverDem, scalar(0.0))) - #Compute corrected slope for river cells - self.Rivercov = ifthenelse(self.River == 1, 1.0, scalar(0.0)) - self.Rivercov = ifthen(self.subcatch > 0, cover(self.Rivercov, scalar(0.0))) - self.nrupcell = upstream(self.TopoLdd, self.Rivercov) - self.RiverSlope = (upstream(self.TopoLdd, self.RiverDem)/self.nrupcell - downstream(self.TopoLdd, self.RiverDem)) \ - / (2 * self.reallength) - self.RiverSlope = max(0.0003, self.RiverSlope) - - #Cover the non river cells with original slope/dem - self.RiverDem = cover(self.RiverDem, self.Altitude) - self.RiverSlope = cover(self.RiverSlope, self.Slope) - - else : - self.RiverDem = self.Altitude - self.RiverSlope = self.Slope - - #Correct slope with drain length - drainlength = detdrainlength(self.TopoLdd, self.xl, self.yl) - riverslopecor = drainlength / self.DCL - self.RiverSlope = self.RiverSlope * riverslopecor - - #Smooth river slope - self.RiverSlope = self.smoothriv(self.River, subcatch, self.streamorder, self.RiverSlope, 6) - -# #Calculate RiverWidth -# upstr = catchmenttotal(1, self.TopoLdd) -# Qmax = 2200 -# Qscale = upstr / mapmaximum(upstr) * Qmax -# alf = 120 -# self.N = self.readtblDefault(self.Dir + "/" + self.intbl + "/N.tbl",self.LandUse,self.TopoId,self.Soil,0.072) -# self.NRiver = self.readtblFlexDefault(self.Dir + "/" + self.intbl + "/N_River.tbl", 0.036, wflow_streamorder) -# self.N = ifthen(subcatch > 0, cover(ifthenelse(self.River, self.NRiver, self.N), self.N)) -# self.W = (alf * (alf + 2.0) ** (0.6666666667)) ** (0.375)* Qscale ** (0.375) \ -# * (max(0.0001, windowaverage(self.Slope, celllength() * 4.0))) ** (-0.1875)* self.N ** (0.375) -# self.WRiv = ifthen(self.River,((alf * (alf + 2.0) ** (0.6666666667)) ** (0.375)* Qscale ** (0.375)* \ -# self.RiverSlope ** (-0.1875)* self.N ** (0.375))) -# #Smooth river width -# self.WRiv = self.smoothriv(self.River, subcatch, self.streamorder, self.WRiv, 4) -# self.W = ifthen(subcatch > 0, cover(ifthenelse(self.River, self.WRiv, self.W), self.W)) -# self.RiverWidth = self.W - - - ''' Determine variables for the soil loss model ''' - - # Canopy gap fraction based on LAI or input table - if hasattr(self,"LAI"): - if not hasattr(self,"Kext"): - logging.error("Kext (canopy extinction coefficient) not defined! Needed becausee LAI is defined.") - logging.error("Please add it to the modelparameters section. e.g.:") - logging.error("Kext=inmaps/clim/LCtoExtinctionCoefficient.tbl,tbl,0.5,1,inmaps/clim/LC.map") - self.CanopyGapFraction = exp(-self.Kext * self.LAI) - else: - self.CanopyGapFraction = self.readtblDefault(self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl", self.LandUse, self.TopoId, self.Soil, 0.1) - - - # Determine sand content - self.PercentSand = 100 - (self.PercentClay + self.PercentSilt) - - # Calculate detachability of the soil (k) [g/J] - self.ErosK = ifthenelse(pcrand(self.PercentClay>=40.,pcrand(self.PercentSand>=20.,self.PercentSand<=45.)),2.0, - ifthenelse(pcrand(self.PercentClay>=27.,pcrand(self.PercentSand>=20.,self.PercentSand<=45.)),1.7, - ifthenelse(pcrand(self.PercentSilt<=40.,self.PercentSand<=20.),2.0, - ifthenelse(pcrand(self.PercentSilt>40.,self.PercentClay>=40.),1.6, - ifthenelse(pcrand(self.PercentClay>=35.,self.PercentSand>=45.),1.9, - ifthenelse(pcrand(self.PercentClay>=27.,self.PercentSand<20.),1.6, - ifthenelse(pcrand(self.PercentClay<=10.,self.PercentSilt>=80.),1.2, - ifthenelse(self.PercentSilt>=50,1.5, - ifthenelse(pcrand(self.PercentClay>=7.,pcrand(self.PercentSand<=52.,self.PercentSilt>=28.)),2.0, - ifthenelse(self.PercentClay>=20.,2.1, - ifthenelse(self.PercentClay>=self.PercentSand-70.,2.6, - ifthenelse(self.PercentClay>=(2.*self.PercentSand)-170.,3,scalar(1.9))))))))))))) - - # Compute USLE K factor - if self.UsleKMethod == 1: - self.UsleC =self.wf_readmap(os.path.join(self.Dir,"staticmaps/USLE_K.map"), 0.1, fail=True) - if self.UsleKMethod == 2: - # Calculate USLE K factor (from Renard et al. 1997, with the geometric mean particle diameter Dg) - self.Dg = exp(0.01*(self.PercentClay*ln(0.001) + self.PercentSilt*ln(0.025) + self.PercentSand*ln(0.999))) #[mm] - self.UsleK = 0.0034 + 0.0405 * exp(-1/2 * ((log10(self.Dg)+1.659)/0.7101)**2) - #Remove possible outliers - self.UsleK = max(0.0,self.UsleK) - - if self.UsleKMethod == 3: - # Calculate USLE K factor (from Williams and Renard 1983, EPIC: a new method for assessing erosion's effect on soil productivity) - self.SN = (1-self.PercentSand)/100 - self.UsleK = (0.2 + 0.3*exp(-0.0256*self.PercentSand * (1 - self.PercentSilt/100))) \ - * (self.PercentSilt/(max(0.01,self.PercentClay+self.PercentSilt)))**0.3 \ - * (1 - (0.25*self.PercentOC)/(self.PercentOC + exp(3.72 - 2.95*self.PercentOC))) \ - * (1 - (0.7*self.SN) / (self.SN + exp(-5.51 + 22.9*self.SN))) - #Remove possible outliers - self.UsleK = max(0.0,self.UsleK) - - #Compute USLE C factor - if self.UsleCMethod == 1: - self.UsleC =self.wf_readmap(os.path.join(self.Dir,"staticmaps/USLE_C.map"), 0.1, fail=True) - if self.UsleCMethod == 2: - # USLE C factor map based on land use (from Gericke 2015, Soil loss estimation and empirical relationships for sediment delivery ratios of European river catchments) - self.UsleC = self.readtblDefault(self.Dir + "/" + self.intbl + "/USLE_C.tbl", self.LandUse, self.TopoId, self.Soil, 1.0) + self.timestepsecs = int( + configget(self.config, "model", "timestepsecs", "86400") + ) + sizeinmetres = int(configget(self.config, "layout", "sizeinmetres", "0")) + self.basetimestep = 86400 + # Reads all parameter from disk + self.wf_updateparameters() - #Parameters for either EUROSEM or ANSWERS rainfall erosion - if self.RainErodMethod == 1: #EUROSEM - #Canopy height [m] - self.CanopyHeight =self.wf_readmap(os.path.join(self.Dir,"staticmaps/canopy_height.map"), 1.0, fail=True) - #Coefficient - self.ErosSpl = self.readtblDefault(self.Dir + "/" + self.intbl + "/eros_spl.tbl", self.LandUse, self.TopoId, self.Soil, 2.0) - if self.RainErodMethod == 2: #ANSWERS - #Coefficient - self.ErosSpl = self.readtblDefault(self.Dir + "/" + self.intbl + "/eros_spl.tbl", self.LandUse, self.TopoId, self.Soil, 0.108) - - - ''' Variables for inland transport model ''' - # If the river model is run, use Yalin's equation with particle differentiation. - # If just the soil loss model is run, use Govers equation without particle differentiation - - if (self.LandTransportMethod == 2 or self.LandTransportMethod == 3): - #Calculation of D50 and fraction of fine and very fine sand (fvfs) from Fooladmand et al, 2006 - self.PercentSand999 = self.PercentSand * (999-25)/(1000-25) #%sand with a mean radius of 999um instead of 1000 - self.vd50 = ln((1/((self.PercentClay+self.PercentSilt)/100) -1) / (1/(self.PercentClay/100) -1)) - self.wd50 = ln((1/((self.PercentClay+self.PercentSilt+self.PercentSand999)/100) -1) / (1/(self.PercentClay/100) -1)) - ad50 = 1 / (-3.727699) # 1 / ln((25-1)/(999-1)) - bd50 = ad50 * 3.17805 # ad50 / ln((25-1)/1) - self.cd50 = ad50 * ln(self.vd50 / self.wd50) - self.ud50 = (-self.vd50)**(1 - bd50) / (-self.wd50)**(-bd50) - - self.D50 = 1 + (-1/self.ud50 * ln(1 / (1/(self.PercentClay/100) -1)))**(1/self.cd50) #[um] - self.D50 = self.D50 / 1000 #[mm] - - # #Fraction of fine/very fine sand (d=125um) and coarse sand - # self.percent_fvfs = 100 * 1 / (1+(1/(self.percent_clay/100)-1)*exp(-self.ud50*(125-1)**self.cd50)) - self.percent_silt - self.percent_clay #[%] - # self.percent_coars = self.percent_sand - self.percent_fvfs - - #Calculate Govers transport capacity coefficients - self.cGovers = ((self.D50*1000 + 5) / 0.32)**(-0.6) - self.nGovers = ((self.D50*1000 + 5) / 300)**0.25 - - elif self.LandTransportMethod == 1: - #Determine sediment size distribution, estimated from primary particle size distribution (Foster et al., 1980) - self.FracClay = 0.20 * self.PercentClay/100 - self.FracSilt = 0.13 * self.PercentSilt/100 - self.FracSand = self.PercentSand/100 * (1 - self.PercentClay/100)**2.4 - self.FracSagg = ifthenelse(self.PercentClay<25, 2.0 * self.PercentClay/100, - ifthenelse(self.PercentClay>50, 0.57, - 0.28*(self.PercentClay/100-0.25)+0.5)) - self.FracLagg = 1.0 - self.FracClay - self.FracSilt - self.FracSand - self.FracSagg - - - ''' Variables for the river transport model ''' - - if self.RunRiverModel == 1: - #River particle size distribution (estimated with SWAT method) !!! D50 can be calibrated - self.FracClayRiv = ifthenelse(self.D50River <= 0.005, 0.65, ifthenelse(self.D50River > 2.0, 0.005, scalar(0.15))) - self.FracSiltRiv = ifthenelse(pcrand(self.D50River > 0.005, self.D50River <= 0.050), 0.65, scalar(0.15)) - self.FracSandRiv = ifthenelse(pcrand(self.D50River > 0.050, self.D50River <= 2.0), 0.65, scalar(0.15)) - self.FracGravRiv = ifthenelse(self.D50River > 2.0, 0.65, scalar(0.05)) - - #Parameters of Bagnold transport formula - if self.RivTransportMethod == 2: - self.cBagnold = self.readtblDefault(self.Dir + "/" + self.intbl + "/c_Bagnold.tbl", self.LandUse, self.TopoId, self.Soil, 0.0015) - self.expBagnold = self.readtblDefault(self.Dir + "/" + self.intbl + "/exp_Bagnold.tbl", self.LandUse, self.TopoId, self.Soil, 1.4) - - #Parameters of Kodatie transport formula - if self.RivTransportMethod == 3: - self.aK = ifthenelse(self.D50River <= 0.05, 281.4, ifthenelse(self.D50River <= 0.25, 2829.6, - ifthenelse(self.D50River <= 2, 2123.4, scalar(431884.8)))) - self.bK = ifthenelse(self.D50River <= 0.05, 2.622, ifthenelse(self.D50River <= 0.25, 3.646, - ifthenelse(self.D50River <= 2, 3.300, scalar(1.0)))) - self.cK = ifthenelse(self.D50River <= 0.05, 0.182, ifthenelse(self.D50River <= 0.25, 0.406, - ifthenelse(self.D50River <= 2, 0.468, scalar(1.0)))) - self.dK = ifthenelse(self.D50River <= 0.05, 0.0, ifthenelse(self.D50River <= 0.25, 0.412, - ifthenelse(self.D50River <= 2, 0.613, scalar(2.0)))) - - - #Critical bed and bank shear stresses and erodibilities [N/m2] [m3/N.s] - #Bank from Julian & Torres 2006 + Hanson & Simon 2001 - self.TCrBank = (0.1 + 0.1779*(100*self.FracClayRiv+100*self.FracSiltRiv)+0.0028*(100*self.FracClayRiv+100*self.FracSiltRiv)**2 \ - -2.34*10**(-5)*(100*self.FracClayRiv+100*self.FracSiltRiv)**3)*self.CovRiver - self.KdBank = 0.2 * self.TCrBank ** (-0.5) * 10**(-6) - #Bed from Shields diagram - E = ((2.65-1)*9.81*(self.D50River*10**(-3))**3/(1*10**(-12)))**(0.33) - self.TCrBed = (2.65-1)*9.81*self.D50River*(0.13*E**(-0.392)*exp(-0.015*E**2) + 0.045*(1-exp(-0.068*E))) - self.KdBed = 0.2 * self.TCrBed**(-0.5) * 10**(-6) - - - ''' Variables for reservoirs model ''' - - if hasattr(self, "ReserVoirSimpleLocs") or hasattr(self, "ReserVoirComplexLocs"): - self.ReserVoirLocs = self.ZeroMap - self.filter_Eros_TC = self.ZeroMap + 1.0 - - if hasattr(self, "ReserVoirSimpleLocs"): - # Check if we have simple and or complex reservoirs - tt_simple = pcr2numpy(self.ReserVoirSimpleLocs, 0.0) - self.nrresSimple = tt_simple.max() - self.ReserVoirLocs = self.ReserVoirLocs + cover(scalar(self.ReserVoirSimpleLocs)) - areamap = self.reallength * self.reallength - res_area = areatotal(spatial(areamap), self.ReservoirSimpleAreas) + """ Read static model parameters/maps from ini file """ - resarea_pnt = ifthen(boolean(self.ReserVoirSimpleLocs), res_area) - self.ResSimpleArea = ifthenelse( - cover(self.ResSimpleArea, scalar(0.0)) > 0, - self.ResSimpleArea, - cover(resarea_pnt, scalar(0.0)), + self.intbl = configget(self.config, "model", "intbl", "intbl") + self.RunRiverModel = int(configget(self.config, "model", "runrivermodel", "1")) + self.slopecorr = int(configget(self.config, "model", "slopecorr", "0")) + self.UsleKMethod = int(configget(self.config, "model", "uslekmethod", "2")) + self.UsleCMethod = int(configget(self.config, "model", "uslecmethod", "2")) + self.RainErodMethod = int( + configget(self.config, "model", "rainerodmethod", "1") ) - self.filter_Eros_TC = ifthenelse( - boolean(cover(res_area, scalar(0.0))), res_area * 0.0, self.filter_Eros_TC + self.LandTransportMethod = int( + configget(self.config, "model", "landtransportmethod", "1") ) - else: - self.nrresSimple = 0 + self.RivTransportMethod = int( + configget(self.config, "model", "rivtransportmethod", "1") + ) + if self.RunRiverModel == 1: + self.LandTransportMethod = 1 - if hasattr(self, "ReserVoirComplexLocs"): - tt_complex = pcr2numpy(self.ReserVoirComplexLocs, 0.0) - self.nrresComplex = tt_complex.max() - self.ReserVoirLocs = self.ReserVoirLocs + cover(scalar(self.ReserVoirComplexLocs)) - res_area = cover(scalar(self.ReservoirComplexAreas), 0.0) - self.filter_Eros_TC = ifthenelse( - res_area > 0, res_area * 0.0, self.filter_Eros_TC + # Static maps to use + wflow_dem = configget( + self.config, "model", "wflow_dem", "staticmaps/wflow_dem.map" ) - else: - self.nrresComplex = 0 - - if (self.nrresSimple + self.nrresComplex) > 0: - self.ReserVoirLocs = ordinal(self.ReserVoirLocs) - self.logger.info( - "A total of " - + str(self.nrresSimple) - + " simple reservoirs and " - + str(self.nrresComplex) - + " complex reservoirs found." + self.Altitude = self.wf_readmap( + os.path.join(self.Dir, wflow_dem), 0.0, fail=True ) - self.ReserVoirDownstreamLocs = downstream(self.TopoLdd, self.ReserVoirLocs) - self.TopoLddOrg = lddrepair(cover(ifthen(boolean(self.ReserVoirLocs), ldd(5)), self.TopoLdd)) + wflow_landuse = configget( + self.config, "model", "wflow_landuse", "staticmaps/wflow_landuse.map" + ) + self.LandUse = self.wf_readmap( + os.path.join(self.Dir, wflow_landuse), 0.0, fail=True + ) + wflow_soil = configget( + self.config, "model", "wflow_soil", "staticmaps/wflow_soil.map" + ) + self.Soil = self.wf_readmap(os.path.join(self.Dir, wflow_soil), 0.0, fail=False) + wflow_subcatch = configget( + self.config, "model", "wflow_subcatch", "staticmaps/wflow_subcatch.map" + ) + self.TopoId = self.wf_readmap( + os.path.join(self.Dir, wflow_subcatch), 0.0, fail=True + ) + wflow_Hype = configget( + self.config, "model", "wflow_Hype", "staticmaps/SUBID-HYPE-Rhine.map" + ) + self.HypeId = self.wf_readmap( + os.path.join(self.Dir, wflow_Hype), 0.0, fail=False + ) + wflow_ldd = configget( + self.config, "model", "wflow_ldd", "staticmaps/wflow_ldd.map" + ) + self.TopoLdd = self.wf_readmap( + os.path.join(self.Dir, wflow_ldd), 0.0, fail=True + ) + wflow_river = configget( + self.config, "model", "wflow_river", "staticmaps/wflow_river.map" + ) + self.River = self.wf_readmap( + os.path.join(self.Dir, wflow_river), 0.0, fail=True + ) + wflow_riverwidth = configget( + self.config, "model", "wflow_riverwidth", "staticmaps/RiverWidth.map" + ) + self.RiverWidth = self.wf_readmap( + os.path.join(self.Dir, wflow_riverwidth), 0.0, fail=True + ) + self.RiverWidth = ifthenelse(self.River == 1, self.RiverWidth, scalar(0.0)) + wflow_dcl = configget(self.config, "model", "wflow_dcl", "staticmaps/DCL.map") + self.DCL = self.wf_readmap( + os.path.join(self.Dir, wflow_dcl), 0.0, fail=True + ) # Drain/River length + wflow_streamorder = configget( + self.config, + "model", + "wflow_streamorder", + "staticmaps/wflow_streamorder.map", + ) + self.streamorder = self.wf_readmap( + os.path.join(self.Dir, wflow_streamorder), 0.0, fail=True + ) - tt_filter = pcr2numpy(self.filter_Eros_TC, 1.0) - self.filterResArea = tt_filter.min() - - - self.logger.info("Starting Dynamic run...") + # Soil + wflow_clay = configget( + self.config, "model", "wflow_clay", "staticmaps/percent_clay.map" + ) + self.PercentClay = self.wf_readmap( + os.path.join(self.Dir, wflow_clay), 0.1, fail=True + ) + wflow_silt = configget( + self.config, "model", "wflow_silt", "staticmaps/percent_silt.map" + ) + self.PercentSilt = self.wf_readmap( + os.path.join(self.Dir, wflow_silt), 0.1, fail=True + ) + wflow_oc = configget( + self.config, "model", "wflow_oc", "staticmaps/percent_oc.map" + ) + self.PercentOC = self.wf_readmap( + os.path.join(self.Dir, wflow_oc), 0.1, fail=False + ) + # wflow_bulk = configget(self.config, "model", "wflow_bulk", "staticmaps/bulk_density.map") + # self.BulkDensity = self.wf_readmap(os.path.join(self.Dir,wflow_bulk), 0.1, fail=False) + """ Read tbl parameters and link them to landuse, soil, subcatch...) """ - def resume(self): - """ + # Soil impervious area + self.PathFrac = self.readtblDefault( + self.Dir + "/" + self.intbl + "/PathFrac.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 0.01, + ) + + # Soil model parameters + self.ErosOv = self.readtblDefault( + self.Dir + "/" + self.intbl + "/eros_ov.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 0.90, + ) + + if self.RunRiverModel == 1: + # River model parameters + self.D50River = self.readtblFlexDefault( + self.Dir + "/" + self.intbl + "/D50_River.tbl", 0.050, wflow_streamorder + ) + self.D50River = ifthenelse(self.River == 1, self.D50River, scalar(0.0)) + self.CovRiver = self.readtblDefault( + self.Dir + "/" + self.intbl + "/cov_River.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 1.0, + ) + + """ Determine global variables """ + # Map with zeros + self.ZeroMap = 0.0 * self.Altitude + + # Subcatchment map + self.subcatch = ordinal(self.TopoId) + self.subcatch = ifthen(self.subcatch > 0, self.subcatch) + + # HYPE subcatchment map + self.HYPEcatch = ordinal(self.HypeId) + self.HYPEcatch = ifthen(self.HYPEcatch > 0, self.HYPEcatch) + + # Determine real slope, cell length and cell area + self.xl, self.yl, self.reallength = pcrut.detRealCellLength( + self.ZeroMap, sizeinmetres + ) + self.Slope = slope(self.Altitude) + self.Slope = max(0.00001, self.Slope * celllength() / self.reallength) + + self.cellareaKm = (self.reallength / 1000.0) ** 2 + self.UpArea = accuflux(self.TopoLdd, self.cellareaKm) + + # Sine of the slope + self.sinSlope = sin(atan(self.Slope)) + + # If necessary reduce dem and slope for river cells + if self.slopecorr == 1: + self.AltitudeMin = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/wflow_demmin.map"), 0.0, fail=True + ) + # In 90m SRTM correct sea cells (-32768) with wflow_dem + self.AltitudeMin = ifthenelse( + self.AltitudeMin == -32768, scalar(0.0), self.AltitudeMin + ) + # Take this minimum value only for river cells (rest is original DEM) + self.RiverDem = ifthenelse(self.River == 1, self.AltitudeMin, scalar(0.0)) + self.RiverDem = ifthen(self.subcatch > 0, cover(self.RiverDem, scalar(0.0))) + # Compute corrected slope for river cells + self.Rivercov = ifthenelse(self.River == 1, 1.0, scalar(0.0)) + self.Rivercov = ifthen(self.subcatch > 0, cover(self.Rivercov, scalar(0.0))) + self.nrupcell = upstream(self.TopoLdd, self.Rivercov) + self.RiverSlope = ( + upstream(self.TopoLdd, self.RiverDem) / self.nrupcell + - downstream(self.TopoLdd, self.RiverDem) + ) / (2 * self.reallength) + self.RiverSlope = max(0.0003, self.RiverSlope) + + # Cover the non river cells with original slope/dem + self.RiverDem = cover(self.RiverDem, self.Altitude) + self.RiverSlope = cover(self.RiverSlope, self.Slope) + + else: + self.RiverDem = self.Altitude + self.RiverSlope = self.Slope + + # Correct slope with drain length + drainlength = detdrainlength(self.TopoLdd, self.xl, self.yl) + riverslopecor = drainlength / self.DCL + self.RiverSlope = self.RiverSlope * riverslopecor + + # Smooth river slope + self.RiverSlope = self.smoothriv( + self.River, subcatch, self.streamorder, self.RiverSlope, 6 + ) + + # #Calculate RiverWidth + # upstr = catchmenttotal(1, self.TopoLdd) + # Qmax = 2200 + # Qscale = upstr / mapmaximum(upstr) * Qmax + # alf = 120 + # self.N = self.readtblDefault(self.Dir + "/" + self.intbl + "/N.tbl",self.LandUse,self.TopoId,self.Soil,0.072) + # self.NRiver = self.readtblFlexDefault(self.Dir + "/" + self.intbl + "/N_River.tbl", 0.036, wflow_streamorder) + # self.N = ifthen(subcatch > 0, cover(ifthenelse(self.River, self.NRiver, self.N), self.N)) + # self.W = (alf * (alf + 2.0) ** (0.6666666667)) ** (0.375)* Qscale ** (0.375) \ + # * (max(0.0001, windowaverage(self.Slope, celllength() * 4.0))) ** (-0.1875)* self.N ** (0.375) + # self.WRiv = ifthen(self.River,((alf * (alf + 2.0) ** (0.6666666667)) ** (0.375)* Qscale ** (0.375)* \ + # self.RiverSlope ** (-0.1875)* self.N ** (0.375))) + # #Smooth river width + # self.WRiv = self.smoothriv(self.River, subcatch, self.streamorder, self.WRiv, 4) + # self.W = ifthen(subcatch > 0, cover(ifthenelse(self.River, self.WRiv, self.W), self.W)) + # self.RiverWidth = self.W + + """ Determine variables for the soil loss model """ + + # Canopy gap fraction based on LAI or input table + if hasattr(self, "LAI"): + if not hasattr(self, "Kext"): + logging.error( + "Kext (canopy extinction coefficient) not defined! Needed becausee LAI is defined." + ) + logging.error("Please add it to the modelparameters section. e.g.:") + logging.error( + "Kext=inmaps/clim/LCtoExtinctionCoefficient.tbl,tbl,0.5,1,inmaps/clim/LC.map" + ) + self.CanopyGapFraction = exp(-self.Kext * self.LAI) + else: + self.CanopyGapFraction = self.readtblDefault( + self.Dir + "/" + self.intbl + "/CanopyGapFraction.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 0.1, + ) + + # Determine sand content + self.PercentSand = 100 - (self.PercentClay + self.PercentSilt) + + # Calculate detachability of the soil (k) [g/J] + self.ErosK = ifthenelse( + pcrand( + self.PercentClay >= 40.0, + pcrand(self.PercentSand >= 20.0, self.PercentSand <= 45.0), + ), + 2.0, + ifthenelse( + pcrand( + self.PercentClay >= 27.0, + pcrand(self.PercentSand >= 20.0, self.PercentSand <= 45.0), + ), + 1.7, + ifthenelse( + pcrand(self.PercentSilt <= 40.0, self.PercentSand <= 20.0), + 2.0, + ifthenelse( + pcrand(self.PercentSilt > 40.0, self.PercentClay >= 40.0), + 1.6, + ifthenelse( + pcrand(self.PercentClay >= 35.0, self.PercentSand >= 45.0), + 1.9, + ifthenelse( + pcrand( + self.PercentClay >= 27.0, self.PercentSand < 20.0 + ), + 1.6, + ifthenelse( + pcrand( + self.PercentClay <= 10.0, + self.PercentSilt >= 80.0, + ), + 1.2, + ifthenelse( + self.PercentSilt >= 50, + 1.5, + ifthenelse( + pcrand( + self.PercentClay >= 7.0, + pcrand( + self.PercentSand <= 52.0, + self.PercentSilt >= 28.0, + ), + ), + 2.0, + ifthenelse( + self.PercentClay >= 20.0, + 2.1, + ifthenelse( + self.PercentClay + >= self.PercentSand - 70.0, + 2.6, + ifthenelse( + self.PercentClay + >= (2.0 * self.PercentSand) + - 170.0, + 3, + scalar(1.9), + ), + ), + ), + ), + ), + ), + ), + ), + ), + ), + ), + ) + + # Compute USLE K factor + if self.UsleKMethod == 1: + self.UsleC = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/USLE_K.map"), 0.1, fail=True + ) + if self.UsleKMethod == 2: + # Calculate USLE K factor (from Renard et al. 1997, with the geometric mean particle diameter Dg) + self.Dg = exp( + 0.01 + * ( + self.PercentClay * ln(0.001) + + self.PercentSilt * ln(0.025) + + self.PercentSand * ln(0.999) + ) + ) # [mm] + self.UsleK = 0.0034 + 0.0405 * exp( + -1 / 2 * ((log10(self.Dg) + 1.659) / 0.7101) ** 2 + ) + # Remove possible outliers + self.UsleK = max(0.0, self.UsleK) + + if self.UsleKMethod == 3: + # Calculate USLE K factor (from Williams and Renard 1983, EPIC: a new method for assessing erosion's effect on soil productivity) + self.SN = (1 - self.PercentSand) / 100 + self.UsleK = ( + ( + 0.2 + + 0.3 + * exp(-0.0256 * self.PercentSand * (1 - self.PercentSilt / 100)) + ) + * (self.PercentSilt / (max(0.01, self.PercentClay + self.PercentSilt))) + ** 0.3 + * ( + 1 + - (0.25 * self.PercentOC) + / (self.PercentOC + exp(3.72 - 2.95 * self.PercentOC)) + ) + * (1 - (0.7 * self.SN) / (self.SN + exp(-5.51 + 22.9 * self.SN))) + ) + # Remove possible outliers + self.UsleK = max(0.0, self.UsleK) + + # Compute USLE C factor + if self.UsleCMethod == 1: + self.UsleC = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/USLE_C.map"), 0.1, fail=True + ) + if self.UsleCMethod == 2: + # USLE C factor map based on land use (from Gericke 2015, Soil loss estimation and empirical relationships for sediment delivery ratios of European river catchments) + self.UsleC = self.readtblDefault( + self.Dir + "/" + self.intbl + "/USLE_C.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 1.0, + ) + + # Parameters for either EUROSEM or ANSWERS rainfall erosion + if self.RainErodMethod == 1: # EUROSEM + # Canopy height [m] + self.CanopyHeight = self.wf_readmap( + os.path.join(self.Dir, "staticmaps/canopy_height.map"), 1.0, fail=True + ) + # Coefficient + self.ErosSpl = self.readtblDefault( + self.Dir + "/" + self.intbl + "/eros_spl.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 2.0, + ) + if self.RainErodMethod == 2: # ANSWERS + # Coefficient + self.ErosSpl = self.readtblDefault( + self.Dir + "/" + self.intbl + "/eros_spl.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 0.108, + ) + + """ Variables for inland transport model """ + # If the river model is run, use Yalin's equation with particle differentiation. + # If just the soil loss model is run, use Govers equation without particle differentiation + + if self.LandTransportMethod == 2 or self.LandTransportMethod == 3: + # Calculation of D50 and fraction of fine and very fine sand (fvfs) from Fooladmand et al, 2006 + self.PercentSand999 = ( + self.PercentSand * (999 - 25) / (1000 - 25) + ) #%sand with a mean radius of 999um instead of 1000 + self.vd50 = ln( + (1 / ((self.PercentClay + self.PercentSilt) / 100) - 1) + / (1 / (self.PercentClay / 100) - 1) + ) + self.wd50 = ln( + ( + 1 + / ( + (self.PercentClay + self.PercentSilt + self.PercentSand999) + / 100 + ) + - 1 + ) + / (1 / (self.PercentClay / 100) - 1) + ) + ad50 = 1 / (-3.727699) # 1 / ln((25-1)/(999-1)) + bd50 = ad50 * 3.17805 # ad50 / ln((25-1)/1) + self.cd50 = ad50 * ln(self.vd50 / self.wd50) + self.ud50 = (-self.vd50) ** (1 - bd50) / (-self.wd50) ** (-bd50) + + self.D50 = 1 + ( + -1 / self.ud50 * ln(1 / (1 / (self.PercentClay / 100) - 1)) + ) ** ( + 1 / self.cd50 + ) # [um] + self.D50 = self.D50 / 1000 # [mm] + + # #Fraction of fine/very fine sand (d=125um) and coarse sand + # self.percent_fvfs = 100 * 1 / (1+(1/(self.percent_clay/100)-1)*exp(-self.ud50*(125-1)**self.cd50)) - self.percent_silt - self.percent_clay #[%] + # self.percent_coars = self.percent_sand - self.percent_fvfs + + # Calculate Govers transport capacity coefficients + self.cGovers = ((self.D50 * 1000 + 5) / 0.32) ** (-0.6) + self.nGovers = ((self.D50 * 1000 + 5) / 300) ** 0.25 + + elif self.LandTransportMethod == 1: + # Determine sediment size distribution, estimated from primary particle size distribution (Foster et al., 1980) + self.FracClay = 0.20 * self.PercentClay / 100 + self.FracSilt = 0.13 * self.PercentSilt / 100 + self.FracSand = self.PercentSand / 100 * (1 - self.PercentClay / 100) ** 2.4 + self.FracSagg = ifthenelse( + self.PercentClay < 25, + 2.0 * self.PercentClay / 100, + ifthenelse( + self.PercentClay > 50, + 0.57, + 0.28 * (self.PercentClay / 100 - 0.25) + 0.5, + ), + ) + self.FracLagg = ( + 1.0 - self.FracClay - self.FracSilt - self.FracSand - self.FracSagg + ) + + """ Variables for the river transport model """ + + if self.RunRiverModel == 1: + # River particle size distribution (estimated with SWAT method) !!! D50 can be calibrated + self.FracClayRiv = ifthenelse( + self.D50River <= 0.005, + 0.65, + ifthenelse(self.D50River > 2.0, 0.005, scalar(0.15)), + ) + self.FracSiltRiv = ifthenelse( + pcrand(self.D50River > 0.005, self.D50River <= 0.050), + 0.65, + scalar(0.15), + ) + self.FracSandRiv = ifthenelse( + pcrand(self.D50River > 0.050, self.D50River <= 2.0), 0.65, scalar(0.15) + ) + self.FracGravRiv = ifthenelse(self.D50River > 2.0, 0.65, scalar(0.05)) + + # Parameters of Bagnold transport formula + if self.RivTransportMethod == 2: + self.cBagnold = self.readtblDefault( + self.Dir + "/" + self.intbl + "/c_Bagnold.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 0.0015, + ) + self.expBagnold = self.readtblDefault( + self.Dir + "/" + self.intbl + "/exp_Bagnold.tbl", + self.LandUse, + self.TopoId, + self.Soil, + 1.4, + ) + + # Parameters of Kodatie transport formula + if self.RivTransportMethod == 3: + self.aK = ifthenelse( + self.D50River <= 0.05, + 281.4, + ifthenelse( + self.D50River <= 0.25, + 2829.6, + ifthenelse(self.D50River <= 2, 2123.4, scalar(431884.8)), + ), + ) + self.bK = ifthenelse( + self.D50River <= 0.05, + 2.622, + ifthenelse( + self.D50River <= 0.25, + 3.646, + ifthenelse(self.D50River <= 2, 3.300, scalar(1.0)), + ), + ) + self.cK = ifthenelse( + self.D50River <= 0.05, + 0.182, + ifthenelse( + self.D50River <= 0.25, + 0.406, + ifthenelse(self.D50River <= 2, 0.468, scalar(1.0)), + ), + ) + self.dK = ifthenelse( + self.D50River <= 0.05, + 0.0, + ifthenelse( + self.D50River <= 0.25, + 0.412, + ifthenelse(self.D50River <= 2, 0.613, scalar(2.0)), + ), + ) + + # Critical bed and bank shear stresses and erodibilities [N/m2] [m3/N.s] + # Bank from Julian & Torres 2006 + Hanson & Simon 2001 + self.TCrBank = ( + 0.1 + + 0.1779 * (100 * self.FracClayRiv + 100 * self.FracSiltRiv) + + 0.0028 * (100 * self.FracClayRiv + 100 * self.FracSiltRiv) ** 2 + - 2.34 + * 10 ** (-5) + * (100 * self.FracClayRiv + 100 * self.FracSiltRiv) ** 3 + ) * self.CovRiver + self.KdBank = 0.2 * self.TCrBank ** (-0.5) * 10 ** (-6) + # Bed from Shields diagram + E = ( + (2.65 - 1) + * 9.81 + * (self.D50River * 10 ** (-3)) ** 3 + / (1 * 10 ** (-12)) + ) ** (0.33) + self.TCrBed = ( + (2.65 - 1) + * 9.81 + * self.D50River + * ( + 0.13 * E ** (-0.392) * exp(-0.015 * E ** 2) + + 0.045 * (1 - exp(-0.068 * E)) + ) + ) + self.KdBed = 0.2 * self.TCrBed ** (-0.5) * 10 ** (-6) + + """ Variables for reservoirs model """ + + if hasattr(self, "ReserVoirSimpleLocs") or hasattr( + self, "ReserVoirComplexLocs" + ): + self.ReserVoirLocs = self.ZeroMap + self.filter_Eros_TC = self.ZeroMap + 1.0 + + if hasattr(self, "ReserVoirSimpleLocs"): + # Check if we have simple and or complex reservoirs + tt_simple = pcr2numpy(self.ReserVoirSimpleLocs, 0.0) + self.nrresSimple = tt_simple.max() + self.ReserVoirLocs = self.ReserVoirLocs + cover( + scalar(self.ReserVoirSimpleLocs) + ) + areamap = self.reallength * self.reallength + res_area = areatotal(spatial(areamap), self.ReservoirSimpleAreas) + + resarea_pnt = ifthen(boolean(self.ReserVoirSimpleLocs), res_area) + self.ResSimpleArea = ifthenelse( + cover(self.ResSimpleArea, scalar(0.0)) > 0, + self.ResSimpleArea, + cover(resarea_pnt, scalar(0.0)), + ) + self.filter_Eros_TC = ifthenelse( + boolean(cover(res_area, scalar(0.0))), + res_area * 0.0, + self.filter_Eros_TC, + ) + else: + self.nrresSimple = 0 + + if hasattr(self, "ReserVoirComplexLocs"): + tt_complex = pcr2numpy(self.ReserVoirComplexLocs, 0.0) + self.nrresComplex = tt_complex.max() + self.ReserVoirLocs = self.ReserVoirLocs + cover( + scalar(self.ReserVoirComplexLocs) + ) + res_area = cover(scalar(self.ReservoirComplexAreas), 0.0) + self.filter_Eros_TC = ifthenelse( + res_area > 0, res_area * 0.0, self.filter_Eros_TC + ) + else: + self.nrresComplex = 0 + + if (self.nrresSimple + self.nrresComplex) > 0: + self.ReserVoirLocs = ordinal(self.ReserVoirLocs) + self.logger.info( + "A total of " + + str(self.nrresSimple) + + " simple reservoirs and " + + str(self.nrresComplex) + + " complex reservoirs found." + ) + self.ReserVoirDownstreamLocs = downstream(self.TopoLdd, self.ReserVoirLocs) + self.TopoLddOrg = lddrepair( + cover(ifthen(boolean(self.ReserVoirLocs), ldd(5)), self.TopoLdd) + ) + + tt_filter = pcr2numpy(self.filter_Eros_TC, 1.0) + self.filterResArea = tt_filter.min() + + self.logger.info("Starting Dynamic run...") + + def resume(self): + """ *Required* This function is required. Read initial state maps (they are output of a previous call to suspend()). The implementation shown here is the most basic setup needed. """ - self.logger.info("Reading initial conditions...") - #: It is advised to use the wf_resume() function - #: here which pick up the variable save by a call to wf_suspend() + self.logger.info("Reading initial conditions...") + #: It is advised to use the wf_resume() function + #: here which pick up the variable save by a call to wf_suspend() - if self.reinit == 1: - self.logger.warn("Setting initial states to default") - if self.RunRiverModel == 1: - self.SedLoad = self.ZeroMap - self.ClayLoad = self.ZeroMap - self.SiltLoad = self.ZeroMap - self.SandLoad = self.ZeroMap - self.SaggLoad = self.ZeroMap - self.LaggLoad = self.ZeroMap - self.GravLoad = self.ZeroMap - - self.OutSedLoad = self.ZeroMap - self.OutClayLoad = self.ZeroMap - self.OutSiltLoad = self.ZeroMap - self.OutSandLoad = self.ZeroMap - self.OutSaggLoad = self.ZeroMap - self.OutLaggLoad = self.ZeroMap - self.OutGravLoad = self.ZeroMap - - self.RivStoreSed = self.ZeroMap - self.RivStoreClay = self.ZeroMap - self.RivStoreSilt = self.ZeroMap - self.RivStoreSand = self.ZeroMap - self.RivStoreSagg = self.ZeroMap - self.RivStoreLagg = self.ZeroMap - self.RivStoreGrav = self.ZeroMap - else: - self.logger.info("Setting initial conditions from state files") - self.wf_resume(os.path.join(self.Dir, "instate")) + if self.reinit == 1: + self.logger.warn("Setting initial states to default") + if self.RunRiverModel == 1: + self.SedLoad = self.ZeroMap + self.ClayLoad = self.ZeroMap + self.SiltLoad = self.ZeroMap + self.SandLoad = self.ZeroMap + self.SaggLoad = self.ZeroMap + self.LaggLoad = self.ZeroMap + self.GravLoad = self.ZeroMap + self.OutSedLoad = self.ZeroMap + self.OutClayLoad = self.ZeroMap + self.OutSiltLoad = self.ZeroMap + self.OutSandLoad = self.ZeroMap + self.OutSaggLoad = self.ZeroMap + self.OutLaggLoad = self.ZeroMap + self.OutGravLoad = self.ZeroMap - def default_summarymaps(self): - """ + self.RivStoreSed = self.ZeroMap + self.RivStoreClay = self.ZeroMap + self.RivStoreSilt = self.ZeroMap + self.RivStoreSand = self.ZeroMap + self.RivStoreSagg = self.ZeroMap + self.RivStoreLagg = self.ZeroMap + self.RivStoreGrav = self.ZeroMap + else: + self.logger.info("Setting initial conditions from state files") + self.wf_resume(os.path.join(self.Dir, "instate")) + + def default_summarymaps(self): + """ *Optional* Return a default list of variables to report as summary maps in the outsum dir. The ini file has more options, including average and sum """ - lst = ['self.PercentSand', 'self.UsleK', - 'self.UsleC', 'self.ErodK', - 'self.D50River','self.FracClayRiv', - 'self.filter_Eros_TC', 'self.RiverSlope', - 'self.RiverWidth', 'self.RiverDem', - 'self.UpArea', 'self.cellareaKm'] - - return lst + lst = [ + "self.PercentSand", + "self.UsleK", + "self.UsleC", + "self.ErodK", + "self.D50River", + "self.FracClayRiv", + "self.filter_Eros_TC", + "self.RiverSlope", + "self.RiverWidth", + "self.RiverDem", + "self.UpArea", + "self.cellareaKm", + ] - def dynamic(self): - """ + return lst + + def dynamic(self): + """ *Required* This is where all the time dependent functions are executed. Time dependent @@ -690,550 +1086,1274 @@ """ - self.wf_updateparameters() - self.Precipitation = max(0.0,self.Precipitation) + self.wf_updateparameters() + self.Precipitation = max(0.0, self.Precipitation) - ########################################################################## - # Soil loss model ############################# - ########################################################################## - - - ''' Splash / Rainfall erosion ''' - #From EUROSEM - if self.RainErodMethod == 1: - # calculate rainfall intensity [mm/h] - self.rintnsty = self.Precipitation/(self.timestepsecs/3600) - - #Kinectic energy of direct throughfall [J/m^2/mm] - #self.KeDirect = max(11.87 + 8.73 * log10(max(0.0001,self.rintnsty)), 0.0) #basis used in USLE - self.KeDirect = max(8.95 + 8.44 * log10(max(0.0001,self.rintnsty)), 0.0) #variant, most used in distributed models - - #Kinetic energy of leaf drainage [J/m^2/mm] - pheff = 0.5 * self.CanopyHeight # [m] - self.KeLeaf = max((15.8 * pheff ** 0.5) - 5.87, 0.0) - - #Depths of rainfall (total, leaf drainage, direct) [mm] - #rdepth_tot = max(self.Precipitation/self.timestepsecs, 0.0) - rDepthTot = max(self.Precipitation, 0.0) - rDepthLeaf = max(rDepthTot * 0.1 * self.CanopyGapFraction, 0.0) #stemflow - rDepthDirect = max(rDepthTot - rDepthLeaf - self.Interception, 0.0) #throughfall - - #Total kinetic energy by rainfall [J/m^2] - self.KeTotal = (rDepthDirect * self.KeDirect + rDepthLeaf * self.KeLeaf) *0.001 - - # Rainfall/Splash erosion - self.SedSpl = self.ErosK * self.KeTotal * exp(-self.ErosSpl * self.WaterLevel) # [g/m^2] - self.Sedspl = self.cellareaKm * self.SedSpl #* self.timestepsecs # [ton/cell/timestep] - - if self.RainErodMethod == 2: - # calculate rainfall intensity [mm/min] - self.rintnsty = self.Precipitation/(self.timestepsecs/60) - - # Splash erosion - self.SedSpl = 0.108 * self.UsleC * self.UsleK * self.cellareaKm * 10**6 * self.rintnsty**2 # [kg/min] - self.SedSpl = self.timestepsecs/60.0 * 10**(-3) * self.SedSpl # [ton/timestep] - - # Remove the impervious areas - self.SedSpl = self.SedSpl * (1.- self.PathFrac) - # Remove nodata values - self.SedSpl = ifthen(self.subcatch > 0, cover(self.SedSpl, 0.0)) - - - ''' Overland flow erosion from ANSWERS''' - # Only calculate overland flow erosion outside of river cells - self.OvRun = cover(ifthenelse(self.River == 1, 0.0, self.SurfaceRunoff), self.SurfaceRunoff) #[m3/s] - self.OvRunRate = self.OvRun * 60 / self.reallength #[m2/min] - - # Overland flow erosion - #For a wide range of slope, it is better to use the sine of slope rather than tangeant - self.SedOv = self.ErosOv * self.UsleC * self.UsleK * self.cellareaKm * 10**6 * self.sinSlope * self.OvRunRate # [kg/min] - self.SedOv = self.timestepsecs/60.0 * 10**(-3) * self.SedOv # [ton/timestep] - - # Remove the impervious areas - self.SedOv = self.SedOv * (1.- self.PathFrac) - # Remove nodata values - self.SedOv = ifthen(self.subcatch > 0, cover(self.SedOv, 0.0)) - - - ''' Total soil detachment ''' - self.SoilLoss = self.SedSpl + self.SedOv #[ton/cell/timestep] - - #Remove land erosion for reservoir cells - if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: + ########################################################################## + # Soil loss model ############################# + ########################################################################## + + """ Splash / Rainfall erosion """ + # From EUROSEM + if self.RainErodMethod == 1: + # calculate rainfall intensity [mm/h] + self.rintnsty = self.Precipitation / (self.timestepsecs / 3600) + + # Kinectic energy of direct throughfall [J/m^2/mm] + # self.KeDirect = max(11.87 + 8.73 * log10(max(0.0001,self.rintnsty)), 0.0) #basis used in USLE + self.KeDirect = max( + 8.95 + 8.44 * log10(max(0.0001, self.rintnsty)), 0.0 + ) # variant, most used in distributed models + + # Kinetic energy of leaf drainage [J/m^2/mm] + pheff = 0.5 * self.CanopyHeight # [m] + self.KeLeaf = max((15.8 * pheff ** 0.5) - 5.87, 0.0) + + # Depths of rainfall (total, leaf drainage, direct) [mm] + # rdepth_tot = max(self.Precipitation/self.timestepsecs, 0.0) + rDepthTot = max(self.Precipitation, 0.0) + rDepthLeaf = max(rDepthTot * 0.1 * self.CanopyGapFraction, 0.0) # stemflow + rDepthDirect = max( + rDepthTot - rDepthLeaf - self.Interception, 0.0 + ) # throughfall + + # Total kinetic energy by rainfall [J/m^2] + self.KeTotal = ( + rDepthDirect * self.KeDirect + rDepthLeaf * self.KeLeaf + ) * 0.001 + + # Rainfall/Splash erosion + self.SedSpl = ( + self.ErosK * self.KeTotal * exp(-self.ErosSpl * self.WaterLevel) + ) # [g/m^2] + self.Sedspl = ( + self.cellareaKm * self.SedSpl + ) # * self.timestepsecs # [ton/cell/timestep] + + if self.RainErodMethod == 2: + # calculate rainfall intensity [mm/min] + self.rintnsty = self.Precipitation / (self.timestepsecs / 60) + + # Splash erosion + self.SedSpl = ( + 0.108 + * self.UsleC + * self.UsleK + * self.cellareaKm + * 10 ** 6 + * self.rintnsty ** 2 + ) # [kg/min] + self.SedSpl = ( + self.timestepsecs / 60.0 * 10 ** (-3) * self.SedSpl + ) # [ton/timestep] + + # Remove the impervious areas + self.SedSpl = self.SedSpl * (1.0 - self.PathFrac) + # Remove nodata values + self.SedSpl = ifthen(self.subcatch > 0, cover(self.SedSpl, 0.0)) + + """ Overland flow erosion from ANSWERS""" + # Only calculate overland flow erosion outside of river cells + self.OvRun = cover( + ifthenelse(self.River == 1, 0.0, self.SurfaceRunoff), self.SurfaceRunoff + ) # [m3/s] + self.OvRunRate = self.OvRun * 60 / self.reallength # [m2/min] + + # Overland flow erosion + # For a wide range of slope, it is better to use the sine of slope rather than tangeant + self.SedOv = ( + self.ErosOv + * self.UsleC + * self.UsleK + * self.cellareaKm + * 10 ** 6 + * self.sinSlope + * self.OvRunRate + ) # [kg/min] + self.SedOv = ( + self.timestepsecs / 60.0 * 10 ** (-3) * self.SedOv + ) # [ton/timestep] + + # Remove the impervious areas + self.SedOv = self.SedOv * (1.0 - self.PathFrac) + # Remove nodata values + self.SedOv = ifthen(self.subcatch > 0, cover(self.SedOv, 0.0)) + + """ Total soil detachment """ + self.SoilLoss = self.SedSpl + self.SedOv # [ton/cell/timestep] + + # Remove land erosion for reservoir cells + if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: self.SedSpl = self.filter_Eros_TC * self.SedSpl self.SedOv = self.filter_Eros_TC * self.SedOv self.SoilLoss = self.filter_Eros_TC * self.SoilLoss - - - ########################################################################## - # Inland sediment routing model ############################# - ########################################################################## - - # If the river model is run, use Yalin's equation with particle differentiation. - # If just the soil loss model is run, use Govers equation without particle differentiation - - ''' Transport of overland flow with no particle differenciation using Govers equation''' - - if (self.LandTransportMethod == 2 or self.LandTransportMethod == 3): - if self.LandTransportMethod == 2: - #Unit stream power - self.velocity = cover(ifthenelse(self.WaterLevel > 0, self.OvRun / (self.reallength * self.WaterLevel), 0.0), 0.0) #[m/s] - self.omega = 10 * self.sinSlope * 100 * self.velocity #[cm/s] - #self.omega = self.sinSlope * 100 * self.velocity #[cm/s] - - #Transport capacity from Govers, 1990 - self.TCf = ifthenelse(self.omega > 0.4, self.cGovers * (self.omega - 0.4)**self.nGovers * 2650, 0.0) #[kg/m3] - #self.TC = self.TCf / (1 - self.TCf/2650) #[kg/m3] - #self.TC = max(self.TC, 2650) - self.TC = self.TCf * self.OvRun * self.timestepsecs * 10**(-3) #[ton/cell/timestep] - # Remove nodata values - self.TC = ifthen(self.subcatch > 0, cover(self.TC, 0.0)) - #Assume that eroded soil on lake cells all reach the river cells of the reservoir - if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: - self.TC = ifthenelse(pcrand(self.filter_Eros_TC == 0, self.River == 0), 10**9, self.TC) - - elif self.LandTransportMethod == 3: - #Transport capacity from Yalin - self.OvLevel = cover(ifthenelse(self.River == 1, 0.0, self.WaterLevel), self.WaterLevel) #[m] - self.delta = max(self.OvLevel * self.sinSlope / (self.D50*10**(-3)*(2650/1000-1)) /0.06 -1, 0.0) - self.TC = self.reallength/self.OvRun * (2650-1000) * self.D50*10**(-3) * (9.81*self.OvLevel*self.sinSlope) * \ - 0.635 * self.delta * (1 - ln(1 + self.delta*2.45/(2650/1000)**0.4*0.06**0.5)/self.delta*2.45/(2650/1000)**0.4*0.06**0.5) #[kg/m3] - self.TC = ifthen(self.subcatch > 0, cover(self.TC*self.OvRun*self.timestepsecs*10**(-3), 0.0)) #[ton/cell/timestep] - #Assume that eroded soil on lake cells all reach the river cells of the reservoir - if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: - self.TC = ifthenelse(pcrand(self.filter_Eros_TC == 0, self.River == 0), 10**9, self.TC) - - #To get total sediment input from land into the river systems, river cells transport all sediment to the output (huge TC) - self.TCRiv = cover(ifthenelse(self.River == 1, 10**9, self.TC), self.TC) - #Transported sediment over the land - self.SedFlux = accucapacityflux(self.TopoLdd, self.SoilLoss, self.TCRiv) #[ton/cell/tinestep] - #Deposited sediment over the land - self.SedDep = accucapacitystate(self.TopoLdd, self.SoilLoss, self.TCRiv) - - #Sediment amount reaching each river cell ''' - self.SedDep2 = accucapacitystate(self.TopoLdd, self.SoilLoss, self.TC) - #Remove inland deposition - self.OvSed = cover(ifthenelse(self.River == 1, self.SedDep2, 0.0), 0.0) - - #Sum the results per subcatchment (input for D-WAQ) [kg/ha/timestep] - self.HYPEOvSedCatch = areatotal(self.OvSed, self.HYPEcatch) #[ton/cell/timestep] - self.HYPEOvSedCatch = self.HYPEOvSedCatch * 1000 / (self.cellareaKm *100) #[kg/ha/timestep] - - ''' Transport of overland flow with particle differenciation using Yalin equation''' - else: - #Determine the eroded amount of clay/silt/sand/aggregates on each cell [ton/cell/timestep] - self.LandErodClay = self.SoilLoss * self.FracClay - self.LandErodSilt = self.SoilLoss * self.FracSilt - self.LandErodSand = self.SoilLoss * self.FracSand - self.LandErodSagg = self.SoilLoss * self.FracSagg - self.LandErodLagg = self.SoilLoss * self.FracLagg - - #Water level of overland flow - self.OvLevel = cover(ifthenelse(self.River == 1, 0.0, self.WaterLevel), self.WaterLevel) #[m] - - #Delta parameter of Yalin for each particle class - self.DClay = max(self.OvLevel * self.sinSlope / (2*10**(-6)*(2650/1000-1)) /0.06 -1, 0.0) - self.DSilt = max(self.OvLevel * self.sinSlope / (10*10**(-6)*(2650/1000-1)) /0.06 -1, 0.0) - self.DSand = max(self.OvLevel * self.sinSlope / (200*10**(-6)*(2650/1000-1)) /0.06 -1, 0.0) - self.DSagg = max(self.OvLevel * self.sinSlope / (30*10**(-6)*(2650/1000-1)) /0.06 -1, 0.0) - self.DLagg = max(self.OvLevel * self.sinSlope / (500*10**(-6)*(2650/1000-1)) /0.06 -1, 0.0) - - #Total transportability - self.Dtot = self.DClay + self.DSilt + self.DSand + self.DSagg + self.DLagg - - #Yalin Transport capacity of overland flow for each particle class - self.TCClay = self.reallength/self.OvRun * (2650-1000) * 2*10**(-6) * (9.81*self.OvLevel*self.sinSlope) * self.DClay/self.Dtot * \ - 0.635 * self.DClay * (1 - ln(1 + self.DClay*2.45/(2650/1000)**0.4*0.06**0.5)/self.DClay*2.45/(2650/1000)**0.4*0.06**0.5) #[kg/m3] - self.TCClay = ifthen(self.subcatch > 0, cover(self.TCClay*self.OvRun*self.timestepsecs*10**(-3), 0.0)) #[ton/cell/timestep] - - self.TCSilt = self.reallength/self.OvRun * (2650-1000) * 10*10**(-6) * (9.81*self.OvLevel*self.sinSlope) * self.DSilt/self.Dtot * \ - 0.635 * self.DSilt * (1 - ln(1 + self.DSilt*2.45/(2650/1000)**0.4*0.06**0.5)/self.DSilt*2.45/(2650/1000)**0.4*0.06**0.5) #[kg/m3] - self.TCSilt = ifthen(self.subcatch > 0, cover(self.TCSilt*self.OvRun*self.timestepsecs*10**(-3), 0.0)) #[ton/cell/timestep] - - self.TCSand = self.reallength/self.OvRun * (2650-1000) * 200*10**(-6) * (9.81*self.OvLevel*self.sinSlope) * self.DSand/self.Dtot * \ - 0.635 * self.DSand * (1 - ln(1 + self.DSand*2.45/(2650/1000)**0.4*0.06**0.5)/self.DSand*2.45/(2650/1000)**0.4*0.06**0.5) #[kg/m3] - self.TCSand = ifthen(self.subcatch > 0, cover(self.TCSand*self.OvRun*self.timestepsecs*10**(-3), 0.0)) #[ton/cell/timestep] - - self.TCSagg = self.reallength/self.OvRun * (2650-1000) * 30*10**(-6) * (9.81*self.OvLevel*self.sinSlope) * self.DSagg/self.Dtot * \ - 0.635 * self.DSagg * (1 - ln(1 + self.DSagg*2.45/(2650/1000)**0.4*0.06**0.5)/self.DSagg*2.45/(2650/1000)**0.4*0.06**0.5) #[kg/m3] - self.TCSagg = ifthen(self.subcatch > 0, cover(self.TCSagg*self.OvRun*self.timestepsecs*10**(-3), 0.0)) #[ton/cell/timestep] - - self.TCLagg = self.reallength/self.OvRun * (2650-1000) * 500*10**(-6) * (9.81*self.OvLevel*self.sinSlope) * self.DLagg/self.Dtot * \ - 0.635 * self.DLagg * (1 - ln(1 + self.DLagg*2.45/(2650/1000)**0.4*0.06**0.5)/self.DLagg*2.45/(2650/1000)**0.4*0.06**0.5) #[kg/m3] - self.TCLagg = ifthen(self.subcatch > 0, cover(self.TCLagg*self.OvRun*self.timestepsecs*10**(-3), 0.0)) #[ton/cell/timestep] - - #Assume that eroded soil on lake cells all reach the river cells of the reservoir - if (self.nrresSimple + self.nrresComplex) > 0: - self.TCClay = cover(ifthenelse(self.filter_Eros_TC == 0, 10**9, self.TCClay),self.TCClay) - self.TCClay = cover(ifthenelse(self.River == 1, 0.0, self.TCClay),self.TCClay) - self.TCSilt = cover(ifthenelse(self.filter_Eros_TC == 0, 10**9, self.TCSilt),self.TCSilt) - self.TCSilt = cover(ifthenelse(self.River == 1, 0.0, self.TCSilt),self.TCSilt) - self.TCSand = cover(ifthenelse(self.filter_Eros_TC == 0, 10**9, self.TCSand),self.TCSand) - self.TCSand = cover(ifthenelse(self.River == 1, 0.0, self.TCSand),self.TCSand) - self.TCSagg = cover(ifthenelse(self.filter_Eros_TC == 0, 10**9, self.TCSagg),self.TCSagg) - self.TCSagg = cover(ifthenelse(self.River == 1, 0.0, self.TCSagg),self.TCSagg) - self.TCLagg = cover(ifthenelse(self.filter_Eros_TC == 0, 10**9, self.TCLagg),self.TCLagg) - self.TCLagg = cover(ifthenelse(self.River == 1, 0.0, self.TCLagg),self.TCLagg) - - #Eroded sediment in overland flow reaching the river system per particle class [ton/cell/timestep] - self.InLandClay = cover(ifthenelse(self.River == 1, accucapacitystate(self.TopoLdd, self.LandErodClay, self.TCClay), 0.0), 0.0) - self.InLandSilt = cover(ifthenelse(self.River == 1, accucapacitystate(self.TopoLdd, self.LandErodSilt, self.TCSilt), 0.0), 0.0) - self.InLandSand = cover(ifthenelse(self.River == 1, accucapacitystate(self.TopoLdd, self.LandErodSand, self.TCSand), 0.0), 0.0) - self.InLandSagg = cover(ifthenelse(self.River == 1, accucapacitystate(self.TopoLdd, self.LandErodSagg, self.TCSagg), 0.0), 0.0) - self.InLandLagg = cover(ifthenelse(self.River == 1, accucapacitystate(self.TopoLdd, self.LandErodLagg, self.TCLagg), 0.0), 0.0) - self.InLandSed = self.InLandClay + self.InLandSilt + self.InLandSand + self.InLandSagg + self.InLandLagg #total - - #Sediment input reaching the river system per particle class and per HYPE subcatchment [kg/ha/timestep] - self.ClayCatch = areatotal(self.InLandClay, self.HYPEcatch) * 1000 / (self.cellareaKm *100) - self.SiltCatch = areatotal(self.InLandSilt, self.HYPEcatch) * 1000 / (self.cellareaKm *100) - self.SandCatch = areatotal(self.InLandSand, self.HYPEcatch) * 1000 / (self.cellareaKm *100) - self.SaggCatch = areatotal(self.InLandSagg, self.HYPEcatch) * 1000 / (self.cellareaKm *100) - self.LaggCatch = areatotal(self.InLandLagg, self.HYPEcatch) * 1000 / (self.cellareaKm *100) - self.HYPESedCatch = self.ClayCatch + self.SiltCatch + self.SandCatch + self.SaggCatch + self.LaggCatch #total - - ########################################################################## - # River transport and processes ############################# - ########################################################################## - - if self.RunRiverModel == 1: - #River sediment loads are separated into different particle class. - #Clay, silt and sand can both come from land, resuspension or river channel erosion. - #Small and large aggregates only come from land erosion or resuspension. - #Gravel only comes from resuspension or river channel erosion. - - ''' Initial concentration and input from land and upstream river cells ''' - #Save the loads from the previous time step that remained in the cell - self.OldSedLoad = self.SedLoad - self.OldClayLoad = self.ClayLoad - self.OldSiltLoad = self.SiltLoad - self.OldSandLoad = self.SandLoad - self.OldSaggLoad = self.SaggLoad - self.OldLaggLoad = self.LaggLoad - self.OldGravLoad = self.GravLoad - - #Input concentration including the old load, the incoming load from upstream river cells and the incoming load from land erosion - self.InSedLoad = self.OldSedLoad + upstream(self.TopoLdd, self.OutSedLoad) + self.InLandSed - self.InClayLoad = self.OldClayLoad + upstream(self.TopoLdd, self.OutClayLoad) + self.InLandClay - self.InSiltLoad = self.OldSiltLoad + upstream(self.TopoLdd, self.OutSiltLoad) + self.InLandSilt - self.InSandLoad = self.OldSandLoad + upstream(self.TopoLdd, self.OutSandLoad) + self.InLandSand - self.InSaggLoad = self.OldSaggLoad + upstream(self.TopoLdd, self.OutSaggLoad) + self.InLandSagg - self.InLaggLoad = self.OldLaggLoad + upstream(self.TopoLdd, self.OutLaggLoad) + self.InLandLagg - self.InGravLoad = self.OldGravLoad + upstream(self.TopoLdd, self.OutGravLoad) - - ''' River erosion ''' - #Hydraulic radius of the river [m] (rectangular channel) - self.HydRad = self.WaterLevel*self.RiverWidth / (self.RiverWidth + 2*self.WaterLevel) - - #Transport capacity - #Engelund and Hansen transport formula - if self.RivTransportMethod == 1 : - #Shields parameter - self.ThetaShields = 1000 * 9.81 * self.HydRad * self.RiverSlope / ((2650-1000) * 9.81 * self.D50River/1000) - self.vmean = ifthenelse(pcrand(self.WaterLevel > 0, self.River ==1), self.SurfaceRunoff/(self.RiverWidth*self.WaterLevel), scalar(0.0)) - self.vshear = (9.81 * self.HydRad * self.RiverSlope)**(0.5) - # self.Cw = 0.05 * (2650/(2650-1000)) * self.vmean * self.RiverSlope \ - # / ((2650-1000)/1000*9.81*self.D50River/1000)**(0.5) * self.ThetaShields**(0.5) # sediment concentration by weight - self.Cw = min( 1.0, ifthenelse(pcrand(self.HydRad > 0, self.River ==1), 2.65 * 0.05 * self.vmean * self.vshear**3 \ - / ((2.65-1)**2 * 9.81**2 * self.D50River * self.HydRad), scalar(0.0))) #concentration by weight - self.MaxSedLoad = max(0.0, self.Cw / (self.Cw + (1-self.Cw)*2.65) * 2.65) #[tons/m3] - self.MaxSedLoad = self.MaxSedLoad * (self.WaterLevel * self.RiverWidth * self.DCL + self.SurfaceRunoff*self.timestepsecs) #[ton] - - #Simplified Bagnold transport formula - if self.RivTransportMethod == 2 : - self.MaxSedLoad = self.cBagnold * (self.SurfaceRunoff / (self.WaterLevel * self.RiverWidth))**self.expBagnold #[ton/m3] - self.MaxSedLoad = self.MaxSedLoad * (self.WaterLevel * self.RiverWidth * self.DCL + self.SurfaceRunoff*self.timestepsecs) #[ton] - - #Kodatie transport formula - if self.RivTransportMethod == 3 : - self.vmean = ifthenelse(pcrand(self.WaterLevel > 0, self.River ==1), self.SurfaceRunoff/(self.RiverWidth*self.WaterLevel), scalar(0.0)) - self.MaxSedLoad = (self.aK * self.vmean**(self.bK) * self.WaterLevel**(self.cK) * self.RiverSlope**(self.dK)) * (self.RiverWidth) #[tons] - - #Yang transport formula - if self.RivTransportMethod == 4 : - self.wsRiv = 411 * self.D50River**2 / 3600 - self.vshear = (9.81 * self.HydRad * self.RiverSlope)**(0.5) - self.var1 = self.vshear * self.D50River/1000 / (1.16 * 10**(-6)) - self.var2 = self.wsRiv * self.D50River/1000 / (1.16 * 10**(-6)) - self.vcr = ifthenelse(self.var1 >= 70, 2.05 * self.wsRiv, self.wsRiv *(2.5 / (log10(self.var1)-0.06) + 0.66)) - - #Sand equation - self.logCppm = 5.435 - 0.286*log10(self.var2) - 0.457*log10(self.vshear/self.wsRiv) + \ - (1.799 - 0.409*log10(self.var2) - 0.314*log10(self.vshear/self.wsRiv)) * \ - log10((self.SurfaceRunoff/(self.RiverWidth*self.WaterLevel) - self.vcr) * self.RiverSlope/self.wsRiv) - #Gravel equation - self.logCppm = ifthenelse(self.D50River < 2.0, self.logCppm, 6.681 - 0.633*log10(self.var2) - 4.816*log10(self.vshear/self.wsRiv) + \ - (2.784 - 0.305*log10(self.var2) - 0.282*log10(self.vshear/self.wsRiv)) * \ - log10((self.SurfaceRunoff/(self.RiverWidth*self.WaterLevel) - self.vcr) * self.RiverSlope/self.wsRiv)) - self.Cw = 10**self.logCppm * 10**(-6) #sediment concentration by weight - self.MaxSedLoad = max(0.0, self.Cw / (self.Cw + (1-self.Cw)*2.65) * 2.65) #[tons/m3] - self.MaxSedLoad = self.MaxSedLoad * (self.WaterLevel * self.RiverWidth * self.DCL + self.SurfaceRunoff*self.timestepsecs) #[ton] - - #Molinas & Wu transport formula - if self.RivTransportMethod == 5 : - self.wsRiv = 411 * self.D50River**2 / 3600 - self.psi = ((2.65-1)*9.81*self.WaterLevel*self.wsRiv*(log10(self.WaterLevel/self.D50River))**2)**(0.5) - self.Cw = 1430 * (0.86 + (self.psi)**(0.5)) * (self.psi)**(1.5) / (0.016 + self.psi) * 10**(-6) #weight - self.MaxSedLoad = max(0.0, self.Cw / (self.Cw + (1-self.Cw)*2.65) * 2.65) #[tons/m3] - self.MaxSedLoad = self.MaxSedLoad * (self.WaterLevel * self.RiverWidth * self.DCL + self.SurfaceRunoff*self.timestepsecs) #[ton] - - - self.MaxSedLoad = ifthen(self.subcatch > 0, cover(self.MaxSedLoad, 0.0)) - - - #Repartition of the effective shear stress between the bank and the bed from Knight et al. 1984 - self.SFBank = ifthenelse(pcrand(self.WaterLevel > 0, self.River == 1), exp(-3.230 * log10(self.RiverWidth / self.WaterLevel + 3) + 6.146), scalar(0.0)) #[%] - #Effective shear stress on river bed and banks [N/m2] - self.TEffBank = ifthenelse(pcrand(self.WaterLevel > 0, self.River == 1), 1000 * 9.81 * self.HydRad * self.RiverSlope * self.SFBank/100 \ - * (1 + self.RiverWidth / (2*self.WaterLevel)), scalar(0.0)) - self.TEffBed = 1000 * 9.81 * self.HydRad * self.RiverSlope * (1 - self.SFBank/100) * (1 + 2*self.WaterLevel/self.RiverWidth) - #Potential erosion rates of the bed and bank [t/cell/timestep] (assuming only one bank is eroding) - self.ERBank = max(0.0, self.KdBank * (self.TEffBank - self.TCrBank) * (self.DCL * self.WaterLevel) * 1.4 * self.timestepsecs) #1.4 is bank default bulk density - self.ERBed = max(0.0, self.KdBed * (self.TEffBed - self.TCrBed) * (self.DCL * self.RiverWidth) * 1.5 * self.timestepsecs) #1.5 is bed default bulk density - #Relative potential erosion rates of the bed and bank [-] - self.RTEBank = ifthenelse(self.ERBank + self.ERBed > 0.0, self.ERBank/(self.ERBank+self.ERBed), 0.0) - self.RTEBed = 1.0 - self.RTEBank - - #Excess transport capacity [ton/cell/timestep] - #Erosion only if the load is below the transport capacity of the flow - self.SedEx = max(self.MaxSedLoad - self.InSedLoad, 0.0) - #Bed and bank are eroded after the previously deposited material - self.EffSedEx = max(self.SedEx - self.RivStoreSed, 0.0) - - #Bank erosion [ton/cell/timestep] - self.BankSedLoad = ifthenelse(self.EffSedEx == 0, 0.0, ifthenelse(self.EffSedEx*self.RTEBank <= self.ERBank, - self.EffSedEx*self.RTEBank, self.ERBank)) - self.BankClay = self.FracClayRiv * self.BankSedLoad - self.BankSilt = self.FracSiltRiv * self.BankSedLoad - self.BankSand = self.FracSandRiv * self.BankSedLoad - self.BankGrav = self.FracGravRiv * self.BankSedLoad - - #Bed erosion [ton/cell/timestep] - self.BedSedLoad = ifthenelse(self.EffSedEx == 0, 0.0, ifthenelse(self.EffSedEx*self.RTEBed <= self.ERBed, - self.EffSedEx*self.RTEBed, self.ERBed)) - self.BedClay = self.FracClayRiv * self.BedSedLoad - self.BedSilt = self.FracSiltRiv * self.BedSedLoad - self.BedSand = self.FracSandRiv * self.BedSedLoad - self.BedGrav = self.FracGravRiv * self.BedSedLoad - - #Erosion/degradation of the previously deposited sediment (from clay to gravel) [ton/cell/timestep] - self.DegStoreClay = ifthenelse(self.RivStoreClay >= self.SedEx, self.SedEx, self.RivStoreClay) - self.RivStoreClay = self.RivStoreClay - self.DegStoreClay #update store - self.SedEx = max(self.SedEx - self.DegStoreClay, 0.0) #update amount of sediment that need to be degraded - - self.DegStoreSilt = ifthenelse(self.RivStoreSilt >= self.SedEx, self.SedEx, self.RivStoreSilt) - self.RivStoreSilt = self.RivStoreSilt - self.DegStoreSilt - self.SedEx = max(self.SedEx - self.DegStoreSilt, 0.0) - - self.DegStoreSagg = ifthenelse(self.RivStoreSagg >= self.SedEx, self.SedEx, self.RivStoreSagg) - self.RivStoreSagg = self.RivStoreSagg - self.DegStoreSagg - self.SedEx = max(self.SedEx - self.DegStoreSagg, 0.0) - - self.DegStoreSand = ifthenelse(self.RivStoreSand >= self.SedEx, self.SedEx, self.RivStoreSand) - self.RivStoreSand = self.RivStoreSand - self.DegStoreSand - self.SedEx = max(self.SedEx - self.DegStoreSand, 0.0) - - self.DegStoreLagg = ifthenelse(self.RivStoreLagg >= self.SedEx, self.SedEx, self.RivStoreLagg) - self.RivStoreLagg = self.RivStoreLagg - self.DegStoreLagg - self.SedEx = max(self.SedEx - self.DegStoreLagg, 0.0) - - self.DegStoreGrav = ifthenelse(self.RivStoreGrav >= self.SedEx, self.SedEx, self.RivStoreGrav) - self.RivStoreGrav = self.RivStoreGrav - self.DegStoreGrav - self.SedEx = max(self.SedEx - self.DegStoreGrav, 0.0) - - self.DegStoreSed = self.DegStoreClay + self.DegStoreSilt + self.DegStoreSand + self.DegStoreSagg + self.DegStoreLagg + self.DegStoreGrav - self.RivStoreSed = self.RivStoreSed - self.DegStoreSed - - - #Sum all erosion sources per particle class - self.RivErodSed = self.BankSedLoad + self.BedSedLoad + self.DegStoreSed - self.RivErodClay = self.BankClay + self.BedClay + self.DegStoreClay - self.RivErodSilt = self.BankSilt + self.BedSilt + self.DegStoreSilt - self.RivErodSand = self.BankSand + self.BedSand + self.DegStoreSand - self.RivErodSagg = self.DegStoreSagg - self.RivErodLagg = self.DegStoreLagg - self.RivErodGrav = self.BankGrav + self.BedGrav + self.DegStoreGrav - - #Assume that there is no erosion/resuspension in reservoir - if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: - self.RivErodSed = self.filter_Eros_TC * self.RivErodSed - self.RivErodClay = self.filter_Eros_TC * self.RivErodClay - self.RivErodSilt = self.filter_Eros_TC * self.RivErodSilt - self.RivErodSand = self.filter_Eros_TC * self.RivErodSand - self.RivErodSagg = self.filter_Eros_TC * self.RivErodSagg - self.RivErodLagg = self.filter_Eros_TC * self.RivErodLagg - self.RivErodGrav = self.filter_Eros_TC * self.RivErodGrav - - - ''' Deposition/settling ''' - #If transport capacity is exceeded, the excess is deposited. - #Else classic deposition with Einstein formula. - self.TCEx = self.MaxSedLoad - self.InSedLoad - self.FracDepEx = ifthenelse(self.TCEx < 0, 1 - self.MaxSedLoad / self.InSedLoad, scalar(0.0)) - - #Fractions of deposited particles in river cells from the Einstein formula [-] - #Particle fall velocity [m/s] from Stokes - self.x = ifthenelse(pcrand(self.SurfaceRunoff > 0, self.River ==1), 1.055*self.DCL / (self.SurfaceRunoff / self.RiverWidth), scalar(0.0)) - self.x = ifthen(self.subcatch > 0, cover(self.x, 0.0)) - - self.FracDepClay = ifthenelse(self.TCEx >= 0, min( 1.0, 1 - 1 / exp(self.x * (411 * 0.002**2 / 3600))), self.FracDepEx) - self.FracDepSilt = ifthenelse(self.TCEx >= 0, min( 1.0, 1 - 1 / exp(self.x * (411 * 0.010**2 / 3600))), self.FracDepEx) - self.FracDepSand = ifthenelse(self.TCEx >= 0, min( 1.0, 1 - 1 / exp(self.x * (411 * 0.200**2 / 3600))), self.FracDepEx) - self.FracDepSagg = ifthenelse(self.TCEx >= 0, min( 1.0, 1 - 1 / exp(self.x * (411 * 0.030**2 / 3600))), self.FracDepEx) - self.FracDepLagg = ifthenelse(self.TCEx >= 0, min( 1.0, 1 - 1 / exp(self.x * (411 * 0.500**2 / 3600))), self.FracDepEx) - self.FracDepGrav = ifthenelse(self.TCEx >= 0, min( 1.0, 1 - 1 / exp(self.x * (411 * 2.000**2 / 3600))), self.FracDepEx) - - #Sediment deposited in the channel [ton/cell/timestep] - self.DepClay = self.FracDepClay * (self.InClayLoad + self.RivErodClay) - self.DepSilt = self.FracDepSilt * (self.InSiltLoad + self.RivErodSilt) - self.DepSand = self.FracDepSand * (self.InSandLoad + self.RivErodSand) - self.DepSagg = self.FracDepSagg * (self.InSaggLoad + self.RivErodSagg) - self.DepLagg = self.FracDepLagg * (self.InLaggLoad + self.RivErodLagg) - self.DepGrav = self.FracDepGrav * (self.InGravLoad + self.RivErodGrav) - self.DepSedLoad = self.DepClay + self.DepSilt + self.DepSand + self.DepSagg + self.DepLagg + self.DepGrav - - #Assume that deposition happens only on reservoir output cells where the reservoir mass balance takes place - if (self.nrresSimple + self.nrresComplex) > 0: - self.DepClay = self.filter_Eros_TC * self.DepClay - self.DepSilt = self.filter_Eros_TC * self.DepSilt - self.DepSand = self.filter_Eros_TC * self.DepSand - self.DepSagg = self.filter_Eros_TC * self.DepSagg - self.DepLagg = self.filter_Eros_TC * self.DepLagg - self.DepGrav = self.filter_Eros_TC * self.DepGrav - self.DepSedLoad = self.filter_Eros_TC * self.DepSedLoad - - #Deposition in reservoirs from Camp 1945 - if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: - self.VcRes = ifthenelse(self.ReserVoirLocs > 0, self.SurfaceRunoff / self.ResArea, 0.0) - self.DepClay = cover(ifthen(self.ReserVoirLocs > 0, self.InClayLoad * min(1.0, (411/3600*0.002**2)/self.VcRes)), self.DepClay) - self.DepSilt = cover(ifthen(self.ReserVoirLocs > 0, self.InSiltLoad * min(1.0, (411/3600*0.010**2)/self.VcRes)), self.DepSilt) - self.DepSand = cover(ifthen(self.ReserVoirLocs > 0, self.InSandLoad * min(1.0, (411/3600*0.200**2)/self.VcRes)), self.DepSand) - self.DepSagg = cover(ifthen(self.ReserVoirLocs > 0, self.InSaggLoad * min(1.0, (411/3600*0.030**2)/self.VcRes)), self.DepSagg) - self.DepLagg = cover(ifthen(self.ReserVoirLocs > 0, self.InLaggLoad * min(1.0, (411/3600*0.500**2)/self.VcRes)), self.DepLagg) - self.DepGrav = cover(ifthen(self.ReserVoirLocs > 0, self.InGravLoad * min(1.0, (411/3600*2.000**2)/self.VcRes)), self.DepGrav) - self.DepSedLoad = cover(ifthen(self.ReserVoirLocs > 0, self.DepClay + self.DepSilt + self.DepSand + self.DepSagg + self.DepLagg + self.DepGrav), self.DepSedLoad) - - #Update the river deposited sediment storage - self.RivStoreSed = self.RivStoreSed + self.DepSedLoad - self.RivStoreClay = self.RivStoreClay + self.DepClay - self.RivStoreSilt = self.RivStoreSilt + self.DepSilt - self.RivStoreSand = self.RivStoreSand + self.DepSand - self.RivStoreSagg = self.RivStoreSagg + self.DepSagg - self.RivStoreLagg = self.RivStoreLagg + self.DepLagg - self.RivStoreGrav = self.RivStoreGrav + self.DepGrav - - - ''' Output loads ''' - #Sediment transported out of the cell during the timestep [ton/cell/timestep] - # 0 in case all sediment are deposited in the cell - # Reduce the fraction so that there is still some sediment staying in the river cell - self.OutFracWat = min(self.SurfaceRunoff*self.timestepsecs / (self.WaterLevel*self.RiverWidth*self.DCL + self.SurfaceRunoff*self.timestepsecs), 0.999999) - - self.OutSedLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InSedLoad+self.RivErodSed-self.DepSedLoad) * self.OutFracWat), 0.0)) - self.OutClayLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InClayLoad+self.RivErodClay-self.DepClay) * self.OutFracWat), 0.0)) - self.OutSiltLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InSiltLoad+self.RivErodSilt-self.DepSilt) * self.OutFracWat), 0.0)) - self.OutSandLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InSandLoad+self.RivErodSand-self.DepSand) * self.OutFracWat), 0.0)) - self.OutSaggLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InSaggLoad+self.RivErodSagg-self.DepSagg) * self.OutFracWat), 0.0)) - self.OutLaggLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InLaggLoad+self.RivErodLagg-self.DepLagg) * self.OutFracWat), 0.0)) - self.OutGravLoad = ifthen(self.subcatch > 0, cover(max(0.0, (self.InGravLoad+self.RivErodGrav-self.DepGrav) * self.OutFracWat), 0.0)) - - - ''' Mass balance, sediment concentration in each river cell ''' - #Sediment load [ton/cell/timestep] - #0 in case all sediment are deposited in the cell - self.SedLoad = max(0.0, self.InSedLoad + self.RivErodSed - self.DepSedLoad - self.OutSedLoad) - self.ClayLoad = max(0.0, self.InClayLoad + self.RivErodClay - self.DepClay - self.OutClayLoad) - self.SiltLoad = max(0.0, self.InSiltLoad + self.RivErodSilt - self.DepSilt - self.OutSiltLoad) - self.SandLoad = max(0.0, self.InSandLoad + self.RivErodSand - self.DepSand - self.OutSandLoad) - self.SaggLoad = max(0.0, self.InSaggLoad + self.RivErodSagg - self.DepSagg - self.OutSaggLoad) - self.LaggLoad = max(0.0, self.InLaggLoad + self.RivErodLagg - self.DepLagg - self.OutLaggLoad) - self.GravLoad = max(0.0, self.InGravLoad + self.RivErodGrav - self.DepGrav - self.OutGravLoad) - - - #Sediment concentration [mg/L] - #The suspended load is composed of clay and silt. - #The bed load is composed of sand, small and large aggregates from overland flow and gravel coming from river bed erosion. - #Conversion from load [ton] to concentration for rivers [mg/L] - #self.ToConc = 10**(6) / (self.RiverWidth * self.WaterLevel * self.DCL) - self.ToConc = 10**(6) / (self.SurfaceRunoff * self.timestepsecs) - self.SedConc = self.OutSedLoad * self.ToConc - self.SSConc = (self.OutClayLoad + self.OutSiltLoad) * self.ToConc - self.BedConc = (self.OutSandLoad + self.OutSaggLoad + self.OutLaggLoad + self.OutGravLoad) * self.ToConc - self.MaxSedConc = self.MaxSedLoad * self.ToConc - - - ''' Check / summary maps ''' - #Checking processes contribution to sediment budget - self.PrecipitationCatch = areatotal(self.Precipitation, self.HYPEcatch) - #self.InLandSedCatch = areatotal(self.InLandSed, self.HYPEcatch) - #self.MaxSedLoadCatch = areatotal(self.MaxSedLoad, self.HYPEcatch) - #self.SedExCatch = areatotal(self.SedEx, self.HYPEcatch) - #self.BankSedLoadCatch = areatotal(self.BankSedLoad, self.HYPEcatch) - #self.BedSedLoadCatch = areatotal(self.BedSedLoad, self.HYPEcatch) - #self.DegStoreSedCatch = areatotal(self.DegStoreSed, self.HYPEcatch) - #self.DepSedLoadCatch = areatotal(self.DepSedLoad, self.HYPEcatch) - #self.OutSedLoadCatch = areatotal(self.OutSedLoad, self.HYPEcatch) - #self.SedLoadCatch = areatotal(self.SedLoad, self.HYPEcatch) - - #Area specific runoff - #annual runoff / upstream (area) - #L/km2/s - #values between 5 and 25 for the Rhine - self.SpecRun = self.SurfaceRunoff * 1000 / accuflux(self.TopoLdd, self.cellareaKm) #[L/km2/s] - - - + + ########################################################################## + # Inland sediment routing model ############################# + ########################################################################## + + # If the river model is run, use Yalin's equation with particle differentiation. + # If just the soil loss model is run, use Govers equation without particle differentiation + + """ Transport of overland flow with no particle differenciation using Govers equation""" + + if self.LandTransportMethod == 2 or self.LandTransportMethod == 3: + if self.LandTransportMethod == 2: + # Unit stream power + self.velocity = cover( + ifthenelse( + self.WaterLevel > 0, + self.OvRun / (self.reallength * self.WaterLevel), + 0.0, + ), + 0.0, + ) # [m/s] + self.omega = 10 * self.sinSlope * 100 * self.velocity # [cm/s] + # self.omega = self.sinSlope * 100 * self.velocity #[cm/s] + + # Transport capacity from Govers, 1990 + self.TCf = ifthenelse( + self.omega > 0.4, + self.cGovers * (self.omega - 0.4) ** self.nGovers * 2650, + 0.0, + ) # [kg/m3] + # self.TC = self.TCf / (1 - self.TCf/2650) #[kg/m3] + # self.TC = max(self.TC, 2650) + self.TC = ( + self.TCf * self.OvRun * self.timestepsecs * 10 ** (-3) + ) # [ton/cell/timestep] + # Remove nodata values + self.TC = ifthen(self.subcatch > 0, cover(self.TC, 0.0)) + # Assume that eroded soil on lake cells all reach the river cells of the reservoir + if ( + self.nrresSimple + self.nrresComplex + ) > 0 and self.filterResArea == 0: + self.TC = ifthenelse( + pcrand(self.filter_Eros_TC == 0, self.River == 0), + 10 ** 9, + self.TC, + ) + + elif self.LandTransportMethod == 3: + # Transport capacity from Yalin + self.OvLevel = cover( + ifthenelse(self.River == 1, 0.0, self.WaterLevel), self.WaterLevel + ) # [m] + self.delta = max( + self.OvLevel + * self.sinSlope + / (self.D50 * 10 ** (-3) * (2650 / 1000 - 1)) + / 0.06 + - 1, + 0.0, + ) + self.TC = ( + self.reallength + / self.OvRun + * (2650 - 1000) + * self.D50 + * 10 ** (-3) + * (9.81 * self.OvLevel * self.sinSlope) + * 0.635 + * self.delta + * ( + 1 + - ln(1 + self.delta * 2.45 / (2650 / 1000) ** 0.4 * 0.06 ** 0.5) + / self.delta + * 2.45 + / (2650 / 1000) ** 0.4 + * 0.06 ** 0.5 + ) + ) # [kg/m3] + self.TC = ifthen( + self.subcatch > 0, + cover(self.TC * self.OvRun * self.timestepsecs * 10 ** (-3), 0.0), + ) # [ton/cell/timestep] + # Assume that eroded soil on lake cells all reach the river cells of the reservoir + if ( + self.nrresSimple + self.nrresComplex + ) > 0 and self.filterResArea == 0: + self.TC = ifthenelse( + pcrand(self.filter_Eros_TC == 0, self.River == 0), + 10 ** 9, + self.TC, + ) + + # To get total sediment input from land into the river systems, river cells transport all sediment to the output (huge TC) + self.TCRiv = cover(ifthenelse(self.River == 1, 10 ** 9, self.TC), self.TC) + # Transported sediment over the land + self.SedFlux = accucapacityflux( + self.TopoLdd, self.SoilLoss, self.TCRiv + ) # [ton/cell/tinestep] + # Deposited sediment over the land + self.SedDep = accucapacitystate(self.TopoLdd, self.SoilLoss, self.TCRiv) + + # Sediment amount reaching each river cell ''' + self.SedDep2 = accucapacitystate(self.TopoLdd, self.SoilLoss, self.TC) + # Remove inland deposition + self.OvSed = cover(ifthenelse(self.River == 1, self.SedDep2, 0.0), 0.0) + + # Sum the results per subcatchment (input for D-WAQ) [kg/ha/timestep] + self.HYPEOvSedCatch = areatotal( + self.OvSed, self.HYPEcatch + ) # [ton/cell/timestep] + self.HYPEOvSedCatch = ( + self.HYPEOvSedCatch * 1000 / (self.cellareaKm * 100) + ) # [kg/ha/timestep] + + """ Transport of overland flow with particle differenciation using Yalin equation""" + else: + # Determine the eroded amount of clay/silt/sand/aggregates on each cell [ton/cell/timestep] + self.LandErodClay = self.SoilLoss * self.FracClay + self.LandErodSilt = self.SoilLoss * self.FracSilt + self.LandErodSand = self.SoilLoss * self.FracSand + self.LandErodSagg = self.SoilLoss * self.FracSagg + self.LandErodLagg = self.SoilLoss * self.FracLagg + + # Water level of overland flow + self.OvLevel = cover( + ifthenelse(self.River == 1, 0.0, self.WaterLevel), self.WaterLevel + ) # [m] + + # Delta parameter of Yalin for each particle class + self.DClay = max( + self.OvLevel + * self.sinSlope + / (2 * 10 ** (-6) * (2650 / 1000 - 1)) + / 0.06 + - 1, + 0.0, + ) + self.DSilt = max( + self.OvLevel + * self.sinSlope + / (10 * 10 ** (-6) * (2650 / 1000 - 1)) + / 0.06 + - 1, + 0.0, + ) + self.DSand = max( + self.OvLevel + * self.sinSlope + / (200 * 10 ** (-6) * (2650 / 1000 - 1)) + / 0.06 + - 1, + 0.0, + ) + self.DSagg = max( + self.OvLevel + * self.sinSlope + / (30 * 10 ** (-6) * (2650 / 1000 - 1)) + / 0.06 + - 1, + 0.0, + ) + self.DLagg = max( + self.OvLevel + * self.sinSlope + / (500 * 10 ** (-6) * (2650 / 1000 - 1)) + / 0.06 + - 1, + 0.0, + ) + + # Total transportability + self.Dtot = self.DClay + self.DSilt + self.DSand + self.DSagg + self.DLagg + + # Yalin Transport capacity of overland flow for each particle class + self.TCClay = ( + self.reallength + / self.OvRun + * (2650 - 1000) + * 2 + * 10 ** (-6) + * (9.81 * self.OvLevel * self.sinSlope) + * self.DClay + / self.Dtot + * 0.635 + * self.DClay + * ( + 1 + - ln(1 + self.DClay * 2.45 / (2650 / 1000) ** 0.4 * 0.06 ** 0.5) + / self.DClay + * 2.45 + / (2650 / 1000) ** 0.4 + * 0.06 ** 0.5 + ) + ) # [kg/m3] + self.TCClay = ifthen( + self.subcatch > 0, + cover(self.TCClay * self.OvRun * self.timestepsecs * 10 ** (-3), 0.0), + ) # [ton/cell/timestep] + + self.TCSilt = ( + self.reallength + / self.OvRun + * (2650 - 1000) + * 10 + * 10 ** (-6) + * (9.81 * self.OvLevel * self.sinSlope) + * self.DSilt + / self.Dtot + * 0.635 + * self.DSilt + * ( + 1 + - ln(1 + self.DSilt * 2.45 / (2650 / 1000) ** 0.4 * 0.06 ** 0.5) + / self.DSilt + * 2.45 + / (2650 / 1000) ** 0.4 + * 0.06 ** 0.5 + ) + ) # [kg/m3] + self.TCSilt = ifthen( + self.subcatch > 0, + cover(self.TCSilt * self.OvRun * self.timestepsecs * 10 ** (-3), 0.0), + ) # [ton/cell/timestep] + + self.TCSand = ( + self.reallength + / self.OvRun + * (2650 - 1000) + * 200 + * 10 ** (-6) + * (9.81 * self.OvLevel * self.sinSlope) + * self.DSand + / self.Dtot + * 0.635 + * self.DSand + * ( + 1 + - ln(1 + self.DSand * 2.45 / (2650 / 1000) ** 0.4 * 0.06 ** 0.5) + / self.DSand + * 2.45 + / (2650 / 1000) ** 0.4 + * 0.06 ** 0.5 + ) + ) # [kg/m3] + self.TCSand = ifthen( + self.subcatch > 0, + cover(self.TCSand * self.OvRun * self.timestepsecs * 10 ** (-3), 0.0), + ) # [ton/cell/timestep] + + self.TCSagg = ( + self.reallength + / self.OvRun + * (2650 - 1000) + * 30 + * 10 ** (-6) + * (9.81 * self.OvLevel * self.sinSlope) + * self.DSagg + / self.Dtot + * 0.635 + * self.DSagg + * ( + 1 + - ln(1 + self.DSagg * 2.45 / (2650 / 1000) ** 0.4 * 0.06 ** 0.5) + / self.DSagg + * 2.45 + / (2650 / 1000) ** 0.4 + * 0.06 ** 0.5 + ) + ) # [kg/m3] + self.TCSagg = ifthen( + self.subcatch > 0, + cover(self.TCSagg * self.OvRun * self.timestepsecs * 10 ** (-3), 0.0), + ) # [ton/cell/timestep] + + self.TCLagg = ( + self.reallength + / self.OvRun + * (2650 - 1000) + * 500 + * 10 ** (-6) + * (9.81 * self.OvLevel * self.sinSlope) + * self.DLagg + / self.Dtot + * 0.635 + * self.DLagg + * ( + 1 + - ln(1 + self.DLagg * 2.45 / (2650 / 1000) ** 0.4 * 0.06 ** 0.5) + / self.DLagg + * 2.45 + / (2650 / 1000) ** 0.4 + * 0.06 ** 0.5 + ) + ) # [kg/m3] + self.TCLagg = ifthen( + self.subcatch > 0, + cover(self.TCLagg * self.OvRun * self.timestepsecs * 10 ** (-3), 0.0), + ) # [ton/cell/timestep] + + # Assume that eroded soil on lake cells all reach the river cells of the reservoir + if (self.nrresSimple + self.nrresComplex) > 0: + self.TCClay = cover( + ifthenelse(self.filter_Eros_TC == 0, 10 ** 9, self.TCClay), + self.TCClay, + ) + self.TCClay = cover( + ifthenelse(self.River == 1, 0.0, self.TCClay), self.TCClay + ) + self.TCSilt = cover( + ifthenelse(self.filter_Eros_TC == 0, 10 ** 9, self.TCSilt), + self.TCSilt, + ) + self.TCSilt = cover( + ifthenelse(self.River == 1, 0.0, self.TCSilt), self.TCSilt + ) + self.TCSand = cover( + ifthenelse(self.filter_Eros_TC == 0, 10 ** 9, self.TCSand), + self.TCSand, + ) + self.TCSand = cover( + ifthenelse(self.River == 1, 0.0, self.TCSand), self.TCSand + ) + self.TCSagg = cover( + ifthenelse(self.filter_Eros_TC == 0, 10 ** 9, self.TCSagg), + self.TCSagg, + ) + self.TCSagg = cover( + ifthenelse(self.River == 1, 0.0, self.TCSagg), self.TCSagg + ) + self.TCLagg = cover( + ifthenelse(self.filter_Eros_TC == 0, 10 ** 9, self.TCLagg), + self.TCLagg, + ) + self.TCLagg = cover( + ifthenelse(self.River == 1, 0.0, self.TCLagg), self.TCLagg + ) + + # Eroded sediment in overland flow reaching the river system per particle class [ton/cell/timestep] + self.InLandClay = cover( + ifthenelse( + self.River == 1, + accucapacitystate(self.TopoLdd, self.LandErodClay, self.TCClay), + 0.0, + ), + 0.0, + ) + self.InLandSilt = cover( + ifthenelse( + self.River == 1, + accucapacitystate(self.TopoLdd, self.LandErodSilt, self.TCSilt), + 0.0, + ), + 0.0, + ) + self.InLandSand = cover( + ifthenelse( + self.River == 1, + accucapacitystate(self.TopoLdd, self.LandErodSand, self.TCSand), + 0.0, + ), + 0.0, + ) + self.InLandSagg = cover( + ifthenelse( + self.River == 1, + accucapacitystate(self.TopoLdd, self.LandErodSagg, self.TCSagg), + 0.0, + ), + 0.0, + ) + self.InLandLagg = cover( + ifthenelse( + self.River == 1, + accucapacitystate(self.TopoLdd, self.LandErodLagg, self.TCLagg), + 0.0, + ), + 0.0, + ) + self.InLandSed = ( + self.InLandClay + + self.InLandSilt + + self.InLandSand + + self.InLandSagg + + self.InLandLagg + ) # total + + # Sediment input reaching the river system per particle class and per HYPE subcatchment [kg/ha/timestep] + self.ClayCatch = ( + areatotal(self.InLandClay, self.HYPEcatch) + * 1000 + / (self.cellareaKm * 100) + ) + self.SiltCatch = ( + areatotal(self.InLandSilt, self.HYPEcatch) + * 1000 + / (self.cellareaKm * 100) + ) + self.SandCatch = ( + areatotal(self.InLandSand, self.HYPEcatch) + * 1000 + / (self.cellareaKm * 100) + ) + self.SaggCatch = ( + areatotal(self.InLandSagg, self.HYPEcatch) + * 1000 + / (self.cellareaKm * 100) + ) + self.LaggCatch = ( + areatotal(self.InLandLagg, self.HYPEcatch) + * 1000 + / (self.cellareaKm * 100) + ) + self.HYPESedCatch = ( + self.ClayCatch + + self.SiltCatch + + self.SandCatch + + self.SaggCatch + + self.LaggCatch + ) # total + + ########################################################################## + # River transport and processes ############################# + ########################################################################## + + if self.RunRiverModel == 1: + # River sediment loads are separated into different particle class. + # Clay, silt and sand can both come from land, resuspension or river channel erosion. + # Small and large aggregates only come from land erosion or resuspension. + # Gravel only comes from resuspension or river channel erosion. + + """ Initial concentration and input from land and upstream river cells """ + # Save the loads from the previous time step that remained in the cell + self.OldSedLoad = self.SedLoad + self.OldClayLoad = self.ClayLoad + self.OldSiltLoad = self.SiltLoad + self.OldSandLoad = self.SandLoad + self.OldSaggLoad = self.SaggLoad + self.OldLaggLoad = self.LaggLoad + self.OldGravLoad = self.GravLoad + + # Input concentration including the old load, the incoming load from upstream river cells and the incoming load from land erosion + self.InSedLoad = ( + self.OldSedLoad + + upstream(self.TopoLdd, self.OutSedLoad) + + self.InLandSed + ) + self.InClayLoad = ( + self.OldClayLoad + + upstream(self.TopoLdd, self.OutClayLoad) + + self.InLandClay + ) + self.InSiltLoad = ( + self.OldSiltLoad + + upstream(self.TopoLdd, self.OutSiltLoad) + + self.InLandSilt + ) + self.InSandLoad = ( + self.OldSandLoad + + upstream(self.TopoLdd, self.OutSandLoad) + + self.InLandSand + ) + self.InSaggLoad = ( + self.OldSaggLoad + + upstream(self.TopoLdd, self.OutSaggLoad) + + self.InLandSagg + ) + self.InLaggLoad = ( + self.OldLaggLoad + + upstream(self.TopoLdd, self.OutLaggLoad) + + self.InLandLagg + ) + self.InGravLoad = self.OldGravLoad + upstream( + self.TopoLdd, self.OutGravLoad + ) + + """ River erosion """ + # Hydraulic radius of the river [m] (rectangular channel) + self.HydRad = ( + self.WaterLevel + * self.RiverWidth + / (self.RiverWidth + 2 * self.WaterLevel) + ) + + # Transport capacity + # Engelund and Hansen transport formula + if self.RivTransportMethod == 1: + # Shields parameter + self.ThetaShields = ( + 1000 + * 9.81 + * self.HydRad + * self.RiverSlope + / ((2650 - 1000) * 9.81 * self.D50River / 1000) + ) + self.vmean = ifthenelse( + pcrand(self.WaterLevel > 0, self.River == 1), + self.SurfaceRunoff / (self.RiverWidth * self.WaterLevel), + scalar(0.0), + ) + self.vshear = (9.81 * self.HydRad * self.RiverSlope) ** (0.5) + # self.Cw = 0.05 * (2650/(2650-1000)) * self.vmean * self.RiverSlope \ + # / ((2650-1000)/1000*9.81*self.D50River/1000)**(0.5) * self.ThetaShields**(0.5) # sediment concentration by weight + self.Cw = min( + 1.0, + ifthenelse( + pcrand(self.HydRad > 0, self.River == 1), + 2.65 + * 0.05 + * self.vmean + * self.vshear ** 3 + / ((2.65 - 1) ** 2 * 9.81 ** 2 * self.D50River * self.HydRad), + scalar(0.0), + ), + ) # concentration by weight + self.MaxSedLoad = max( + 0.0, self.Cw / (self.Cw + (1 - self.Cw) * 2.65) * 2.65 + ) # [tons/m3] + self.MaxSedLoad = self.MaxSedLoad * ( + self.WaterLevel * self.RiverWidth * self.DCL + + self.SurfaceRunoff * self.timestepsecs + ) # [ton] + + # Simplified Bagnold transport formula + if self.RivTransportMethod == 2: + self.MaxSedLoad = ( + self.cBagnold + * (self.SurfaceRunoff / (self.WaterLevel * self.RiverWidth)) + ** self.expBagnold + ) # [ton/m3] + self.MaxSedLoad = self.MaxSedLoad * ( + self.WaterLevel * self.RiverWidth * self.DCL + + self.SurfaceRunoff * self.timestepsecs + ) # [ton] + + # Kodatie transport formula + if self.RivTransportMethod == 3: + self.vmean = ifthenelse( + pcrand(self.WaterLevel > 0, self.River == 1), + self.SurfaceRunoff / (self.RiverWidth * self.WaterLevel), + scalar(0.0), + ) + self.MaxSedLoad = ( + self.aK + * self.vmean ** (self.bK) + * self.WaterLevel ** (self.cK) + * self.RiverSlope ** (self.dK) + ) * ( + self.RiverWidth + ) # [tons] + + # Yang transport formula + if self.RivTransportMethod == 4: + self.wsRiv = 411 * self.D50River ** 2 / 3600 + self.vshear = (9.81 * self.HydRad * self.RiverSlope) ** (0.5) + self.var1 = self.vshear * self.D50River / 1000 / (1.16 * 10 ** (-6)) + self.var2 = self.wsRiv * self.D50River / 1000 / (1.16 * 10 ** (-6)) + self.vcr = ifthenelse( + self.var1 >= 70, + 2.05 * self.wsRiv, + self.wsRiv * (2.5 / (log10(self.var1) - 0.06) + 0.66), + ) + + # Sand equation + self.logCppm = ( + 5.435 + - 0.286 * log10(self.var2) + - 0.457 * log10(self.vshear / self.wsRiv) + + ( + 1.799 + - 0.409 * log10(self.var2) + - 0.314 * log10(self.vshear / self.wsRiv) + ) + * log10( + ( + self.SurfaceRunoff / (self.RiverWidth * self.WaterLevel) + - self.vcr + ) + * self.RiverSlope + / self.wsRiv + ) + ) + # Gravel equation + self.logCppm = ifthenelse( + self.D50River < 2.0, + self.logCppm, + 6.681 + - 0.633 * log10(self.var2) + - 4.816 * log10(self.vshear / self.wsRiv) + + ( + 2.784 + - 0.305 * log10(self.var2) + - 0.282 * log10(self.vshear / self.wsRiv) + ) + * log10( + ( + self.SurfaceRunoff / (self.RiverWidth * self.WaterLevel) + - self.vcr + ) + * self.RiverSlope + / self.wsRiv + ), + ) + self.Cw = 10 ** self.logCppm * 10 ** ( + -6 + ) # sediment concentration by weight + self.MaxSedLoad = max( + 0.0, self.Cw / (self.Cw + (1 - self.Cw) * 2.65) * 2.65 + ) # [tons/m3] + self.MaxSedLoad = self.MaxSedLoad * ( + self.WaterLevel * self.RiverWidth * self.DCL + + self.SurfaceRunoff * self.timestepsecs + ) # [ton] + + # Molinas & Wu transport formula + if self.RivTransportMethod == 5: + self.wsRiv = 411 * self.D50River ** 2 / 3600 + self.psi = ( + (2.65 - 1) + * 9.81 + * self.WaterLevel + * self.wsRiv + * (log10(self.WaterLevel / self.D50River)) ** 2 + ) ** (0.5) + self.Cw = ( + 1430 + * (0.86 + (self.psi) ** (0.5)) + * (self.psi) ** (1.5) + / (0.016 + self.psi) + * 10 ** (-6) + ) # weight + self.MaxSedLoad = max( + 0.0, self.Cw / (self.Cw + (1 - self.Cw) * 2.65) * 2.65 + ) # [tons/m3] + self.MaxSedLoad = self.MaxSedLoad * ( + self.WaterLevel * self.RiverWidth * self.DCL + + self.SurfaceRunoff * self.timestepsecs + ) # [ton] + + self.MaxSedLoad = ifthen(self.subcatch > 0, cover(self.MaxSedLoad, 0.0)) + + # Repartition of the effective shear stress between the bank and the bed from Knight et al. 1984 + self.SFBank = ifthenelse( + pcrand(self.WaterLevel > 0, self.River == 1), + exp(-3.230 * log10(self.RiverWidth / self.WaterLevel + 3) + 6.146), + scalar(0.0), + ) # [%] + # Effective shear stress on river bed and banks [N/m2] + self.TEffBank = ifthenelse( + pcrand(self.WaterLevel > 0, self.River == 1), + 1000 + * 9.81 + * self.HydRad + * self.RiverSlope + * self.SFBank + / 100 + * (1 + self.RiverWidth / (2 * self.WaterLevel)), + scalar(0.0), + ) + self.TEffBed = ( + 1000 + * 9.81 + * self.HydRad + * self.RiverSlope + * (1 - self.SFBank / 100) + * (1 + 2 * self.WaterLevel / self.RiverWidth) + ) + # Potential erosion rates of the bed and bank [t/cell/timestep] (assuming only one bank is eroding) + self.ERBank = max( + 0.0, + self.KdBank + * (self.TEffBank - self.TCrBank) + * (self.DCL * self.WaterLevel) + * 1.4 + * self.timestepsecs, + ) # 1.4 is bank default bulk density + self.ERBed = max( + 0.0, + self.KdBed + * (self.TEffBed - self.TCrBed) + * (self.DCL * self.RiverWidth) + * 1.5 + * self.timestepsecs, + ) # 1.5 is bed default bulk density + # Relative potential erosion rates of the bed and bank [-] + self.RTEBank = ifthenelse( + self.ERBank + self.ERBed > 0.0, + self.ERBank / (self.ERBank + self.ERBed), + 0.0, + ) + self.RTEBed = 1.0 - self.RTEBank + + # Excess transport capacity [ton/cell/timestep] + # Erosion only if the load is below the transport capacity of the flow + self.SedEx = max(self.MaxSedLoad - self.InSedLoad, 0.0) + # Bed and bank are eroded after the previously deposited material + self.EffSedEx = max(self.SedEx - self.RivStoreSed, 0.0) + + # Bank erosion [ton/cell/timestep] + self.BankSedLoad = ifthenelse( + self.EffSedEx == 0, + 0.0, + ifthenelse( + self.EffSedEx * self.RTEBank <= self.ERBank, + self.EffSedEx * self.RTEBank, + self.ERBank, + ), + ) + self.BankClay = self.FracClayRiv * self.BankSedLoad + self.BankSilt = self.FracSiltRiv * self.BankSedLoad + self.BankSand = self.FracSandRiv * self.BankSedLoad + self.BankGrav = self.FracGravRiv * self.BankSedLoad + + # Bed erosion [ton/cell/timestep] + self.BedSedLoad = ifthenelse( + self.EffSedEx == 0, + 0.0, + ifthenelse( + self.EffSedEx * self.RTEBed <= self.ERBed, + self.EffSedEx * self.RTEBed, + self.ERBed, + ), + ) + self.BedClay = self.FracClayRiv * self.BedSedLoad + self.BedSilt = self.FracSiltRiv * self.BedSedLoad + self.BedSand = self.FracSandRiv * self.BedSedLoad + self.BedGrav = self.FracGravRiv * self.BedSedLoad + + # Erosion/degradation of the previously deposited sediment (from clay to gravel) [ton/cell/timestep] + self.DegStoreClay = ifthenelse( + self.RivStoreClay >= self.SedEx, self.SedEx, self.RivStoreClay + ) + self.RivStoreClay = self.RivStoreClay - self.DegStoreClay # update store + self.SedEx = max( + self.SedEx - self.DegStoreClay, 0.0 + ) # update amount of sediment that need to be degraded + + self.DegStoreSilt = ifthenelse( + self.RivStoreSilt >= self.SedEx, self.SedEx, self.RivStoreSilt + ) + self.RivStoreSilt = self.RivStoreSilt - self.DegStoreSilt + self.SedEx = max(self.SedEx - self.DegStoreSilt, 0.0) + + self.DegStoreSagg = ifthenelse( + self.RivStoreSagg >= self.SedEx, self.SedEx, self.RivStoreSagg + ) + self.RivStoreSagg = self.RivStoreSagg - self.DegStoreSagg + self.SedEx = max(self.SedEx - self.DegStoreSagg, 0.0) + + self.DegStoreSand = ifthenelse( + self.RivStoreSand >= self.SedEx, self.SedEx, self.RivStoreSand + ) + self.RivStoreSand = self.RivStoreSand - self.DegStoreSand + self.SedEx = max(self.SedEx - self.DegStoreSand, 0.0) + + self.DegStoreLagg = ifthenelse( + self.RivStoreLagg >= self.SedEx, self.SedEx, self.RivStoreLagg + ) + self.RivStoreLagg = self.RivStoreLagg - self.DegStoreLagg + self.SedEx = max(self.SedEx - self.DegStoreLagg, 0.0) + + self.DegStoreGrav = ifthenelse( + self.RivStoreGrav >= self.SedEx, self.SedEx, self.RivStoreGrav + ) + self.RivStoreGrav = self.RivStoreGrav - self.DegStoreGrav + self.SedEx = max(self.SedEx - self.DegStoreGrav, 0.0) + + self.DegStoreSed = ( + self.DegStoreClay + + self.DegStoreSilt + + self.DegStoreSand + + self.DegStoreSagg + + self.DegStoreLagg + + self.DegStoreGrav + ) + self.RivStoreSed = self.RivStoreSed - self.DegStoreSed + + # Sum all erosion sources per particle class + self.RivErodSed = self.BankSedLoad + self.BedSedLoad + self.DegStoreSed + self.RivErodClay = self.BankClay + self.BedClay + self.DegStoreClay + self.RivErodSilt = self.BankSilt + self.BedSilt + self.DegStoreSilt + self.RivErodSand = self.BankSand + self.BedSand + self.DegStoreSand + self.RivErodSagg = self.DegStoreSagg + self.RivErodLagg = self.DegStoreLagg + self.RivErodGrav = self.BankGrav + self.BedGrav + self.DegStoreGrav + + # Assume that there is no erosion/resuspension in reservoir + if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: + self.RivErodSed = self.filter_Eros_TC * self.RivErodSed + self.RivErodClay = self.filter_Eros_TC * self.RivErodClay + self.RivErodSilt = self.filter_Eros_TC * self.RivErodSilt + self.RivErodSand = self.filter_Eros_TC * self.RivErodSand + self.RivErodSagg = self.filter_Eros_TC * self.RivErodSagg + self.RivErodLagg = self.filter_Eros_TC * self.RivErodLagg + self.RivErodGrav = self.filter_Eros_TC * self.RivErodGrav + + """ Deposition/settling """ + # If transport capacity is exceeded, the excess is deposited. + # Else classic deposition with Einstein formula. + self.TCEx = self.MaxSedLoad - self.InSedLoad + self.FracDepEx = ifthenelse( + self.TCEx < 0, 1 - self.MaxSedLoad / self.InSedLoad, scalar(0.0) + ) + + # Fractions of deposited particles in river cells from the Einstein formula [-] + # Particle fall velocity [m/s] from Stokes + self.x = ifthenelse( + pcrand(self.SurfaceRunoff > 0, self.River == 1), + 1.055 * self.DCL / (self.SurfaceRunoff / self.RiverWidth), + scalar(0.0), + ) + self.x = ifthen(self.subcatch > 0, cover(self.x, 0.0)) + + self.FracDepClay = ifthenelse( + self.TCEx >= 0, + min(1.0, 1 - 1 / exp(self.x * (411 * 0.002 ** 2 / 3600))), + self.FracDepEx, + ) + self.FracDepSilt = ifthenelse( + self.TCEx >= 0, + min(1.0, 1 - 1 / exp(self.x * (411 * 0.010 ** 2 / 3600))), + self.FracDepEx, + ) + self.FracDepSand = ifthenelse( + self.TCEx >= 0, + min(1.0, 1 - 1 / exp(self.x * (411 * 0.200 ** 2 / 3600))), + self.FracDepEx, + ) + self.FracDepSagg = ifthenelse( + self.TCEx >= 0, + min(1.0, 1 - 1 / exp(self.x * (411 * 0.030 ** 2 / 3600))), + self.FracDepEx, + ) + self.FracDepLagg = ifthenelse( + self.TCEx >= 0, + min(1.0, 1 - 1 / exp(self.x * (411 * 0.500 ** 2 / 3600))), + self.FracDepEx, + ) + self.FracDepGrav = ifthenelse( + self.TCEx >= 0, + min(1.0, 1 - 1 / exp(self.x * (411 * 2.000 ** 2 / 3600))), + self.FracDepEx, + ) + + # Sediment deposited in the channel [ton/cell/timestep] + self.DepClay = self.FracDepClay * (self.InClayLoad + self.RivErodClay) + self.DepSilt = self.FracDepSilt * (self.InSiltLoad + self.RivErodSilt) + self.DepSand = self.FracDepSand * (self.InSandLoad + self.RivErodSand) + self.DepSagg = self.FracDepSagg * (self.InSaggLoad + self.RivErodSagg) + self.DepLagg = self.FracDepLagg * (self.InLaggLoad + self.RivErodLagg) + self.DepGrav = self.FracDepGrav * (self.InGravLoad + self.RivErodGrav) + self.DepSedLoad = ( + self.DepClay + + self.DepSilt + + self.DepSand + + self.DepSagg + + self.DepLagg + + self.DepGrav + ) + + # Assume that deposition happens only on reservoir output cells where the reservoir mass balance takes place + if (self.nrresSimple + self.nrresComplex) > 0: + self.DepClay = self.filter_Eros_TC * self.DepClay + self.DepSilt = self.filter_Eros_TC * self.DepSilt + self.DepSand = self.filter_Eros_TC * self.DepSand + self.DepSagg = self.filter_Eros_TC * self.DepSagg + self.DepLagg = self.filter_Eros_TC * self.DepLagg + self.DepGrav = self.filter_Eros_TC * self.DepGrav + self.DepSedLoad = self.filter_Eros_TC * self.DepSedLoad + + # Deposition in reservoirs from Camp 1945 + if (self.nrresSimple + self.nrresComplex) > 0 and self.filterResArea == 0: + self.VcRes = ifthenelse( + self.ReserVoirLocs > 0, self.SurfaceRunoff / self.ResArea, 0.0 + ) + self.DepClay = cover( + ifthen( + self.ReserVoirLocs > 0, + self.InClayLoad + * min(1.0, (411 / 3600 * 0.002 ** 2) / self.VcRes), + ), + self.DepClay, + ) + self.DepSilt = cover( + ifthen( + self.ReserVoirLocs > 0, + self.InSiltLoad + * min(1.0, (411 / 3600 * 0.010 ** 2) / self.VcRes), + ), + self.DepSilt, + ) + self.DepSand = cover( + ifthen( + self.ReserVoirLocs > 0, + self.InSandLoad + * min(1.0, (411 / 3600 * 0.200 ** 2) / self.VcRes), + ), + self.DepSand, + ) + self.DepSagg = cover( + ifthen( + self.ReserVoirLocs > 0, + self.InSaggLoad + * min(1.0, (411 / 3600 * 0.030 ** 2) / self.VcRes), + ), + self.DepSagg, + ) + self.DepLagg = cover( + ifthen( + self.ReserVoirLocs > 0, + self.InLaggLoad + * min(1.0, (411 / 3600 * 0.500 ** 2) / self.VcRes), + ), + self.DepLagg, + ) + self.DepGrav = cover( + ifthen( + self.ReserVoirLocs > 0, + self.InGravLoad + * min(1.0, (411 / 3600 * 2.000 ** 2) / self.VcRes), + ), + self.DepGrav, + ) + self.DepSedLoad = cover( + ifthen( + self.ReserVoirLocs > 0, + self.DepClay + + self.DepSilt + + self.DepSand + + self.DepSagg + + self.DepLagg + + self.DepGrav, + ), + self.DepSedLoad, + ) + + # Update the river deposited sediment storage + self.RivStoreSed = self.RivStoreSed + self.DepSedLoad + self.RivStoreClay = self.RivStoreClay + self.DepClay + self.RivStoreSilt = self.RivStoreSilt + self.DepSilt + self.RivStoreSand = self.RivStoreSand + self.DepSand + self.RivStoreSagg = self.RivStoreSagg + self.DepSagg + self.RivStoreLagg = self.RivStoreLagg + self.DepLagg + self.RivStoreGrav = self.RivStoreGrav + self.DepGrav + + """ Output loads """ + # Sediment transported out of the cell during the timestep [ton/cell/timestep] + # 0 in case all sediment are deposited in the cell + # Reduce the fraction so that there is still some sediment staying in the river cell + self.OutFracWat = min( + self.SurfaceRunoff + * self.timestepsecs + / ( + self.WaterLevel * self.RiverWidth * self.DCL + + self.SurfaceRunoff * self.timestepsecs + ), + 0.999999, + ) + + self.OutSedLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InSedLoad + self.RivErodSed - self.DepSedLoad) + * self.OutFracWat, + ), + 0.0, + ), + ) + self.OutClayLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InClayLoad + self.RivErodClay - self.DepClay) + * self.OutFracWat, + ), + 0.0, + ), + ) + self.OutSiltLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InSiltLoad + self.RivErodSilt - self.DepSilt) + * self.OutFracWat, + ), + 0.0, + ), + ) + self.OutSandLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InSandLoad + self.RivErodSand - self.DepSand) + * self.OutFracWat, + ), + 0.0, + ), + ) + self.OutSaggLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InSaggLoad + self.RivErodSagg - self.DepSagg) + * self.OutFracWat, + ), + 0.0, + ), + ) + self.OutLaggLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InLaggLoad + self.RivErodLagg - self.DepLagg) + * self.OutFracWat, + ), + 0.0, + ), + ) + self.OutGravLoad = ifthen( + self.subcatch > 0, + cover( + max( + 0.0, + (self.InGravLoad + self.RivErodGrav - self.DepGrav) + * self.OutFracWat, + ), + 0.0, + ), + ) + + """ Mass balance, sediment concentration in each river cell """ + # Sediment load [ton/cell/timestep] + # 0 in case all sediment are deposited in the cell + self.SedLoad = max( + 0.0, + self.InSedLoad + self.RivErodSed - self.DepSedLoad - self.OutSedLoad, + ) + self.ClayLoad = max( + 0.0, + self.InClayLoad + self.RivErodClay - self.DepClay - self.OutClayLoad, + ) + self.SiltLoad = max( + 0.0, + self.InSiltLoad + self.RivErodSilt - self.DepSilt - self.OutSiltLoad, + ) + self.SandLoad = max( + 0.0, + self.InSandLoad + self.RivErodSand - self.DepSand - self.OutSandLoad, + ) + self.SaggLoad = max( + 0.0, + self.InSaggLoad + self.RivErodSagg - self.DepSagg - self.OutSaggLoad, + ) + self.LaggLoad = max( + 0.0, + self.InLaggLoad + self.RivErodLagg - self.DepLagg - self.OutLaggLoad, + ) + self.GravLoad = max( + 0.0, + self.InGravLoad + self.RivErodGrav - self.DepGrav - self.OutGravLoad, + ) + + # Sediment concentration [mg/L] + # The suspended load is composed of clay and silt. + # The bed load is composed of sand, small and large aggregates from overland flow and gravel coming from river bed erosion. + # Conversion from load [ton] to concentration for rivers [mg/L] + # self.ToConc = 10**(6) / (self.RiverWidth * self.WaterLevel * self.DCL) + self.ToConc = 10 ** (6) / (self.SurfaceRunoff * self.timestepsecs) + self.SedConc = self.OutSedLoad * self.ToConc + self.SSConc = (self.OutClayLoad + self.OutSiltLoad) * self.ToConc + self.BedConc = ( + self.OutSandLoad + + self.OutSaggLoad + + self.OutLaggLoad + + self.OutGravLoad + ) * self.ToConc + self.MaxSedConc = self.MaxSedLoad * self.ToConc + + """ Check / summary maps """ + # Checking processes contribution to sediment budget + self.PrecipitationCatch = areatotal(self.Precipitation, self.HYPEcatch) + # self.InLandSedCatch = areatotal(self.InLandSed, self.HYPEcatch) + # self.MaxSedLoadCatch = areatotal(self.MaxSedLoad, self.HYPEcatch) + # self.SedExCatch = areatotal(self.SedEx, self.HYPEcatch) + # self.BankSedLoadCatch = areatotal(self.BankSedLoad, self.HYPEcatch) + # self.BedSedLoadCatch = areatotal(self.BedSedLoad, self.HYPEcatch) + # self.DegStoreSedCatch = areatotal(self.DegStoreSed, self.HYPEcatch) + # self.DepSedLoadCatch = areatotal(self.DepSedLoad, self.HYPEcatch) + # self.OutSedLoadCatch = areatotal(self.OutSedLoad, self.HYPEcatch) + # self.SedLoadCatch = areatotal(self.SedLoad, self.HYPEcatch) + + # Area specific runoff + # annual runoff / upstream (area) + # L/km2/s + # values between 5 and 25 for the Rhine + self.SpecRun = ( + self.SurfaceRunoff * 1000 / accuflux(self.TopoLdd, self.cellareaKm) + ) # [L/km2/s] + + # The main function is used to run the program from the command line -def main(argv=None): + +def main(argv=None): """ *Optional but needed it you want to run the model from the command line* Perform command line execution of the model. This example uses the getopt module to parse the command line options. The user can set the caseName, the runDir, the timestep and the configfile. - """ + """ global multpars caseName = "default_sediment" runId = "run_default" - configfile="wflow_sediment.ini" + configfile = "wflow_sediment.ini" _lastTimeStep = 0 _firstTimeStep = 0 - timestepsecs=86400 - wflow_cloneMap = 'wflow_subcatch.map' - - # This allows us to use the model both on the command line and to call + timestepsecs = 86400 + wflow_cloneMap = "wflow_subcatch.map" + + # This allows us to use the model both on the command line and to call # the model using main function from another python script. - + if argv is None: argv = sys.argv[1:] if len(argv) == 0: usage() - return + return - opts, args = getopt.getopt(argv, 'C:S:T:c:s:R:') - + opts, args = getopt.getopt(argv, "C:S:T:c:s:R:") + for o, a in opts: - if o == '-C': caseName = a - if o == '-R': runId = a - if o == '-c': configfile = a - if o == '-s': timestepsecs = int(a) - if o == '-T': _lastTimeStep=int(a) - if o == '-S': _firstTimeStep=int(a) - - if (len(opts) <=1): + if o == "-C": + caseName = a + if o == "-R": + runId = a + if o == "-c": + configfile = a + if o == "-s": + timestepsecs = int(a) + if o == "-T": + _lastTimeStep = int(a) + if o == "-S": + _firstTimeStep = int(a) + + if len(opts) <= 1: usage() - - myModel = WflowModel(wflow_cloneMap, caseName,runId,configfile) - dynModelFw = wf_DynamicFramework(myModel, _lastTimeStep,firstTimestep=_firstTimeStep) - dynModelFw.createRunId(NoOverWrite=False,level=logging.DEBUG) + + myModel = WflowModel(wflow_cloneMap, caseName, runId, configfile) + dynModelFw = wf_DynamicFramework( + myModel, _lastTimeStep, firstTimestep=_firstTimeStep + ) + dynModelFw.createRunId(NoOverWrite=False, level=logging.DEBUG) dynModelFw._runInitial() dynModelFw._runResume() - #dynModelFw._runDynamic(0,0) + # dynModelFw._runDynamic(0,0) dynModelFw._runDynamic(_firstTimeStep, _lastTimeStep) dynModelFw._runSuspend() dynModelFw._wf_shutdown() - + if __name__ == "__main__": main() Index: wflow/wflow_sphy.py =================================================================== diff -u -r2624202607c719d443d22173d308230fd8196fd3 -r11471548f856592485a77ca37df4ad838eed2eaf --- wflow/wflow_sphy.py (.../wflow_sphy.py) (revision 2624202607c719d443d22173d308230fd8196fd3) +++ wflow/wflow_sphy.py (.../wflow_sphy.py) (revision 11471548f856592485a77ca37df4ad838eed2eaf) @@ -43,6 +43,7 @@ updateCols = [] #: columns used in updating """ Column used in updating """ + def usage(*args): """ Print usage information @@ -72,7 +73,7 @@ self.Dir = os.path.abspath(Dir) self.configfile = configfile self.SaveDir = os.path.join(self.Dir, self.runId) - + def stateVariables(self): """ returns a list of state variables that are essential to the model. @@ -100,13 +101,8 @@ baseflow waters in lakes and reservoirs [m3] """ - states = [ - "RootWater", - "SubWater", - "CapRise", - "RootDrain" - ] - + states = ["RootWater", "SubWater", "CapRise", "RootDrain"] + if self.DynVegFLAG == 1: states.append("Scanopy") if self.GroundFLAG == 1: @@ -117,7 +113,7 @@ states.extend(["SnowStore", "SnowWatStore"]) if self.GlacFLAG == 1: states.append("GlacFrac") - if (self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1): + if self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1: states.append("QRAold") if self.RainRA_FLAG == 1: states.append("RainRAold") @@ -127,7 +123,7 @@ states.append("GlacRAold") if self.BaseRA_FLAG == 1: states.append("BaseRAold") - if (self.LakeFLAG == 1 or self.ResFLAG == 1): + if self.LakeFLAG == 1 or self.ResFLAG == 1: states.append("StorRES") if self.RainRA_FLAG == 1: states.append("RainRAstor") @@ -137,10 +133,9 @@ states.append("GlacRAstor") if self.BaseRA_FLAG == 1: states.append("BaseRAstor") - - + return states - + # The following are made to better connect to deltashell/openmi def supplyCurrentTime(self): # - this may not be required """ @@ -152,7 +147,7 @@ return self.currentTimeStep() * int( configget(self.config, "model", "timestepsecs", "86400") ) - + def parameters(self): """ Define all model parameters here that the framework should handle for the model @@ -171,11 +166,14 @@ self.config, "inputmapstacks", "Prec", "/inmaps/prec" ) # timeseries for rainfall self.Tair_mapstack = self.Dir + configget( - self.config, "inputmapstacks", "Tair", "/inmaps/tavg") + self.config, "inputmapstacks", "Tair", "/inmaps/tavg" + ) self.Tmax_mapstack = self.Dir + configget( - self.config, "inputmapstacks", "Tmax", "/inmaps/tmax") + self.config, "inputmapstacks", "Tmax", "/inmaps/tmax" + ) self.Tmin_mapstack = self.Dir + configget( - self.config, "inputmapstacks", "Tmin", "/inmaps/tmin") + self.config, "inputmapstacks", "Tmin", "/inmaps/tmin" + ) # Meteo and other forcing modelparameters.append( @@ -198,7 +196,7 @@ lookupmaps=[], ) ) - + modelparameters.append( self.ParamType( name="Tmax", @@ -219,10 +217,9 @@ lookupmaps=[], ) ) - + return modelparameters - - + def suspend(self): """ Suspends the model to disk. All variables needed to restart the model @@ -239,8 +236,7 @@ if self.fewsrun: self.logger.info("Saving initial conditions for FEWS...") self.wf_suspend(os.path.join(self.Dir, "outstate")) - - + def initial(self): global statistics @@ -250,22 +246,22 @@ setglobaloption("unittrue") self.thestep = scalar(0) - + self.logger.info("running for " + str(self.nrTimeSteps()) + " timesteps") # Set and get defaults from ConfigFile here ################################### self.reinit = int(configget(self.config, "run", "reinit", "0")) self.fewsrun = int(configget(self.config, "run", "fewsrun", "0")) self.OverWriteInit = int(configget(self.config, "model", "OverWriteInit", "0")) - + # Print model info print( "The Spatial Processes in HYdrology (SPHY) model is " "developed and owned by FutureWater, Wageningen, The Netherlands" ) print("Version 2.1") print(" ") - + # Read the modules to be used self.GlacFLAG = int(configget(self.config, "model", "GlacFLAG", "0")) self.SnowFLAG = int(configget(self.config, "model", "SnowFLAG", "0")) @@ -275,18 +271,18 @@ self.DynVegFLAG = int(configget(self.config, "model", "DynVegFLAG", "0")) self.GroundFLAG = int(configget(self.config, "model", "GroundFLAG", "0")) # Optional modules - if (self.GroundFLAG == 0): + if self.GroundFLAG == 0: self.SeepStatFLAG = int(configget(self.config, "model", "SeepStatic", "0")) if self.DynVegFLAG == 0: self.KcStatFLAG = int(configget(self.config, "model", "KCstatic", "0")) self.ETREF_FLAG = int(configget(self.config, "model", "ETREF_FLAG", "0")) - - if (self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1): + + if self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1: self.SnowRA_FLAG = int(configget(self.config, "model", "SnowRA_FLAG", "0")) self.RainRA_FLAG = int(configget(self.config, "model", "RainRA_FLAG", "0")) self.GlacRA_FLAG = int(configget(self.config, "model", "GlacRA_FLAG", "0")) self.BaseRA_FLAG = int(configget(self.config, "model", "BaseRA_FLAG", "0")) - + # import the required modules import datetime, calendar from wflow.sphy import reporting as reporting @@ -309,7 +305,7 @@ self.rootzone = rootzone self.subzone = subzone del datetime, calendar, pi, reporting, timecalc, ET, rootzone, subzone - + # -import additional modules if required if self.GlacFLAG == 1: self.SnowFLAG = 1 @@ -334,12 +330,12 @@ self.lakes = lakes del lakes if self.ResFLAG == 1: - import wflow.sphy.reservoirs as reservoirs # import reservoir module + import wflow.sphy.reservoirs as reservoirs # import reservoir module self.reservoirs = reservoirs del reservoirs if self.LakeFLAG == 1 or self.ResFLAG == 1: - import wflow.sphy.advanced_routing as advanced_routing # overwrite the simple routing scheme + import wflow.sphy.advanced_routing as advanced_routing # overwrite the simple routing scheme self.routing = advanced_routing del advanced_routing @@ -354,8 +350,7 @@ self.groundwater = groundwater del groundwater - - + # -set the global options setglobaloption("radians") # -set the 2000 julian date number @@ -364,24 +359,18 @@ self.mm_rep_FLAG = int(configget(self.config, "model", "mm_rep_FLAG", "1")) # -convert flow from m3/s to mm self.ToMM = 1000 * 3600 * 24 / cellarea() - + # static maps to use (normally default) - wflow_dem = configget( - self.config, "model", "dem", "staticmaps/dem.map" - ) - wflow_slope = configget( - self.config, "model", "slope", "staticmaps/slope.map" - ) + wflow_dem = configget(self.config, "model", "dem", "staticmaps/dem.map") + wflow_slope = configget(self.config, "model", "slope", "staticmaps/slope.map") wflow_locs = configget( self.config, "model", "locations", "staticmaps/outlets.map" ) wflow_landuse = configget( self.config, "model", "landuse", "staticmaps/landuse.map" ) - wflow_ldd = configget( - self.config, "model", "flowdir", "staticmaps/ldd.map" - ) - + wflow_ldd = configget(self.config, "model", "flowdir", "staticmaps/ldd.map") + wflow_rootF = configget( self.config, "model", "RootFieldMap", "staticmaps/root_field.map" ) @@ -406,11 +395,9 @@ wflow_subK = configget( self.config, "model", "SubKsat", "staticmaps/sub_ksat.map" ) - + # 2: Input base maps ######################################################## - self.DEM = self.wf_readmap( - os.path.join(self.Dir, wflow_dem), 0.0, fail=True - ) + self.DEM = self.wf_readmap(os.path.join(self.Dir, wflow_dem), 0.0, fail=True) self.Slope = self.wf_readmap( os.path.join(self.Dir, wflow_slope), 0.0, fail=True ) @@ -423,7 +410,7 @@ self.FlowDir = self.wf_readmap( os.path.join(self.Dir, wflow_ldd), 0.0, fail=False ) - + self.RootFieldMap = self.wf_readmap( os.path.join(self.Dir, wflow_rootF), 0.35, fail=False ) @@ -448,26 +435,24 @@ self.SubKsat = self.wf_readmap( os.path.join(self.Dir, wflow_subK), 10, fail=False ) - - + # Set static initial values here ######################################### self.ZeroMap = 0.0 * scalar(self.DEM) # map with only zero's self.Latitude = ycoordinate(boolean(self.ZeroMap)) - self.Longitude = xcoordinate(boolean(self.ZeroMap)) + self.Longitude = xcoordinate(boolean(self.ZeroMap)) # Read parameters NEW Method self.logger.info("Linking parameters to landuse, catchment and soil...") self.wf_updateparameters() self.wf_multparameters() - - + # Set static initial variables self.RootDrainVel = self.RootKsat * self.Slope - - if (self.GroundFLAG == 0): + + if self.GroundFLAG == 0: self.SubDrainVel = self.SubKsat * self.Slope - + # -calculate soil properties self.RootField = self.RootFieldMap * self.RootDepthFlat self.RootSat = self.RootSatMap * self.RootDepthFlat @@ -477,17 +462,18 @@ self.SubField = self.SubFieldMap * self.SubDepthFlat self.RootTT = (self.RootSat - self.RootField) / self.RootKsat self.SubTT = (self.SubSat - self.SubField) / self.SubKsat - + # soil max and soil min for scaling of gwl if groundwater module is not used if self.GroundFLAG == 0: self.SoilMax = self.RootSat + self.SubSat self.SoilMin = self.RootDry + self.SubField - + if self.ETREF_FLAG == 0: from wflow.sphy import hargreaves + self.Hargreaves = hargreaves del hargreaves - + setglobaloption("matrixtable") # -read lake maps and parameters if lake module is used if self.LakeFLAG == 1: @@ -574,7 +560,6 @@ self.QFRAC = ifthenelse(self.ResID != 0, scalar(0), self.QFRAC) else: self.QFRAC = ifthenelse(self.ResID != 0, scalar(0), 1) - def default_summarymaps(self): ##-maybe not needed. check later """ @@ -590,13 +575,10 @@ if self.GlacFLAG == 1: lst = ["self.GlacFrac"] else: - lst =[] + lst = [] - return lst - + return lst - - def resume(self): """ read initial state maps (they are output of a previous call to suspend()) """ @@ -614,7 +596,7 @@ # -initial canopy storage self.Scanopy = self.ZeroMap # -initial ndvi if first map is not provided - #self.ndviOld = scalar((self.NDVImax + self.NDVImin) / 2) + # self.ndviOld = scalar((self.NDVImax + self.NDVImin) / 2) # -initial groundwater properties if self.GroundFLAG == 1: self.GwRecharge = self.ZeroMap + 2.0 @@ -630,56 +612,97 @@ # -initial glacier properties if self.GlacFLAG == 1: self.GlacFrac = self.ZeroMap - # -initial routed total runoff - if (self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1): + # -initial routed total runoff + if self.RoutFLAG == 1 or self.LakeFLAG == 1 or self.ResFLAG == 1: self.QRAold = self.ZeroMap self.RainRAold = self.ZeroMap self.SnowRAold = self.ZeroMap self.GlacRAold = self.ZeroMap self.BaseRAold = self.ZeroMap # -initial storage in lakes and reservoirs - if (self.LakeFLAG == 1 or self.ResFLAG == 1): + if self.LakeFLAG == 1 or self.ResFLAG == 1: # -Read initial storages from table/reservoir file - self.StorRES = self.ZeroMap - self.RainRAstor = self.ZeroMap + self.StorRES = self.ZeroMap + self.RainRAstor = self.ZeroMap self.SnowRAstor = self.ZeroMap - self.GlacRAstor = self.ZeroMap - self.BaseRAstor = self.ZeroMap - + self.GlacRAstor = self.ZeroMap + self.BaseRAstor = self.ZeroMap + # -Read initial storages from table/reservoir file if self.LakeFLAG == 1: LakeStor_Tab = os.path.join(self.Dir, "lake_id.tbl") if os.path.exists(LakeStor_Tab): - self.StorRES = self.StorRES + cover(lookupscalar(LakeStor_Tab, 1, self.LakeID), 0) * 10 ** 6 # convert to m3 - self.RainRAstor = self.RainRAstor + cover(lookupscalar(LakeStor_Tab, 2, self.LakeID), 0) * 10 ** 6 - self.SnowRAstor = self.SnowRAstor + cover(lookupscalar(LakeStor_Tab, 3, self.LakeID), 0) * 10 ** 6 - self.GlacRAstor = self.GlacRAstor + cover(lookupscalar(LakeStor_Tab, 4, self.LakeID), 0) * 10 ** 6 - self.BaseRAstor = self.BaseRAstor + cover(lookupscalar(LakeStor_Tab, 5, self.LakeID), 0) * 10 ** 6 + self.StorRES = ( + self.StorRES + + cover(lookupscalar(LakeStor_Tab, 1, self.LakeID), 0) + * 10 ** 6 + ) # convert to m3 + self.RainRAstor = ( + self.RainRAstor + + cover(lookupscalar(LakeStor_Tab, 2, self.LakeID), 0) + * 10 ** 6 + ) + self.SnowRAstor = ( + self.SnowRAstor + + cover(lookupscalar(LakeStor_Tab, 3, self.LakeID), 0) + * 10 ** 6 + ) + self.GlacRAstor = ( + self.GlacRAstor + + cover(lookupscalar(LakeStor_Tab, 4, self.LakeID), 0) + * 10 ** 6 + ) + self.BaseRAstor = ( + self.BaseRAstor + + cover(lookupscalar(LakeStor_Tab, 5, self.LakeID), 0) + * 10 ** 6 + ) else: - self.logger.debug("Initial default state lake_id.tbl not found returning 0.0") + self.logger.debug( + "Initial default state lake_id.tbl not found returning 0.0" + ) if self.ResFLAG == 1: ResStor_Tab = os.path.join(self.Dir, "reservoir_id.tbl") if os.path.exists(ResStor_Tab): - self.StorRES = self.StorRES + cover(lookupscalar(ResStor_Tab, 2, self.ResID), 0) * 10 ** 6 - self.RainRAstor = self.RainRAstor + cover(lookupscalar(ResStor_Tab, 3, self.ResID), 0) * 10 ** 6 - self.SnowRAstor = self.SnowRAstor + cover(lookupscalar(ResStor_Tab, 4, self.ResID), 0) * 10 ** 6 - self.GlacRAstor = self.GlacRAstor + cover(lookupscalar(ResStor_Tab, 5, self.ResID), 0) * 10 ** 6 - self.BaseRAstor = self.BaseRAstor + cover(lookupscalar(ResStor_Tab, 6, self.ResID), 0) * 10 ** 6 + self.StorRES = ( + self.StorRES + + cover(lookupscalar(ResStor_Tab, 2, self.ResID), 0) + * 10 ** 6 + ) + self.RainRAstor = ( + self.RainRAstor + + cover(lookupscalar(ResStor_Tab, 3, self.ResID), 0) + * 10 ** 6 + ) + self.SnowRAstor = ( + self.SnowRAstor + + cover(lookupscalar(ResStor_Tab, 4, self.ResID), 0) + * 10 ** 6 + ) + self.GlacRAstor = ( + self.GlacRAstor + + cover(lookupscalar(ResStor_Tab, 5, self.ResID), 0) + * 10 ** 6 + ) + self.BaseRAstor = ( + self.BaseRAstor + + cover(lookupscalar(ResStor_Tab, 6, self.ResID), 0) + * 10 ** 6 + ) else: - self.logger.debug("Initial default state res_id.tbl not found returning 0.0") - + self.logger.debug( + "Initial default state res_id.tbl not found returning 0.0" + ) else: self.wf_resume(os.path.join(self.Dir, "instate")) - + # -initial water storage in rootzone + subsoil self.SoilWater = self.RootWater + self.SubWater - + if self.SnowFLAG == 1: self.TotalSnowStore = self.SnowStore + self.SnowWatStore - - def dynamic(self): self.wf_updateparameters() # read forcing an dynamic parameters @@ -690,12 +713,12 @@ if self.SnowFLAG == 0: self.SnowStore = scalar(0) SnowFrac = ifthenelse(self.SnowStore > 0, scalar(1 - self.GlacFrac), 0) - RainFrac = ifthenelse(self.SnowStore == 0, scalar(1 - self.GlacFrac), 0) - + RainFrac = ifthenelse(self.SnowStore == 0, scalar(1 - self.GlacFrac), 0) + # -Read the precipitation time-series self.Precip = self.Prec self.PrecipF = self.Prec * (1 - self.GlacFrac) - + # -Temperature and reference evapotranspiration Temp = self.Tair if self.ETREF_FLAG == 0: @@ -704,10 +727,10 @@ self.ETref = self.Hargreaves.Hargreaves( pcr, self.Hargreaves.extrarad(self, pcr), Temp, TempMax, TempMin ) - + # -Interception and effective precipitation # -Update canopy storage - if self.DynVegFLAG == 1: + if self.DynVegFLAG == 1: # -fill missing ndvi values with ndvi base self.NDVI = ifthenelse(defined(self.NDVI) == 1, self.NDVI, self.NDVIbase) # -calculate the vegetation parameters @@ -741,29 +764,38 @@ # -canopy storage self.Scanopy = intercep[2] - # Snow and rain if self.SnowFLAG == 1: # -Snow and rain differentiation self.Snow = ifthenelse(Temp >= self.Tcrit, 0, self.Precip) self.SnowF = self.Snow * (1 - self.GlacFrac) self.Rain = ifthenelse(Temp < self.Tcrit, 0, self.Precip) - + # -Snow melt PotSnowMelt = self.snow.PotSnowMelt(pcr, Temp, self.DDFS) self.ActSnowMelt = self.snow.ActSnowMelt(pcr, self.SnowStore, PotSnowMelt) self.ActSnowMeltF = self.ActSnowMelt * SnowFrac # -Update snow store self.SnowStore = self.snow.SnowStoreUpdate( - pcr, self.SnowStore, self.Snow, self.ActSnowMelt, Temp, self.SnowWatStore + pcr, + self.SnowStore, + self.Snow, + self.ActSnowMelt, + Temp, + self.SnowWatStore, ) # -Caclulate the maximum amount of water that can be stored in snowwatstore MaxSnowWatStore = self.snow.MaxSnowWatStorage(self.SnowSC, self.SnowStore) OldSnowWatStore = self.SnowWatStore # -Calculate the actual amount of water stored in snowwatstore self.SnowWatStore = self.snow.SnowWatStorage( - pcr, Temp, MaxSnowWatStore, self.SnowWatStore, self.ActSnowMelt, self.Rain + pcr, + Temp, + MaxSnowWatStore, + self.SnowWatStore, + self.ActSnowMelt, + self.Rain, ) # -Changes in total water storage in snow (SnowStore and SnowWatStore) OldTotalSnowStore = self.TotalSnowStore @@ -786,7 +818,7 @@ OldTotalSnowStore = 0 self.TotalSnowStore = 0 self.RainF = self.Rain * (1 - self.GlacFrac) - + # -Glacier calculations if self.GlacFLAG == 1: # -Glacier melt from clean ice glaciers @@ -804,8 +836,8 @@ else: self.GlacR = 0 self.GlacMelt = 0 - self.GlacPerc = 0 - + self.GlacPerc = 0 + # -Potential evapotranspiration (THIS SHOULD STILL BE IMPROVED WITH DYNAMIC VEGETATION MODULE) self.ETpot = self.ET.ETpot(self.ETref, self.Kc) self.ETpotF = self.ETpot * RainFrac @@ -848,8 +880,8 @@ ) self.rootpercF = self.rootperc * (1 - self.GlacFrac) # -Update rootwater content - self.RootWater = self.RootWater - self.rootperc - + self.RootWater = self.RootWater - self.rootperc + # -Sub soil calculations self.SubWater = self.SubWater + self.rootperc if self.GroundFLAG == 0: @@ -865,7 +897,7 @@ self.RootField, ) self.CapRiseF = self.CapRise * (1 - self.GlacFrac) - + # -Update sub soil water content self.SubWater = self.SubWater - self.CapRise if ( @@ -889,7 +921,7 @@ ) # -Update sub soil water content self.SubWater = self.SubWater - self.SubDrain - + # -Changes in soil water storage OldSoilWater = self.SoilWater self.SoilWater = (self.RootWater + self.SubWater) * (1 - self.GlacFrac) @@ -900,7 +932,7 @@ self.RootD = self.RootDrain * (1 - self.GlacFrac) # -Rain runoff self.RainR = self.RootR + self.RootD - + # -Groundwater calculations if self.GroundFLAG == 1: GwOld = self.Gw @@ -920,7 +952,10 @@ self.H_gw = self.groundwater.HLevel( pcr, self.H_gw, self.alphaGw, self.GwRecharge, self.YieldGw ) - self.GWL = ((self.SubDepthFlat + self.RootDepthFlat + self.GwDepth) / 1000 - self.H_gw) * (-1) + self.GWL = ( + (self.SubDepthFlat + self.RootDepthFlat + self.GwDepth) / 1000 + - self.H_gw + ) * (-1) else: # -Use drainage from subsoil as baseflow self.BaseR = self.SubDrain @@ -931,8 +966,8 @@ ) # scale between 0 (dry) and 1 (wet) self.GWL = self.GWL_base - (SoilRel - 0.5) * self.GWL_base - self.TotRF = self.BaseR + self.RainR + self.SnowR + self.GlacR - + self.TotRF = self.BaseR + self.RainR + self.SnowR + self.GlacR + # -Water balance if self.GroundFLAG == 1: self.waterbalance = ( @@ -964,7 +999,9 @@ # -Update storage in lakes/reservoirs (m3) with specific runoff self.StorRES = self.StorRES + ifthenelse( self.QFRAC == 0, - 0.001 * cellarea() * (self.BaseR + self.RainR + self.GlacR + self.SnowR), + 0.001 + * cellarea() + * (self.BaseR + self.RainR + self.GlacR + self.SnowR), 0, ) OldStorage = self.StorRES @@ -984,13 +1021,15 @@ self.StorRES = tempvar[1] self.Qout = self.lakes.QLake(self, pcr, LakeLevel) else: - self.Qout = self.reservoirs.QRes(self, pcr) - + self.Qout = self.reservoirs.QRes(self, pcr) + # -Calculate volume available for routing (=outflow lakes/reservoir + cell specific runoff) RunoffVolume = upstream(self.FlowDir, self.Qout) + ifthenelse( self.QFRAC == 0, 0, - 0.001 * cellarea() * (self.BaseR + self.RainR + self.GlacR + self.SnowR), + 0.001 + * cellarea() + * (self.BaseR + self.RainR + self.GlacR + self.SnowR), ) # -Routing of total flow tempvar = self.routing.ROUT( @@ -999,8 +1038,8 @@ self.StorRES = tempvar[0] self.Q = tempvar[1] self.Qin = tempvar[2] - self.QRAold = self.Q - + self.QRAold = self.Q + # -Routing of individual contributers # -Snow routing if self.SnowRA_FLAG == 1 and self.SnowFLAG == 1: @@ -1013,13 +1052,18 @@ self.QFRAC == 0, 0, 0.001 * cellarea() * self.SnowR ) tempvar = self.routing.ROUT( - self, pcr, cRunoffVolume, self.SnowRAold, self.cQout, self.SnowRAstor + self, + pcr, + cRunoffVolume, + self.SnowRAold, + self.cQout, + self.SnowRAstor, ) self.SnowRAstor = tempvar[0] self.SnowRA = tempvar[1] self.cQin = tempvar[2] - self.SnowRAold = self.SnowRA - + self.SnowRAold = self.SnowRA + # -Rain routing if self.RainRA_FLAG == 1: self.RainRAstor = self.RainRAstor + ifthenelse( @@ -1031,14 +1075,19 @@ self.QFRAC == 0, 0, 0.001 * cellarea() * self.RainR ) tempvar = self.routing.ROUT( - self, pcr, cRunoffVolume, self.RainRAold, self.cQout, self.RainRAstor + self, + pcr, + cRunoffVolume, + self.RainRAold, + self.cQout, + self.RainRAstor, ) self.RainRAstor = tempvar[0] self.RainRA = tempvar[1] self.cQin = tempvar[2] self.RainRAold = self.RainRA - - # -Glacier routing + + # -Glacier routing if self.GlacRA_FLAG == 1 and self.GlacFLAG == 1: self.GlacRAstor = self.GlacRAstor + ifthenelse( self.QFRAC == 0, self.GlacR * 0.001 * cellarea(), 0 @@ -1049,13 +1098,18 @@ self.QFRAC == 0, 0, 0.001 * cellarea() * self.GlacR ) tempvar = self.routing.ROUT( - self, pcr, cRunoffVolume, self.GlacRAold, self.cQout, self.GlacRAstor + self, + pcr, + cRunoffVolume, + self.GlacRAold, + self.cQout, + self.GlacRAstor, ) self.GlacRAstor = tempvar[0] self.GlacRA = tempvar[1] self.cQin = tempvar[2] self.GlacRAold = self.GlacRA - + # -Baseflow routing if self.BaseRA_FLAG == 1: self.BaseRAstor = self.BaseRAstor + ifthenelse( @@ -1067,13 +1121,18 @@ self.QFRAC == 0, 0, 0.001 * cellarea() * self.BaseR ) tempvar = self.routing.ROUT( - self, pcr, cRunoffVolume, self.BaseRAold, self.cQout, self.BaseRAstor + self, + pcr, + cRunoffVolume, + self.BaseRAold, + self.cQout, + self.BaseRAstor, ) self.BaseRAstor = tempvar[0] self.BaseRA = tempvar[1] self.cQin = tempvar[2] self.BaseRAold = self.BaseRA - + # -Normal routing module elif self.RoutFLAG == 1: # -Rout total runoff @@ -1091,7 +1150,7 @@ pcr, self.SnowR, self.SnowRAold, self.FlowDir, self.kx ) self.SnowRAold = self.SnowRA - + # -Rain routing if self.RainRA_FLAG == 1: self.RainRA = self.routing.ROUT( @@ -1109,10 +1168,9 @@ self.BaseRA = self.routing.ROUT( pcr, self.BaseR, self.BaseRAold, self.FlowDir, self.kx ) - self.BaseRAold = self.BaseRA - - - + self.BaseRAold = self.BaseRA + + # The main function is used to run the program from the command line @@ -1258,4 +1316,3 @@ if __name__ == "__main__": main() - \ No newline at end of file