subroutine sd_channel_sediment3 use sd_channel_module use hydrograph_module use time_module implicit none integer :: iob ! |object number real :: trap_eff !frac |trap efficiency in the flood plain real :: cohesion ! |soil bank cohesion real :: b_exp ! |exponent for bank erosion equation real :: vel_fall !m/s |fall velocity of sediment particles in channel real :: dep_fall !m |fall depth of sediment particles in channel real :: del_rto !frac |fraction of sediment deposited in channel real :: conc_chng ! |change in concentration (and mass) in channel sol and org N and P real :: ebtm_m !m |erosion of bottom of channel real :: ebank_m !m |meander cut on one side real :: ebtm_t !tons |bottom erosion real :: ebank_t !tons |bank erosion real :: shear_btm_cr ! | real :: shear_btm ! | real :: bf_flow !m3/s |bankfull flow rate * adjustment factor real :: pk_rto !ratio |peak to mean flow rate ratio real :: bd_fac ! |bulk density factor for critical velocity calculation real :: cohes_fac ! |cohesion factor for critical velocity calculation !real :: qman !m^3/s or m/s |flow rate or flow velocity real :: vel, veg, vel_cr, rad_curv, vel_bend, vel_rch, arc_len, prot_len, h_rad ich = isdch iob = sp_ob1%chandeg + jrch - 1 ebtm_m = 0. ebank_m = 0. ebtm_t = 0. ebank_t = 0. fp_dep = hz ch_dep = hz bank_ero = hz bed_ero = hz ch_trans = hz !! calculate channel sed and nutrient processes if inflow > 0 if (ht1%flo > 1.e-6) then !! Another eq from Peter - Qmax=Qmean*(1+2.66*Drainage Area^-.3) pk_rto = 0.2 + 0.5 / 250. * ob(icmd)%area_ha pk_rto = Max (1., pk_rto) pk_rto = 1. + 2.66 * (ob(icmd)%area_ha / 100.) ** -.3 !pk_rto = 1. + 1.33 * (ob(icmd)%area_ha / 100.) ** -.3 !pk_rto = 1. + 2. * (ob(icmd)%area_ha / 100.) ** -.3 !pk_rto = Min (2., pk_rto) !pk_rto = 1.5 peakrate = pk_rto * ht1%flo / 86400. !m3/s !! interpolate rating curve using peak rate call rcurv_interp_flo (ich, peakrate) !! use peakrate as flow rate h_rad = rcurv%xsec_area / rcurv%wet_perim sd_ch(ich)%chn = 0.5 + 0.2 * (ob(icmd)%area_ha / 100.) ** -.3 sd_ch(ich)%chn = Min (0.15, sd_ch(ich)%chn) sd_ch(ich)%chn = 0.39 * sd_ch(ich)%chs ** 0.38 * h_rad ** -0.16 sd_ch(ich)%chn = Min (0.15, sd_ch(ich)%chn) sd_ch(ich)%chn = Max (0.03, sd_ch(ich)%chn) vel = h_rad ** .6666 * Sqrt(sd_ch(ich)%chs) / (sd_ch(ich)%chn + .001) !vel = peakrate / rcurv%xsec_area !! compute flood plain deposition sd_ch(ich)%bankfull_flo = 1.75 bf_flow = sd_ch(ich)%bankfull_flo * ch_rcurv(ich)%elev(2)%flo_rate if (peakrate > bf_flow) then trap_eff = 0.24 * log(sd_ch(ich)%fp_inun_days) + 0.1 fp_dep%sed = trap_eff * ht1%sed !! deposit Particulate P and N in the floodplain fp_dep%orgn = trap_eff * sd_ch(ich)%n_dep_enr * ht1%orgn fp_dep%sedp = trap_eff * sd_ch(ich)%p_dep_enr * ht1%sedp !! trap nitrate and sol P in flood plain - when not simulating flood plain interactions? fp_dep%no3 = trap_eff * ht1%no3 fp_dep%solp = trap_eff * ht1%solp end if ch_morph(ich)%fp_mm = ch_morph(ich)%fp_mm + fp_dep%sed ht2 = ht1 - fp_dep !! calculate channel deposition based on fall velocity - SWRRB book !! assume particle size = 0.03 mm -- median silt size vel_fall = 411. * sd_ch(ich)%part_size ** 2 ! m/h dep_fall = vel_fall * rcurv%ttime !! assume bankfull flow depth if (dep_fall < sd_ch(ich)%chd) then del_rto = 1. - .5 * dep_fall / sd_ch(ich)%chd else del_rto = .5 * sd_ch(ich)%chd / dep_fall end if ch_dep%sed = (1. - del_rto) * ht2%sed ch_dep%orgn = sd_ch(ich)%n_dep_enr * (1. - del_rto) * ht2%orgn ch_dep%sedp = sd_ch(ich)%p_dep_enr * (1. - del_rto) * ht2%sedp ht2 = ht2 - ch_dep !! calculate in channel nutrient losses - input delivery ratio !! calculate in channel nutrient transformations conc_chng = 1. - exp(-sd_ch(ich)%n_sol_part * rcurv%ttime) ch_trans%no3 = conc_chng * ht2%no3 ch_trans%orgn = -ch_trans%no3 conc_chng = 1. - exp(-sd_ch(ich)%p_sol_part * rcurv%ttime) ch_trans%solp = conc_chng * ht2%solp ch_trans%sedp = -ch_trans%solp !! calc bank erosion cohesion = (15. / (15.3 - 0.438 * sd_ch(ich)%ch_clay + 0.0044 * sd_ch(ich)%ch_clay ** 2.)) * 1000. cohesion = amax1 (1000., cohesion) !min 1000 for sandy soils veg = 0. !Pa 200-10000. !function of cover factor bd_fac = Max (0.001, 0.03924 * sd_ch(ich)%ch_bd * 1000. - 1000.) cohes_fac = 0.021 * cohesion + veg vel_cr = log10 (2200. * sd_ch(ich)%chd) * (0.0004 * (bd_fac + 0.6 * cohes_fac)) ** 0.5 !cohesion = (-87.1 + (42.82 * sd_ch(ich)%ch_clay) - (0.267 * sd_ch(ich)%ch_clay ** 2.) & ! + (0.029 * sd_ch(ich)%ch_clay ** 3.)) * 2. !cohesion = (15. / (15.3 - 0.438 * sd_ch(ich)%ch_clay + 0.0044 * sd_ch(ich)%ch_clay ** 2.)) * 1000. !cohesion = amax1 (1000., cohesion) !min 1000 for sandy soils !veg = exp (-5. * sd_ch(ich)%chd) * cha1(icha)%dat%cov !function of cover factor !vel_cr = log10 (8.8 * sd_ch(ich)%chd / 0.004) * (0.0004 * ((sd_ch(ich)%ch_bd * & ! 1000. - 1000.) * 9.81 * 0.004 + 0.021 * cohesion + veg)) ** 0.5 !! calculate radius of curvature rad_curv = ((12. * sd_ch(ich)%chw) * sd_ch(ich)%sinu ** 1.5) / (13. * (sd_ch(ich)%sinu -1.) ** 0.5) vel_bend = vel * (1. / rad_curv + 1.) vel_rch = 0.33 * vel_bend + 0.66 * vel b_exp = 12.3 / sqrt (sd_ch(ich)%ch_clay + 1.) b_exp = min (3.5, b_exp) if (vel_rch > vel_cr) then ! .and. rcurv%dep / sd_ch(ich)%chd > 0.1) then ebank_m = 0.0024 * (vel_rch / vel_cr) ** b_exp !bank erosion m/yr else ebank_m = 0. end if ch_morph(ich)%w_yr = ch_morph(ich)%w_yr + ebank_m !! calc mass of sediment eroded -> t = bankcut (m) * depth (m) * lengthcut (m) * bd (t/m3) !! arc length = 0.33 * meander wavelength * sinuosity -> protected length arc_len = 0.33 * (12. * sd_ch(ich)%chw) * sd_ch(ich)%sinu prot_len = arc_len * sd_ch(ich)%arc_len_fr prot_len = 0.2 * sd_ch(ich)%chl * 1000. ebank_t = ebank_m * sd_ch(ich)%chd * prot_len * sd_ch(ich)%ch_bd ebank_t = 0.8 * ebank_t !assume 80% wash load and 20% bed deposition ebank_t = max (0., ebank_t) bank_ero%sed = ebank_t !! calculate associated nutrients bank_ero%orgn = bank_ero%sed * sd_ch(ich)%n_conc bank_ero%sedp = (1. - sd_ch(ich)%p_bio) * bank_ero%sed * sd_ch(ich)%p_conc bank_ero%no3 = 0. bank_ero%solp = sd_ch(ich)%p_bio * bank_ero%sed * sd_ch(ich)%p_conc bank_ero%no2 = 0. ht2 = ht2 + bank_ero !! calculate bed erosion !! no downcutting below equilibrium slope if (sd_ch(ich)%chs > 0.000001) then !sd_ch(ich)%chseq) then !! calc critical shear and shear on bottom of channel shear_btm_cr = sd_ch(ich)%d50 shear_btm = 9800. * rcurv%dep * sd_ch(ich)%chs !! Pa = N/m^2 * m * m/m !! critical shear function of d50 vel_cr = 0.293 * (sd_ch(ich)%d50) ** 0.5 vel_cr = 1. if (vel > vel_cr) then ebtm_m = 0.0001 * (vel_rch / vel_cr) ** 1.5 !bed erosion m/yr end if !! calc mass of sediment eroded -> t = m * width (m) * length (km) * 1000 m/km * bd (t/m3) ebtm_t = 1000. * ebtm_m * sd_ch(ich)%chw * sd_ch(ich)%chl * sd_ch(ich)%ch_bd end if ch_morph(ich)%d_yr = ch_morph(ich)%d_yr + ebtm_m bed_ero%sed = sd_ch(ich)%wash_bed_fr * ebtm_t !! calculate associated nutrients bed_ero%orgn = bed_ero%sed * sd_ch(ich)%n_conc bed_ero%sedp = (1. - sd_ch(ich)%p_bio) * bed_ero%sed * sd_ch(ich)%p_conc bed_ero%no3 = 0. bed_ero%solp = sd_ch(ich)%p_bio * bed_ero%sed * sd_ch(ich)%p_conc bed_ero%no2 = 0. ht2 = ht2 + bed_ero end if ! inflow>0 !! channel sediment budget for output ch_sed_bud(ich)%in_sed = ht1%sed ch_sed_bud(ich)%out_sed = ht2%sed ch_sed_bud(ich)%fp_dep = fp_dep%sed ch_sed_bud(ich)%ch_dep = ch_dep%sed ch_sed_bud(ich)%bank_ero = bank_ero%sed ch_sed_bud(ich)%bed_ero = bed_ero%sed !! channel nutrient budget for output ch_sed_bud(ich)%in_no3 = ht1%no3 ch_sed_bud(ich)%in_orgn = ht1%orgn ch_sed_bud(ich)%out_no3 = ht2%no3 ch_sed_bud(ich)%out_orgn = ht2%orgn ch_sed_bud(ich)%fp_no3 = fp_dep%no3 ch_sed_bud(ich)%bank_no3 = bank_ero%no3 ch_sed_bud(ich)%bed_no3 = bed_ero%no3 ch_sed_bud(ich)%fp_orgn = fp_dep%orgn ch_sed_bud(ich)%ch_orgn = ch_dep%orgn ch_sed_bud(ich)%bank_orgn = bank_ero%orgn ch_sed_bud(ich)%bed_orgn = bed_ero%orgn ch_sed_bud(ich)%in_solp = ht1%solp ch_sed_bud(ich)%in_orgp = ht1%sedp ch_sed_bud(ich)%out_solp = ht2%solp ch_sed_bud(ich)%out_orgp = ht2%sedp ch_sed_bud(ich)%fp_solp = fp_dep%solp ch_sed_bud(ich)%bank_solp = bank_ero%solp ch_sed_bud(ich)%bed_solp = bed_ero%solp ch_sed_bud(ich)%fp_orgp = fp_dep%sedp ch_sed_bud(ich)%ch_orgp = ch_dep%sedp ch_sed_bud(ich)%bank_orgp = bank_ero%sedp ch_sed_bud(ich)%bed_orgp = bed_ero%sedp ch_sed_bud(ich)%no3_orgn = ch_trans%no3 ch_sed_bud(ich)%solp_orgp = ch_trans%solp return end subroutine sd_channel_sediment3