function Output = UBCRM_base(S, Q, D50, D84, Mu, H, Tstar, RL, TL)
% UBCRM_core is a function to isolate the maximum efficency channel using the
% UBCRMH model where bank strength is represented by modified friction
% angle, varying from the basic friction angle, Phi to just under 90
% degrees
%
%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Q Formative discharge (cumecs)
%
% S Energy gradient (m/m)
%
% D50 Characteristic grain size for sediment transport(mm)
%
% D84 grain size for bed stability analysis and flow resistance (mm)
%
% Mu bank strength expressed relative to the bed entrainment
% threshold (dimensionless) NOTE: if H > 0, then Mu should be set
% to 1.0, so that bank strength is attributed to the cohesive
% term, H
%
% H bank strength as effective rooting depth (m) (optional) NOTE
% that H should be set to 0 if Mu > 1, so that bank strength is
% attributed to a modified friction angle determined by Mu
%
% Tstar critical Shield's number for D84 (optional)
%
% RL flow resistance law selected (optional)
% 1 = Ferguson (2007)
% 2 = Keulegan
% n (you can also simply enter a value of Manning's n)
%
% TL transport law to be selected (optional)
% 1 = Parker (1990)
% 2 = Eaton and Church (2011)
%
%%%% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1 W solution water surface width [m]
%
% 2 d mean hydraulic depth [m]
%
% 3 U mean flow velocity [m]
%
% 4 QB bed material transport rate [me/s]
%
% 5 Qcrit critical flow for entrainment of surface D50 (assuming
% Shields number of 0.06)
%
% 6 bed_str shear stress acting on the bed [Pa]
%
% 7 bank_str shear stres action on the banks [Pa]
%
% 8 thetaSI side slope angle [degrees]
%
% 9 Y height of the trapazoid portion of the channel cross section
%
%
% Copyright 2014 Brett Eaton
%
% Department of Geography, University of British Columbia
%
% Modified: July 2014
%
% This file is part of UBCRM.
%
% UBCRM is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% UBCRM is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Foobar. If not, see .
Phi = 35; %set Phi to default of 35 if no value selected
if nargin < 6
H = 0; %set H to 0 if no cohsive term is entered
end
if nargin < 7 %set Tstar to default of 0.02 if no value selected (following Andrews, 1983)
Tstar = 0.02;
end
ModPhi = atand(Mu*tand(Phi)); %translate Mu to a friction angle
if nargin < 8 %set RL to default of 1 if no value selected
RL = 1;
end
if nargin < 9 %set TL to default of 1 if no value selected
TL = 1;
end
%%%%%%%%%% STEP 1: CHOOSE AN INITIAL VALUE AND PERFORM INITIAL CALCS %%%%%
Ptest = 3.0*Q^0.5; % use hydraulic geometry equation to guess the first bed width
% calculate the local gradient in sediment transport at Ptest
Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
dQ1(1) = Y_plus - Y_minus;
% now move in direction of gradient towards solution
if dQ1 > 0
Ptest = Ptest +0.25*Ptest;
else
Ptest = Ptest - 0.25*Ptest;
end
% perform a second calculation
Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H,ModPhi, Phi, Tstar, RL, TL);
Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
dQ1(2) = Y_plus - Y_minus;
%%%%%%%%%% STEP 2: SET UP A WHILE LOOP TO QUICKLY FIND THE MAXIMUM %%%%%%%
i=2;
while dQ1(i-1)/dQ1(i) > 0 % continue loop until we get a reversal of the signs
if dQ1(i) > 0
Ptest = Ptest +0.25*Ptest;
else
Ptest = Ptest - 0.25*Ptest;
end
i = i + 1;
Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
dQ1(i) = Y_plus - Y_minus;
end
%%%%%%%%%%%% STEP 3: NOW THAT WE HAVE IDENTIFIED THE BOUNDS WITHIN WHICH
%%%%%%%%%%%% THE SOLUTION LIES, WE SET THE BOUNDS AND FIND THE SOLUTION
% set the gates where P_upper is above the solution and P_lower is below it
if dQ1(i) < 0
P_upper = Ptest;
P_lower = Ptest - 0.25*Ptest;%move back one step to other side of max
else
P_lower = Ptest;
P_upper = Ptest + 0.25*Ptest; %move forward one step to other side
end
% seek the maximum efficiency solution using Newtonian method
Tol = 0.001; %set the convergence toleranace, as a proportion of the solution bottom width
while abs((P_upper - P_lower)/Ptest) > Tol
%calculate the transport rate for a stable channel on either side of
%the test bed width
Y_minus = FindTransp(Ptest*0.99999, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
Y_plus = FindTransp(Ptest*1.00001, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
dQ2 = Y_plus - Y_minus;
%MOVE THE GATES BASED ON GRADIENT dQ/dP
if dQ2 < 0
P_upper = Ptest;
else
P_lower = Ptest;
end
%recalculate the test bottom width
Ptest = 0.5*(P_upper + P_lower);
end
%%%%%%%%%%% STEP 4: CALCULATE THE FINAL SOLUTION %%%%%%%%%%%%%%%%%%%%%%%%
Pbed = Ptest; %record the solution bed width
[QB, thetaSI] = ... %extract the solution transport rate and geometry
FindTransp(Pbed, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL);
[~, Y, bed_str, bank_str] = ... %extract trap depth and stresses
FindSS(thetaSI, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL);
[~, W, d, U, ~, ~] = ... %extract channel dimensions
FindQ(Y, H, Pbed, thetaSI, S, D84, RL);
Ycrit = 1.65 .* 0.06 .* (D50/1000) ./ S; %estimate the flow at which the bed becomes entrained
if Ycrit > Y
Hcrit = Ycrit - Y;
Ycrit = Y;
else
Hcrit = 0;
end
[Qcrit] = ... %extract channel dimensions
FindQ(Ycrit, Hcrit, Pbed, thetaSI, S, D84, RL);
%generate output including shear stress calculations
Output = [W, d, U, QB, Qcrit, bed_str, bank_str, thetaSI, Y];
end
%%%%%%%%%%%%%%% END OF MAIN FUNCTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% SUB FUNCTIONS LISTED BELOW %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% SUB FUNCTION FINDTRANSP, CALLED BY MAIN %%%%%%%%%%%%%%%%%
function [QB, thetaSI] = FindTransp(Pbed, S, Q, D50, D84, H, ModPhi, Phi, Tstar, RL, TL)
%FindTransp.m is a function that takes a given channel width, slope,
%discharge, bed textures, and bank strength parameters and then solves for
%the bank angle for a critically stable bank, then calculates the sediment
%transport rate associated with that geometry. The function calls FindSS.m
%to calculate the stable channel geometry, which in turn calls FindQ.m to
%ensure that all geometries considered are consistent with continuity
%
%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pbed the bottom width of the trapezoidal channel [m]
% S energy gradient [m/m]
% Q discharge carried by the channel [m3/s]
% D50 grain size used for sediment transport calculation [mm]
% D84 grain size used for flow resistance and stability analysis [mm]
% H bank strength parameter, expressed as effective rooting depth [m]
% (ranges from phi to <90 degrees)
% Phi friction angle for bank and bed sediment
% Tstar Critical Shields number for D84
% RL choice of flow resistance laws (see FindQ.m for details)
%%%% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% QB transport rate for specified conditions [m3/s]
% thetaSI angle of the channel banks at critical stability [degrees]
% YQ the trapezoid depth associated with specified Q [m]
%determine stable bank angle for the given bed width associated with continuity
Trange = [0.001 ModPhi - 0.001];
%Trange = 0.75*Phi; %test alternate solving routine
options = optimset('TolX', 10e-5);
thetaSI = fzero(@(x) FindSS(x, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL) - 1.0,...
Trange,options);
%extract the solution shear stress values and critical shear stresses
[~, Y, bed_str] = FindSS(thetaSI, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL);
%calculate the sediment transport rate using a sediment transport law
%specified below
[~,~,~,U] = FindQ(Y, H, Pbed, thetaSI, S, D84, RL);
QB = TranspRate(D50, bed_str, U, S, Pbed, TL);
end
%%%%%%%%%%%%%%% SUB FUNCTION FINDSS, CALLED BY MAIN, FINDTRANSP %%%%%%%%%%%
function [SI, YQ, bed_str, bank_str, bed_crit, bank_crit] = FindSS(theta, Pbed, S, Q, D84, H, ModPhi, Phi, Tstar, RL)
%FindSS is a function that calculates the ratio of the critical bank
%strength and the shear stress acting onthe bank for a trapezoidal channel
%with known bed width, slope and discharge
%
%%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% theta side slope angle [degrees]
% Pbed the bottom width of the trapezoidal channel [m]
% S energy gradient [m/m]
% Q discharge carried by the channel [m3/s]
% D84 grain size used for flow resistance and stability analysis [mm]
% H bank strength parameter, expressed as effective rooting depth
% [m]
% (ranges from phi to <90 degrees)
% Phi friction angle for bank and bed sediment
% Tstar Critical Shields number for D84
% RL choice of flow resistance laws (see FindQ.m for details)
%%%% OUTPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SI bank stability index (critical/applied stress
% YQ the trapezoid depth associated with S1 [m]
% bed_str shear stress acting on the bank [Pa]
% bank_str shear stress acting on the bed [Pa]
% bed_crit critical shear stress for bed entrainment [Pa]
%
%%%%%%%%%%%%%%% determine depth associated with continuity %%%%%%%%%%%%%%%%
% Yrange = [(0.0001*Q^0.5 - H) (5.0*Q^0.5)]; %allow Y to approach -H,
% options = optimset('TolX', 10e-5);
% YQ = fzero(@(x) FindQ(x, H, Pbed, theta, S, D84, RL) - Q, Yrange, options);
%this sets the range of depths within which the solution lies
%try using the bisection method
% Yu = 100;
% Yl = - H;
% Yt = (Yu + Yl) / 2;
%
% TestQt = (FindQ(Yt, H, Pbed, theta, S, D84, RL) - Q)/Q;
% Tol = 10e-6;
%
% while abs(TestQt) > Tol
% if TestQt > 0
% Yu = Yt;
% elseif TestQt < 0
% Yl = Yt;
% end
% Yt = (Yu + Yl) / 2;
% TestQt = (FindQ(Yt, H, Pbed, theta, S, D84, RL) - Q)/Q;
% end
% YQ = Yt;
%try using Newton Raphson method
Yt = 0.33*Q^(1/3);
Tol = 10e-6;
Qt = FindQ(Yt, H, Pbed, theta, S, D84, RL);
TestQt = (Qt - Q) / Q;
while abs(TestQt) > Tol
GradQ = (FindQ(1.0001*Yt, H, Pbed, theta, S, D84, RL) - ...
FindQ(0.9999*Yt, H, Pbed, theta, S, D84, RL))...
/(1.0001*Yt - 0.9999*Yt);
Intercept = Qt - GradQ*Yt; %estimate the linear equation at Yt,Qt
Yt = (Q-Intercept)/GradQ; %use the equation to recalculate Yt
Qt = FindQ(Yt, H, Pbed, theta, S, D84, RL); %check the Q value there
TestQt = (Qt - Q) / Q;
end
YQ = Yt;
%%%%%%%%%%%%%%% partitiion the shear stress %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%dimensions when YQ <= 0 are based on a rectangle
if YQ <= 0
Ht = H + YQ;
Pbank = 2 * (Ht);
W = Pbed;
end
%dimensions when YQ > 0 are a composite trapezoid and rectangle
if YQ > 0
Ht = H + YQ;
Pbank = 2*H + 2*YQ/sin(theta*pi/180);
W = Pbed + 2*YQ/tan(theta*pi/180);
end
totstress = 9.8*1000*(Ht)*S;
SFbank = 10^(-1.4026*log10(Pbed/Pbank + 1.5) + 0.247);
%partioning equation
bed_str = totstress*(1 - SFbank)*(W/(2*Pbed) + 0.5);
%stress acting on the bed
bank_str = totstress*SFbank*(W + Pbed)*sin(theta*pi/180)/(4*(Ht));
%stress acting on the bank
%this is a set of empirical equations derrived by Knight and others
%%%%%%%%%%%%%%% calculate the critical shear stress values %%%%%%%%%%%%%%
bed_crit = Tstar*9.8*1650*(D84/1000);
bank_crit = 9.8*1650*(D84/1000)*...
Tstar*(tan(ModPhi*pi/180) / tan(Phi*pi/180))*...
(1 - (sin(theta*pi/180)^2 / sin(ModPhi*pi/180)^2))^0.5;
%%%%%%%%%%%%%%% calculate a bank stability index %%%%%%%%%%%%%%%%%%%%%%%%
% if YQ <= 0 %in this case, root cohesion dominates the bank strength
%
% SI = 1; %and we therefore set SI = 1, regardless of theta
%
% else %if YQ > 0, then proceed with calculation, per normal
SI = bank_crit/ bank_str;
% end
end
%%%%%%%%%%%%%%% SUB FUNCTION FINDQ, CALLED BY MAIN, FINDSS %%%%%%%%%%%%%%%
function [Q, W, d, U, P, R] = FindQ(Y, H, Pbed, Theta, S, D84, RL)
%FindQ is a function to solve for the discharge associated with a set of
%input geometry parameters, a bed surface grain size, and a choice of flow
%resistance law to use, all assuming a trapezoidal channel geometry
%
%%% INPUTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Y flow depth [m]
% P bottom width [m]
% Theta angle of the channel banks [degrees]
% S energy gradient [m/m]
% RL choice of resistance law to use
% 1 = Ferguson (2007)
% 2 = Keulegan
%
%%%% OUTPUT = Q [m3/s] %%%%%%%%%%%%%%%%%%%%%%%%
%apply default flow resistance law
if nargin < 7
RL = 1;
end
%convert to useful units
Theta = Theta * pi / 180; %convert to radians
D84 = D84 / 1000; %convert to m
%define constants
g = 9.81;
a1 = 5.5; %constant for RL = 1
a2 = 2.2; %constant for RL = 1
%calculate geometric parameters
%when Y is negative, then the cross section is rectangular
if Y <= 0
dW = 0;
Area = (H+Y)*Pbed; %cross sectional area of flow
P = Pbed + 2*(H+Y); %wetted perimeter
R = Area / P; %hydraulic radius
%when Y is positive, add a trapezoid to the dimensions of the upper
%rectangle of height H
elseif Y > 0
dW = Y / tan(Theta); %width of the bank at the base
dP = sqrt(Y^2 + dW^2); %length of bank
Area = Y*(Pbed + dW)... %cross section of lower trapezod
+ H*(Pbed + 2*dW); %cross sectional of upper rectangle
P = Pbed + 2*dP + 2*H; %wetted perimeter
R = Area / P; %hydraulic radius
end
%choose the appropriate flow resistance law
if RL == 1
%Use Ferguson's continuously varying power law
Res = a1*a2*(R/D84) / (a1^2 + a2^2*(R/D84)^(5/3))^(1/2);
elseif RL ==2
%Keulegan equation
Res = (1/0.4)*log(12.2*R/(D84));
else
%here is a backdoor to use a value of n for RL
Res = R^(1/6) / (RL * sqrt(g));
end
%calculate velocity, U and discharge, Q
U = Res*sqrt(g*R*S);
Q = Area * U;
W = Pbed + 2*dW;
d = Area / W;
end
%%%%%%%%%%%%%%% SUB FUNCTION TRANSPRATE, CALLED BY FINDTRANSP %%%%%%%%%%%%
function Qb = TranspRate(D50, tau, U, S, W, TL)
%D50 the surface D50, used as the characteristic grain size
%tau the bed shear stress
%U the mean flow velocity, used to estimate stream power
%W the width of the sediment transport zone
if TL == 1; %calculate sediment transport using Parker 1990 equation
%define constants used in the equation
rho = 1000;
rho_sed = 2650;
g = 9.81;
% calculate dimensionless shear stress
Tstar = tau/(g*(rho_sed-rho)*(D50/1000));
Tref=0.0386;
phi=Tstar/Tref;
%determine which equation to use
if phi > 1.59
G=5474*(1-0.853/phi)^4.5;
%t=1
elseif phi < 1
G=phi^14.2;
%t=2
else
G=exp(14.2*(phi-1)-9.28*(phi-1).^2);
%t=3
end
%Calculate the transport rate
Qb=W*G*(0.0025*(tau/rho)^(3/2)/(g*(rho_sed/rho-1)));
end
if TL == 2; %calculate using Eaton and Church 2011 relation
%define constants used in the equation
rho = 1000;
rho_sed = 2650;
g = 9.81;
R = (rho_sed - rho)/rho;
Tcrit = 0.045; %this is the critical shields number
%calculate stream power
omega = tau * U;
om_star = omega / (rho*(g*R*(D50/1000))^(3/2));
dcrit = Tcrit * R * (D50/1000) / S;
%Use Ferguson's continuously varying power law
a1 = 6.5; %constant for RL = 1
a2 = 2.5; %constant for RL = 1
Res = a1*a2*(dcrit/(D50/1000)) / (a1^2 + a2^2*(dcrit/(D50/1000))^(5/3))^(1/2);
om_crit = Res * Tcrit^(3/2);
%calculate the transport rate
E_star = (0.92 - 0.25*sqrt(om_crit/om_star))^9;
qb_star = E_star * om_star;
Qb = W * qb_star * sqrt(R*g)*(D50/1000)^(3/2);
end
end