1 蓄满产流模型原理
1.1流域蒸散发
流域蒸散发在流域水量平衡中起着重要作用。植物截流、地面填洼水量及张力土壤蓄水量的消退都耗于蒸散发，蒸散发计算成果直接影响模型产流计算成果。
在新安江模型中，流域蒸散发计算按土壤垂向分布的不均匀性将土层分为三层，用三层蒸散发模型计算蒸散发量。各参数如表1-1所示。
表1-1 三层蒸散发模型参数及含义
| 参数 | 含义 | 单位 |
| W | 总的张力水蓄量 | mm |
| WU | 上层张力水蓄量 | mm |
| WL | 下层张力水蓄量 | mm |
| WD | 深层张力水蓄量 | mm |
| WM | 流域平均张力水容量 | mm |
| UM | 上层张力水容量 | mm |
| LM | 下层张力水容量 | mm |
| DM | 深层张力水容量 | mm |
| E | 总的蒸散发量 | mm |
| EU | 上层蒸散发量 | mm |
| EL | 下层蒸散发量 | mm |
| ED | 深层蒸散发量 | mm |
| EP | 蒸散发能力 | mm |
| KC | 蒸散发折算系数 | / |
| C | 深层蒸散发扩散系数 | / |
计算公式：
WM+UM+LM+DM
W=WU+WL+WD
E=EU+EL+ED
EP=KC·EM
WU+P≧EP时：
EU=EP，EL=0，ED=0
WU+P<EP，WL≧C·WLM时：
EU=WU+P，EL=（EP -EU）WL/WLM，ED=0
WU+P<EP，C（EP -EU）≦WL≦C·WLM时：
EU=WU+P，EL=C（EP -EU），ED=0
WU+P<EP，WL<C（EP -EU）时：
EU=WU+P，EL=WL，ED=C（EP -EU）- EL
1.2产流计算
根据蓄满产流的概念，采用蓄水容量-面积分配曲线来考虑土壤缺水量分布不均匀的问题，如图1所示，具体计算公式如下：
若PE+A<WMM，即流域部分产流时：
若PE+A≧WMM，即全流域产流时：
在不透水面积上的产流计算如下：
各参数的意义如表1-2所示。
表1-2 产流计算参数及含义
| 参数 | 含义 | 单位 |
| B | 蓄水容量-面积分配曲线的指数 | / |
| WMM | 流域单点最大蓄水量 | mm |
| A | 流域平均的最大初始土壤含水量对应的值 | mm |
| R | 总径流量 | mm |
| IMP | 不透水面积比例 | / |
图1 蓄水容量-面积分配曲线
1.3水源划分
依据当地下垫面和土壤特性，采用三水源的水源划分结构，同时利用自由水蓄水库结构解决水源划分问题。
自由水蓄水库结构考虑了包气带的垂向调蓄作用，如图2所示。按蓄满产流模型计算出的总径流量R，先进入自由水蓄水库调蓄，再划分水源，从图1中可见，产流面积上自由水蓄水库设置了两个出口，一个是旁侧出口，形成壤中流RS；另一个为向下出口，形成地下径流RG。根据蓄满产流的概念，只有在产流面积FR上才可能产生径流，而产流面积是变化的。所以，自由蓄水库的底宽FR也是变化的。
图2 自由水蓄水库结构图
由于饱和坡面流的产流面积是不断变化的，所以在产流面积FR上自由水蓄水容量分布是不均匀的。三水源划分结构是采用类似于流域蓄水容量-面积分配曲线的流域自由水蓄水容量-面积分配曲线来考虑流域内自由水蓄水容量分布不均匀的问题，如图3所示。
计算公式：
当PE≦0时：
当PE>0时：
（1）如果PE+AU<SMM
（2）如果PE+AU≧SMM
各参数的意义如表1-3所示。
表1-3 三水源划分模型参数及含义
| 参数 | 含义 | 单位 |
| S | 流域自由水蓄水量 | mm |
| SM | 流域单点最大自由水蓄水容量 | mm |
| EX | 流域自由水蓄水容量-面积分配曲线的方次 | / |
| AU | 流域初始自由水蓄水量对应的值 | mm |
| FR | 产流面积比例 | / |
| KG | 自由水对地下水的日出流系数 | / |
| KKG | 地下水消退系数 | / |
| KSS | 自由水对壤中流的日出流系数 | / |
| KKSS | 壤中流消退系数 | / |
| FE | 初始土壤含水容量折算系数 | / |
图3 流域自由水蓄水容量-面积分配曲线
2 汇流计算原理
在新安江模型蓄满产流和三水源划分计算的基础上，对流域进行汇流计算。汇流计算又可分为三部分，分别为地表汇流、壤中流汇流和地下径流汇流，不考虑河网汇流。新安江模型计算流程如图4所示。
图4 新安江模型计算流程示意图
2.1地表汇流计算
地表汇流的坡地汇流汇流采用单位线，计算公式：
式中：
QS：地面径流，m3/s；
RS：地面径流，m3/s；
UH：时段单位线，m3/s；
Φ：卷积运算符
2.2壤中流汇流计算
表层自由水侧向流动，出流后成为表层壤中流进入河网。若土层较厚，表层自由水还可以渗入到深层中，经过深层土的调蓄作用才进入河网。壤中流汇流可采用线性水库法模拟，计算公式：
式中：
QSS：壤中流，m3/s；
KKSS：壤中流消退系数；
RSS：壤中流径流量，mm
2.3地下径流汇流计算
地下径流同样采用线性水库法进行计算，计算公式：
式中：
QG：地下径流，m3/s；
KKG：地下径流消退系数；
RG：地下径流量，mm
3 Matlab代码
clc,clear
Q=xlsread('');%读取流量数据
P=xlsread('');%读取雨量数据
E=xlsread('');%读取蒸发数据
KC=0.71;
B=3.0;
n=size(P,1);
EP=E.*KC;
FE=0.8;
WM=126;
WUM=63;
WLM=13;
WDM=50;
SM=36;
KSS=0.06;
KG=0.05;
KKG=0.995;
KKSS=0.83;%自行输入各参数值，名称与参数名称一一对应，上述仅为示例
WMM=WM*(1+B);
WU=zeros(n,1);
WU(1)=WUM*FE;
WL=zeros(n,1);
WL(1)=WLM*FE;
WD=zeros(n,1);
WD(1)=WDM*FE;%初始化各层的土壤含水量
W=zeros(n,1);
S=zeros(n,1);
S(1)=SM*FE;
W(1)=WU(1)+WL(1)+WD(1);
C=0.17;
EX=0.46;
SMM=SM*(1+EX);
IMP=0.001;
S=;%集水面积，单位为km2，上述参数值自行输入
U=S./(3.6*24);
UH=[0,460.5,140.5,67.8,35.2,18.9,10.3,5.7,3.2,1.8,1.0,0.6,0.3,0.2,0.1,0.1,0];%24h10mm单位线，自行定义，这里仅作为示例
for i=1:n
if WU(i)+P(i)>=EP(i)
EU(i)=EP(i);
EL(i)=0;
ED(i)=0;
else
EU(i)=WU(i)+P(i);
if WL(i)>=C.*WLM
EL(i)=(EP(i)-EU(i)).*WL(i)./WLM;
ED(i)=0;
elseif WL(i)<C.*WLM && WL(i)>=C.*(EP(i)-EU(i))
EL(i)=(EP(i)-EU(i)).*C;
ED(i)=0;
else
EL(i)=WL(i);
ED(i)=C.*(EP(i)-EU(i))-EL(i);
end%三层蒸发模型的计算公式
end
a(i)=WMM.*(1-(1-(W(i)./WM)).^(1/(1+B)));%计算A
PE(i)=P(i)-EU(i)-EL(i)-ED(i);%计算净雨
if PE(i)>0
if a(i)+PE(i)<=WMM;
R(i)=PE(i)+W(i)-WM+WM.*((1-(PE(i)+a(i))./WMM).^(1+B));
else
R(i)=PE(i)-(WM-W(i));%产流量计算
end
AU(i)=SMM*(1-(1-S(i)./SM)^(1/(1+EX)));
FR(i)=R(i)./PE(i);
if PE(i)+AU(i)<SMM
RS(i)=(PE(i)-SM+S(i)+SM.*(1-(PE(i)+AU(i))./SMM).^(1+EX)).*FR(i);
RSS(i)=(SM-SM.*(1-(PE(i)+AU(i))./SMM).^(1+EX)).*KSS.*FR(i);
RG(i)=(SM-SM.*(1-(PE(i)+AU(i))./SMM).^(1+EX))*KG.*FR(i);
S(i+1)=(1-KSS-KG)*(SM-SM.*(1-(PE(i)+AU(i))./SMM).^(1+EX));
else
RS(i)=(PE(i)-SM+S(i)).*FR(i);
RSS(i)=SM*KSS.*FR(i);
RG(i)=SM*KG.*FR(i);
S(i+1)=(1-KSS-KG)*SM;%三层水源划分，分别对地面径流、壤中流和地下径流进行计算
end
else
R(i)=0;
FRt(i)=1-(1-W(i)./WM).^(B/(1+B));
FR(i)=R(i)./PE(i);
RS(i)=0;
RSS(i)=S(i).*KSS.*FRt(i);
RG(i)=S(i).*KG.*FRt(i);
S(i+1)=(1-KSS-KG).*S(i);
end
WU(i+1)=P(i)+WU(i)-EU(i)-R(i);
if WU(i+1)>=WUM
WU(i+1)=WUM;
elseif WU(i+1)<=0
WU(i+1)=0;
end
WL(i+1)=WL(i)-EL(i)+P(i)-R(i)-EU(i)-(WU(i+1)-WU(i));
if WL(i+1)>=WLM
WL(i+1)=WLM;
elseif WL(i+1)<=0
WL(i+1)=0;
end
WD(i+1)=WD(i)-EU(i)-EL(i)-ED(i)+P(i)-R(i)-(WU(i+1)-WU(i))-(WL(i+1)-WL(i));
if WD(i+1)>=WDM
WD(i+1)=WDM;
elseif WD(i+1)<=0
WD(i+1)=0;%计算土壤含水量的变化
end
W(i+1)=WU(i+1)+WL(i+1)+WD(i+1);
RIM(i)=(RS(i)+IMP*P(i))./10;%不透水面积上的产流
if i==1
QRSS(i)=(1-KKSS).*RSS(i).*U;
else
QRSS(i)=KKSS.*QRSS(i-1)+(1-KKSS)*RSS(i).*U;
end
if i==1
QRG(i)=(1-KKG).*RG(i).*U;
else
QRG(i)=KKG.*QRG(i-1)+(1-KKG).*RG(i).*U; %三水源的汇流计算
end
end
for i=1:(size(UH,2)-1)
for s=1:n
UUH(i+s,s)=UH(i+1).*RIM(s);
end
end
QS=sum(UUH,2);
QC=((QS(1:n))'+QRSS+QRG)';
E=EU+EL+ED;
Rz=sum(RIM).*10+sum(RSS)+sum(RG);
%水量平衡检验
x1=sum(PE);
x2=W(n+1,:)-W(1)+S(n+1)-S(1);
x3=Rz;
x4=x1-x2;
本文内容仅供参考，如有问题，欢迎批评指正！
标签：WL,mm,预报,WU,Matlab,产流,新安江,EU,EP From： https://blog.csdn.net/JoeyChiu_/article/details/139422009