clc;clear all;
% 读取降雨P,蒸发皿蒸发EI
P_EI=xlsread('pe.xlsx');
P=P_EI(:,2);EI=P_EI(:,3);
% 读取参数
F=537;%流域面积
WUM=20;WLM=75;WDM=80;%各水层最大蓄水容量
K=0.65;%流域蒸散发能力与蒸发皿蒸发量之比
C=0.11;%深层蒸散发系数
IMP=0;%不透水面积比
B=0.3;%总蓄水容量曲线抛物线指数
SM=20;%流域平均自由水库最大蓄水深度
EX=1;%自由水蓄水容量曲线抛物线指数
KG=0.3;%地下水从自由蓄水库中的出流系数
KSS=0.41;%壤中流出流系数
KKG=0.99;%地下水日退水系数
KKSS=0.6;%壤中流日退水系数
DT=2;%计算时间步长
UH=[0.3 0.6 0.1];
% 计算迭代次数
D=24/DT;%一日内时段数
n=2*D;
% 转换自由水库相关系数时间尺度
KSSD=KSS;
KSS=(1-(1-(KG+KSS))^(1/D))/(1+KG/KSS);%壤中流时段出流系数
KG=KG*KSS/KSSD;%地下水从自由蓄水库中的时段出流系数
% 分配计算空间
EU=zeros(n,1);EL=zeros(n,1);ED=zeros(n,1);
R=zeros(n,1);RG=zeros(n,1);RS=zeros(n,1);RSS=zeros(n,1);
QRG=zeros(n,1);QRSS=zeros(n,1);QRS=zeros(n,1);
% 读取初始状态
WU(1)=0;WL(1)=70;WD(1)=80;
FR(1)=0.1;S(1)=20;
QRS(1)=0;QRSS(1)=40;QRG(1)=20;
% 参数计算
% 土壤最大蓄水容量及深度
WM=WUM+WLM+WDM;
WWMM=WM*(1+B)/(1-IMP);
% 自由水库最大蓄水容量及深度
SSM=(1+EX)*SM;
% 径流计算
% 计算蒸发能力EM
EM=K*EI;
% 时间推进
for t=1:n
W(t)=WU(t)+WL(t)+WD(t);
A(t)=WWMM*(1-(1-W(t)/WM).^(1/(1+B)));
% 计算实际蒸发E
if WU(t)+P(t)>=EM(t)
EU(t)=EM(t);
EL(t)=0;
ED(t)=0;
else
EU(t)=EM(t);
EL(t)=(EM(t)-EU(t))*WL(t)/WLM;
if EL(t)>=C*(EM(t)-EU(t))
EL(t)=C*(EM(t)-EU(t));
ED(t)=0;
else if WL(t)<C*(EM(t)-EU*(t))
EL(t)=WL(t);
ED(t)=C*(EM(t)-EU(t))-EL(t);
end
end
end
E(t)=EU(t)+EL(t)+ED(t);
% 计算产流量R
PE(t)=P(t)-E(t);
if PE(t)<=0
R(t)=0;
else if PE(t)+A(t)>=WWMM
R(t)=PE(t)+W(t)-WM;
else
R(t)=PE(t)+W(t)-WM+WM*(1-(PE(t)+A(t))/WWMM).^(1+B);
end
end
% 划分地面径流RS、壤中流RSS和地下径流RG
AU(t)=SSM*(1-(1-S(t)/SM).^(1/(1+EX)));
FR(t)=1-(1-W(t)/WM).^(B/(1+B));
if PE(t)<=0
R(t)=0;
RS(t)=0;
RSS(t)=S(t)*KSS*FR(t);
RG(t)=S(t)*KG*FR(t);
S(t+1)=(1-KSS-KG)*S(t);
else if PE(t)+AU(t)<SSM
RS(t)=(PE(t)+S(t)-SM+SM*(1-(PE(t)+AU(t))/SSM).^(1+EX))*FR(t);
RSS(t)=(SM-SM*(1-(PE(t)+AU(t))/SSM).^(1+EX))*KSS*FR(t);
RG(t)=(SM-SM*(1-(PE(t)+AU(t))/SSM).^(1+EX))*KG*FR(t);
S(t+1)=(1-KSS-KG)*(SM-SM*(1-(PE(t)+AU(t))/SSM).^(1+EX));
else
RS(t)=(PE(t)+S(t)-SM)*FR(t);
RSS(t)=SM*KSS*FR(t);
RG(t)=SM*KG*FR(t);
S(t+1)=(1-KSS-KG)*SM;
end
DRS(t)=P(t)*IMP;
RS(t)=RS(t)+DRS(t);
end
% 计算时段末土壤蓄水量
inWU(t+1)=WU(t)+PE(t)-R(t);
if inWU(t+1)>WUM
WU(t+1)=WUM;
inWL(t+1)=WL(t)+(inWU(t+1)-WUM);
if inWL(t+1)>WLM
WL(t+1)=WLM;
WD(t+1)=WD(t)+(inWL(t+1)-WLM);
else
WL(t+1)=inWL(t+1);
WD(t+1)=WD(t);
end
else
WU(t+1)=inWU(t+1);
WL(t+1)=WL(t);
WD(t+1)=WD(t);
end
end
% 计算地面径流QRS、壤中流QRSS和地下径流流域出口流量QRG
% 利用径流单位线法计算地面径流
for t=4:n
% 地面径流
QRS(2)=UH(1)*RS(2)*F/(3.6*DT);
QRS(3)=UH(1)*RS(3)*F/(3.6*DT)+UH(2)*RS(2)*F/(3.6*DT);
QRS(t)=UH(1)*RS(t)*F/(3.6*DT)+UH(2)*RS(t-1)*F/(3.6*DT)+UH(3)*RS(t-2)*F/(3.6*DT);
end
% 计算壤中流总入流QRSS和地下径流汇流QRG
for t=1:n
QRSS(t+1)= QRSS(t)*KKSS^(1/D)+RSS(t)*(1-KKSS^(1/D))*F/(3.6*DT);
QRG(t+1)=QRG(t)*KKG^(1/D)+RG(t)*(1-KKG^(1/D))*F/(3.6*DT);
end
% 计算输出总流量Q
for t=1:n
Q(t)=QRS(t)+QRSS(t)+QRG(t);
end
% 绘制流量过程线
hold on;
plot(Q,'B');plot(QRS,'--R');
plot(QRSS,'--G');plot(QRG,'--K');
xlabel('时间 h');ylabel('流量 m^3/s');
legend('总径流','地表径流','壤中流','地下径流');