# RegimeFunctions.r contains the functions that are used by the # UBC regime model # written by B. Eaton, 30 April 2015 #################################################################### RegimeModel = function(Q, mu, S, D50, D84, D90, fs){ # UBCRM.r is a function that runs the UBC Regime Model # Inputs are: # Q formative discharge (m3/s) # mu bank strength, relative to bed (afer Millar, 2005) # S energy gradient for the stream reach (m/m) # D50 median size of the bed surface (mm) influencing transport rate # D84 84th percentile of the bed surface (mmm) influencing resistance # D90 90th percentile of the bed surface (mm) influencing stability # fs percent of the bed surface that is covered by sand (%) # Outputs are: # W wetted width for the channel (m) # d mean hydraulic depth of the channel (m) # U mean flow velocity for the channel (m/s) # Q total discharge (m3/s) # Tbed shear stress acting on the channel bed (Pa) # Tbank shear stress acting on the channel banks (Pa) # Qb bed material sediment transport capacity (m3/s) # C maximum competent grain size (mm) # b trapezoid sideslope angle (degrees) # p bed width (m) # guess the convergence tolerance and initial bed width Tol = 0.00001 p = 4 * Q ^ 0.5 test_plus = FindStable(Q, mu, p*1.001, S, D50, D84, D90, fs) test_minus = FindStable(Q, mu, p*0.999, S, D50, D84, D90, fs) gradient_1 = test_plus$Qb - test_minus$Qb p1 = p #now move in the direction of the gradient if(gradient_1 > 0){p = p + 0.25*p} else {p = p - 0.25*p} test_plus = FindStable(Q, mu, p*1.001, S, D50, D84, D90, fs) test_minus = FindStable(Q, mu, p*0.999, S, D50, D84, D90, fs) gradient_2 = test_plus$Qb - test_minus$Qb p2 = p while(gradient_1/gradient_2 > 0){ gradient_1 = gradient_2 p1 = p if(gradient_1 > 0){p = p + 0.25*p} else {p = p - 0.25*p} test_plus = FindStable(Q, mu, p*1.001, S, D50, D84, D90, fs) test_minus = FindStable(Q, mu, p*0.999, S, D50, D84, D90, fs) gradient_2 = test_plus$Qb - test_minus$Qb p2 = p } p_upper = max(c(p1,p2)) p_lower = min(c(p1,p2)) p = 0.5*(p_upper + p_lower) converg = (p_upper-p_lower)/p while(converg > Tol){ test_plus = FindStable(Q, mu, p*1.001, S, D50, D84, D90, fs) test_minus = FindStable(Q, mu, p*0.999, S, D50, D84, D90, fs) gradient = test_plus$Qb - test_minus$Qb if(gradient > 0){p_lower = p} else {p_upper = p} p = 0.5*(p_upper + p_lower) converg = (p_upper-p_lower)/p } solution = FindStable(Q, mu, p, S, D50, D84, D90, fs) return(solution) } #################################################################### ChannelState = function(p, y, b, S, D50, D84, fs){ # calculates the characteristics of a trapezoidal channel with # specified dimensions and bed material. Inputs are: # p bed width (m) # y maximum flow depth for trapezoid (m) # b sideslope angle of the banks (degrees) # S energy gradient for the stream reach (m/m) # D50 median size of the bed surface (mm) # D84 84th percentile of the bed surface (mmm) # fs percent of the bed surface that is covered by sand (%) # Outputs are: # W wetted width for the channel (m) # d mean hydraulic depth of the channel (m) # U mean flow velocity for the channel (m/s) # Q total discharge (m3/s) # Tbed shear stress acting on the channel bed (Pa) # Tbank shear stress acting on the channel banks (Pa) # Qb bed material sediment transport capacity (m3/s) # C maximum competent grain size (mm) # specify the constants and make unit conversions g <<- 9.81 rho <<- 1000 # water density Gs <<- 1.65 # submerged specific gravity a1 = 6.5 a2 = 2.5 b = b * pi / 180 D50 = D50 / 1000 D84 = D84 / 1000 fs = fs / 100 # calculate the area of flow dW = y / tan(b) dP = y / sin(b) A = y * (p + dW) R = A / (p + 2 * dP) # area / perimeter # use Ferguson 2007 to calculate the stream velocity Res = a1 * a2 * (R/D84) / (a1^2 + a2^2 * (R/D84)^(5/3))^(1/2) # use the Keulegan Equation # Res = (1/0.4)*log(12.2*R/(D84)) U = Res * (g * R * S)^(1/2) # use the equations from Knight and others to partition stress SFbank = 10^(-1.4026 * log10(p/(2 * dP) + 1.5) + 0.247) bed_stress = g * rho * y * S * (1 - SFbank) * ((p + 2 * dW) / (2 * p) + 0.5) bank_stress = g * rho * y * S * SFbank * (2 * p + 2 * dW)*sin(b)/(4*y); # estimate the largest grain that the flow can move C = bed_stress / (0.02 * g * rho * Gs) #estimate the division between key stones and bed material load # (this corresponds to the approximate limit of full mobility) K = bed_stress / (0.04 * g * rho * Gs) # use Wilcock and Crowe to estimate the sediment transport rate tau_star_ref = 0.021 + 0.015 * exp (-20 * fs) tau_ref = tau_star_ref * g * Gs * rho * D50 X = bed_stress / tau_ref if(X < 1.35){ W_star = 0.002 * X^(7.5) }else{ W_star = 14 * (1 - (0.894/(X^0.5)))^(4.5) } Qb = p * (W_star / (1.65 * g)) * (bed_stress/rho)^(3/2) #create a list to store data state = list() state$W = p + 2 * dW state$d = A / state$W # area / width state$U = U state$Q = A * U state$Tbed = bed_stress state$Tbank = bank_stress state$Qb = Qb state$C = C * 1000 state$K = K * 1000 state$b = b * 180 /pi # save the value of sideslope and bed width state$p = p return(state) } #################################################################### FindQ = function(Q, p, b, S, D50, D84, fs){ # find the appropriate depth to satisfy continuity with the # specified design flow, Q # guess the initial channel depth, set up variables y = 0.3 * Q ^ 0.3 deltaX = 0.001 * Q ^ 0.3 tol = 0.00001 # first test test = ChannelState(p, y, b, S, D50, D84, fs) converg = (test$Q - Q) / Q while(abs(converg) > tol){ test2 = ChannelState(p,y+deltaX, b, S, D50, D84, fs) M = (test2$Q - test$Q)/deltaX B = (test$Q - Q) - M*y y = -B/M test = ChannelState(p, y, b, S, D50, D84, fs) converg = (test$Q - Q) / Q } return(y) } #################################################################### FindStable = function(Q, mu, p, S, D50, D84, D90, fs){ # find the stable channel shape for the specified Q and # relative bank strength, mu # specify constants and set mu to an equivalent mod_phi phi <<- 40 mod_phi <<- atan(mu * tan(phi * pi / 180)) * 180/pi Tol = 0.00001 deltaX = 0.00001*mod_phi tau_star <<- 0.02 # set the upper and lower angle limits b_upper = mod_phi - deltaX b_lower = deltaX b = 0.67 * mod_phi # calculate the bank stability index y = FindQ(Q,p,b,S,D50,D84,fs) test = ChannelState(p,y,b,S,D50,D84,fs) bank_crit = g*rho*Gs*(D90/1000) * tau_star * (tan(mod_phi*pi/180) / tan(phi*pi/180)) * (1 - (sin(b*pi/180)^2 / sin(mod_phi*pi/180)^2))^0.5 converg = (test$Tbank - bank_crit) / bank_crit while(abs(converg) > Tol){ if(converg > 0){b_upper = b} else {b_lower = b} b = 0.5*(b_upper + b_lower) y = FindQ(Q,p,b,S,D50,D84,fs) test = ChannelState(p,y,b,S,D50,D84,fs) bank_crit = g*rho*Gs*(D90/1000) * tau_star * (tan(mod_phi*pi/180) / tan(phi*pi/180)) * (1 - (sin(b*pi/180)^2 / sin(mod_phi*pi/180)^2))^0.5 converg = (test$Tbank - bank_crit) / bank_crit } return(test) }