Index: wflow-py/Sandbox/wflow_usle.py =================================================================== diff -u -r8f72da8614533d8676e5cb914d1e52fc08569ed8 -r8310be73f40010119b0079f9d7dbe38bfb53016f --- wflow-py/Sandbox/wflow_usle.py (.../wflow_usle.py) (revision 8f72da8614533d8676e5cb914d1e52fc08569ed8) +++ wflow-py/Sandbox/wflow_usle.py (.../wflow_usle.py) (revision 8310be73f40010119b0079f9d7dbe38bfb53016f) @@ -16,37 +16,30 @@ -c name of the config file (in the case directory) -Steps: -1) Determie soil factor K (alsready in script) +This script assumes monthly timesteps. -2) make monthly P values and determine P factor -This script runs on monthly timesteps!!! And apply Fournier Index, +usle_r, Rainfall erosivity is determined according to: +Beskow, S., Mello, C.R., Norton, L.D., Curi, N., Viola, M.R., Avanzi, J.C., 2009. +Soil erosion prediction in the Grande River Basin, Brazil using distributed modeling. +CATENA 79, 49–59. doi:10.1016/j.catena.2009.05.010 +Determine m using: +------------------ -3) 2.3.3. Cover-management factor (C) +Slope gradient Value of m +1% 0.2 +1–3% 0.3 +3–4.5% 0.4 +4.5% or more 0.5 -4) 2.3.4. Topographic factor (LS) +Determien L using the cell-length -ad 2) +self.usle_l = (self.Lambda/22.13)**self.m # The L factor +# Now determime S using: S=(0.43+0.30s+0.043s^2)/6.613 # s = slope +# Wischmeier and Smith (1978), +self.usle_s = (0.43+ 0.30*(self.Slope*100) + 0.043 * (self.Slope * 100)**2)/6.613 -known as Fournier Index, was applied. This Index has been used -widely in several studies of soil loss and erosivity mapping, as in Irvem -et al. (2007), Mello et al. (2007) and Pandey et al. (2007). -EIi = -125:92 × r2 -iP -  -0:603 -+ 111:173 × r2 -iP -  -0:691 -+ 68:73 × r2 -iP -  -0:841 -3 """ @@ -119,9 +112,10 @@ # These should be linked to soil type modelparameters.append(self.ParamType(name="percent_clay",stack="intbl/percent_clay.tbl",type="statictbl",default=0.1, verbose=False,lookupmaps=[])) modelparameters.append(self.ParamType(name="percent_silt",stack="intbl/percent_silt.tbl",type="statictbl",default=0.1, verbose=False,lookupmaps=[])) + modelparameters.append(self.ParamType(name="percent_oc",stack="intbl/percent_oc.tbl",type="statictbl",default=0.01, verbose=False,lookupmaps=[])) # Sediment delivery ratio - modelparameters.append(self.ParamType(name="dratio",stack="intbl/dratio.tbl",type="statictbl",default=1.0, verbose=False,lookupmaps=[])) +# modelparameters.append(self.ParamType(name="dratio",stack="intbl/dratio.tbl",type="statictbl",default=1.0, verbose=False,lookupmaps=[])) #modelparameters.append(self.ParamType(name="usle_k", stack="intbl/usle_k.tbl", type="statictbl", default=1.0, verbose=False,lookupmaps=[])) modelparameters.append(self.ParamType(name="usle_c", stack="intbl/usle_c.tbl", type="statictbl", default=1.0, verbose=False,lookupmaps=[])) modelparameters.append(self.ParamType(name="usle_p", stack="intbl/usle_p.tbl", type="statictbl", default=1.0, verbose=False,lookupmaps=[])) @@ -184,8 +178,8 @@ """ 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 + #: 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.Dir + "/outstate/") @@ -228,6 +222,8 @@ # Calulate slope taking into account that x,y may be in lat,lon self.Slope = slope(self.Altitude) self.Slope = max(0.00001, self.Slope * celllength() / self.reallength) + # limit slope + self.DUSlope = min(1.0,max(0.005,self.Slope)) """ First determine m exponent based on Slope (https://www.researchgate.net/publication/226655635_Estimation_of_Soil_Erosion_for_a_Himalayan_Watershed_Using_GIS_Technique) @@ -242,26 +238,45 @@ # Now determime S using: S=(0.43+0.30s+0.043s^2)/6.613 # s = slope self.usle_s = (0.43+ 0.30*(self.Slope*100) + 0.043 * (self.Slope * 100)**2)/6.613 - self.percent_sand= 100-self.percent_clay-self.percent_silt - self.usle_k = ifthenelse(pcrand(self.percent_clay>=40.,pcrand(self.percent_sand>=20.,self.percent_sand<=45.)),2.0, - ifthenelse(pcrand(self.percent_clay>=27.,pcrand(self.percent_sand>=20.,self.percent_sand<=45.)),1.7, - ifthenelse(pcrand(self.percent_silt<=40.,self.percent_sand<=20.),2.0, - ifthenelse(pcrand(self.percent_silt>40.,self.percent_clay>=40.),1.6, - ifthenelse(pcrand(self.percent_clay>=35.,self.percent_sand>=45.),1.9, - ifthenelse(pcrand(self.percent_clay>=27.,self.percent_sand<20.),1.6, - ifthenelse(pcrand(self.percent_clay<=10.,self.percent_silt>=80.),1.2, - ifthenelse(self.percent_silt>=50,1.5, - ifthenelse(pcrand(self.percent_clay>=7.,pcrand(self.percent_sand<=52.,self.percent_silt>=28.)),2.0, - ifthenelse(self.percent_clay>=20.,2.1, - ifthenelse(self.percent_clay>=self.percent_sand-70.,2.6, - ifthenelse(self.percent_clay>=(2.*self.percent_sand)-170.,3,scalar(1.9))))))))))))) + self.percent_sand = 100-self.percent_clay-self.percent_silt + """ + Calculation of USLE K factor based on: + Maeda et al. (2010) - 'Potential impacts of agricultural expansion and climate change on soil erosion in the Eastern Arc Mountains of Kenya', doi:10.1016/j.geomorph.2010.07.019 + """ + self.SN = 1 - self.percent_sand / 100 + self.usle_k = (0.2 + 0.3*exp(-0.0256*self.percent_sand * (1 - self.percent_silt/100))) \ + * (self.percent_silt/(self.percent_clay+self.percent_silt))**0.3 \ + * (1 - (0.25*self.percent_oc)/(self.percent_oc + exp(3.72 - 2.95*self.percent_oc))) \ + * (1 - (0.7*self.SN) / (self.SN + exp(-5.51 + 22.9*self.SN))) + """ + Calculation of sediment delivery ratio based on channel: + https://peerj.com/preprints/2227.pdf + http://data.naturalcapitalproject.org/nightly-build/invest-users-guide/html/sdr.html + """ + self.unitareaupstr = catchmenttotal(1, self.TopoLdd) + self.Dup = catchmenttotal(self.usle_c,self.TopoLdd)/self.unitareaupstr * catchmenttotal(self.DUSlope,self.TopoLdd)/self.unitareaupstr *\ + sqrt(catchmenttotal(self.reallength,self.TopoLdd)) + + self.drainlength = detdrainlength(self.TopoLdd, self.xl, self.yl) + self.Ddn = self.drainlength/(self.usle_c * self.DUSlope) + self.IC = log10(self.Dup/self.Ddn) + self.SDRMax = 0.8 # (Vigiak et al., 2012) + self.IC0 = 0.5 + self.k = 2.0 + expfact = exp((self.IC0 - self.IC)/self.k) + self.SDR = self.SDRMax/(1 + expfact) + + self.SDR_area = 0.472 * catchmenttotal(self.reallength/1000.0,self.TopoLdd)**(-0.125) + + + self.logger.info("Starting Dynamic run...") self.thestep = 0 def resume(self): - """ + """ *Required* This function is required. Read initial state maps (they are output of a @@ -270,7 +285,7 @@ """ self.logger.info("Reading initial conditions...") - #: It is advised to use the wf_resume() function + #: It is advised to use the wf_resume() function #: here which pick up the variable save by a call to wf_suspend() try: self.wf_resume(self.Dir + "/instate/") @@ -304,22 +319,26 @@ rintnsty = self.Precipitation/self.timestepsecs - # Weird assumption for now, shoudl be a lookuptabel of LAI and landuse type... + # Weird assumption for now, should be a lookuptabel of LAI and landuse type... self.canopy_cover = min(1.0,self.LAI) # Determine erosivity from monthly rainfall and average yearly sum + # taken from Beskow, S., Mello, C.R., Norton, L.D., Curi, N., Viola, M.R., Avanzi, J.C., 2009. + # Soil erosion prediction in the Grande River Basin, Brazil using distributed modeling. + # CATENA 79, 49–59. doi:10.1016/j.catena.2009.05.010 + self.usle_r = (125.92 * (self.Precipitation/self.Pmean)**0.603 + \ 111.173 * (self.Precipitation/self.Pmean) ** 0.691 + \ 68.73 * (self.Precipitation/self.Pmean)** 0.841) / 3.0 self.SoilLoss = self.usle_l * self.usle_s * self.usle_k * self.usle_r *self.usle_c * self.usle_p + self.SoilLossAvgUpstr = catchmenttotal(self.SoilLoss,self.TopoLdd)/catchmenttotal(1,self.TopoLdd) + self.SoilLossTotalUpstr = catchmenttotal(self.SoilLoss,self.TopoLdd)/catchmenttotal((self.reallength * self.reallength)/10000.0,self.TopoLdd) + self.SedimentYieldUpstrAvg = self.SDR * self.SoilLossAvgUpstr # Average upstream of each pixel + self.SedimentYieldUpstrTot = self.SDR * self.SoilLossTotalUpstr # Total upstream of each pixel + self.SedimentYield = self.SDR * self.SoilLoss - #self.SedRunoff = accuflux(self.TopoLdd,self.hhsedy) - # limit downstream flow by surface runoff erosion rate - #self.SedStore = self.SedStore + self.hhsedy - #self.SedRunoff = accucapacityflux(self.TopoLdd, self.SedStore,self.sedov *20.0) - #self.SedStore = accucapacitystate(self.TopoLdd, self.SedStore,self.sedov * 20.0 )