水文预报课程设计.doc

《水文预报》
课程设计报告

学 院：_____水利与环境学院_____
专 业：____水文与水资源工程____
班 级： 200905201
姓 名：________马天玉__________
学 号：______20090520115___
指导教师：________胡彩虹________
精品资料
第一章 基本任务
1.1蒸发折算系数Kc的优选
根据已给数据资料及参数（本报告采用89-92年的历史数据），将流域作为整体：
（1）进行日模型产流量计算；
（2）比较计算年径流与实测年径流；
（3）通过误差分析，优选蒸发折算系数Kc；
（4）89～90年的历时数据作为率定参数，91～92年的数据作为模型检验。
1.2暴雨预报
根据已给的设计暴雨资料和任务一率定的Kc，将流域作为整体进行如下计算：
（1）次洪产流量计算，划分水源；
（2）直接径流汇流，地下径流汇流的计算。
（3）采用2004年暴雨数据进行预报。
根据已给的资料、参数及做过的习题，自己编写程序，将流域作为整体进行产流量计算；将计算年径流与实测年径流进行比较。

第二章 基本资料
2.1流域概况
白盆珠水库位于广东省东江一级支流西枝江的上游，坝址以上集雨面积856 km2。流域地处粤东沿海的西部，海洋性气候显著，气候温和，雨量丰沛。暴雨成因主要是锋面雨和台风雨，常受热带风暴影响。降雨年际间变化大，年内分配不均，多年平均降雨量为1800mm，实测年最大降雨量为3417mm，汛期4～9月降雨量占年降雨量的81％左右：径流系数0.5～0.7。
流域内地势平缓，土壤主要有黄壤和砂壤，具有明显的腐殖层，淀积层和母质土等层次结构，透水性好。台地、丘陵多生长松、杉、樟等高大乔木；平原则以种植农作物和经济作物为主，植被良好。
流域上游有宝口水文站，流域面积553km2，占白盆珠水库坝址以上集雨面积的64.6%。白盆珠水库有10年逐日入库流量资料、逐日蒸发资料和时段入库流量资料：流域内有7个雨量站，其中宝口以上有4个。雨量站分布较均匀．有10年逐日降水资料和时段降水资料；宝口水文站具有10年以上水位、流量资料；流域属山区性小流域且受到地形、地貌等下垫面条件影响，洪水陡涨缓落，汇流时间一般2—3h，有时更短；一次洪水总历时2～5d。

图2-1 该地区水文站分布图

2.2数据资料
(1) 计算流域面积为553km2。
(2) 流域内有四个雨量站，权重系数分别为0.33、0.14、0.33、0.20。
(3) 日产流模型数据，具体见资料文件名：87-92data.xls，数据格式为：
T(i) Q(i) E(i) P1(i) P2(i) P3(i) P4(i)
(4)暴雨预报的数据，见表2-1
表2-1 2004年暴雨过程数据表
时间
蒸散发(mm)
降雨量（mm）
禾多布
马山
高潭
宝口
T（i）
E（i）
P1（i）
P2（i）
P3（i）
P4（i）
2004-9-23 12:00
1.3
6.2
9.9
21.6
17.3
2004-9-23 15:00
1.3
7.6
16
20.6
12.6
2004-9-23 18:00
1.3
6.2
6.4
14.9
15.9
2004-9-23 21:00
1.3
8.8
17.2
29.4
18.5
2004-9-24 0:00
1.2
25
34.8
35.3
24.6
2004-9-24 3:00
0.9
29.9
29.2
43.9
37.8
2004-9-24 6:00
0.9
38.6
24.8
46.9
33
2004-9-24 9:00
0.9
6.9
7.5
6.1
12.3
2004-9-24 12:00
0.9
28.3
29.9
34.2
28.5
2004-9-24 15:00
0.9
25.6
42.7
39.8
75.4
2004-9-24 18:00
0.9
93.9
137.6
124
13.2
2004-9-24 21:00
0.9
85.3
90.8
85
75.9
2004-9-25 0:00
0.8
51.5
47.7
49.2
38.5
2004-9-25 3:00
1.1
39.8
70.3
42.1
97.7
2004-9-25 6:00
1.1
43.2
47.3
61.5
45.9
2004-9-25 9:00
1.1
20.5
13.3
15.8
13.1
2004-9-25 12:00
1.1
10.5
8
1.8
3.3
2004-9-25 15:00
1.1
7.4
8.4
7.6
10.9
2004-9-25 18:00
1.1
1.8
2.8
2.1
4.6
2004-9-25 21:00
1.1
0.2
0
0.3
0
2004-9-26 0:00
1.2
0
0
0
0
2004-9-26 3:00
2.1
0
0
0
0
2004-9-26 6:00
2.1
0
0
0
0
2004-9-26 9:00
2.1
0
0
0
0
2004-9-26 12:00
2.1
0
0
0
0
2004-9-26 15:00
2.1
0
0
0
0
2004-9-26 18:00
2.1
0
0
0
0
2004-9-26 21:00
2
0
0
0
0

（5）计算参数数据，见表2-2
表2-2 计算参数表
计算年份
参数
初始张力水蓄量
1989—1990
Wm
Um
Lm
Dm
W
WU
WL
WD
140
20
60
60
110
10
40
60
B
C
Fc
IM
0.2
0.16
22
0.001
（6）流域单位线
单位线过程（m3/s）为：0，40，80，130，100，80，48，20，10，5，0
（7）地下径流汇流
Cg=0.978，Qg=55.3m3/s

第三章 计算公式
该流域海洋性气候显著、气候温和、雨量丰沛，多年平均降雨量为1800mm，径流系数0.5-0.7，土壤主要有黄壤和砂壤，层次结构明显，透水性好，植被覆盖度高，地势平坦，由此可初步判定该流域的产流机制为蓄满产流模式。
3.1产流计算
3.1.1蒸散发计算
根据流域蓄满产流特点，蒸散发计算采用的是三层蒸散发计算模式。三层蒸发模式的具体计算如下：
1)当WU+P≥EP,
EU=Ep,EL=0,ED=0;
2)当WU+P<EP, WL≥C×WLM,
EU=WU+P,EL=(EP-EU)×WL/WLM,ED=0;
3)当WU+P<EP, C(EP-EU)≤ WL<C×WLM,
EU=WU+P,EL=C(EP-EU),ED=0;
4)当WU+P<EP，WL<C(EP-EU)，
EU=WU+P,EL=WL,ED=C(EP-EU)-EL.
R
ΔW
图2-1 包气带蓄水容量曲线

式中：WU为上层土壤蓄水量，WL为下层土壤蓄水量，EU为上层土壤蒸发量，EL为下层土壤蒸发量，ED为深层土壤蒸发量，P为流域平均降雨量，Ep为流域平均蒸发能力，C为深层蒸散发扩散系数，WLM为下层张力水蓄水容量。
3.1.2产流量计算
根据流域特点，产流量计算系根据蓄满产流理论得出的。蓄满产流，即任一地点上，土壤含水量达田间持水量前，降雨量全部补充土壤含水量，不产流；当土壤蓄满后，其后续降雨量全部产生径流。流域内各点包气带的蓄水容量是不同的，将各点包气带蓄水容量从小到大排列，以包气带达到田间持水量时的土壤含水量WM′为纵坐标，以流域内小于等于该WM′的面积占全流域的面积比α为横坐标，所绘的曲线称为流域蓄水容量曲线。
a=WMM× (1-(1-W/WM) 1/(b+1)
PE>0，则产流；否则不产流。产流时：
1)当PE+a≤WMM：
R=PE+W-WM+WM×(1-(PE+a)/WMM)b+1
2)当PE+a>WMM：
R=PE+W-WM
式中：PE为扣除蒸发量后的降雨量，a为土壤含水量W对应的土壤水深，WM为流域平均蓄水容量，WMM为流域各地点包气带蓄水容量的最大值，b为流域包气带蓄水容量分布的不均匀指数，R为流域产流量。
3.1.3二水源划分
流域坡地上的降雨产流量因产流过程的条件和运动路径不同，受流域的调蓄作用不同，各径流成分在流量过程线上的反应是不一样的。在实际工作中，常需按各种径流成分分别计算或模拟，因为要对产流量进行水源划分。
直接径流和地下径流水源划分如下：
1）当PE<=FC时：RS=0.0
RG=R
2）当PE＞FC时：RG=FC*R/PE
RS=R-RG
式中：FC为稳定下渗率，RS为直接径流，RG为地下径流。
3.1.4各层蓄水量计算
降雨补充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，补充来源为降雨减去蒸散发减去径流量，顺序为上、下、深层依此补充。
三层蓄水量变化的具体计算如下：
1）WU[i]+P[i]-EU[i]-R[i]<=UM，
WU[i+1]=WU[i]+P[i]-EU[i]-R[i];
WL[i+1]=WL[i]-EL[i];
WD[i+1]=WD[i]-ED[i];
2）WU[i]+P[i]-EU[i]-R[i]>UM，WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM)<=LM，
WU[i+1]=UM;
WL[i+1]=WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM);
WD[i+1]=WD[i]-ED[i];
3)WU[i]+P[i]-EU[i]-R[i]>UM，WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM)>LM，
WD[i]-ED[i]+WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM)-LM<=DM时，
WU[i+1]=UM;
WL[i+1]=LM;
WD[i+1]=WD[i]-ED[i]+WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM)-LM;
4)WU[i]+P[i]-EU[i]-R[i]>UM，WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM)>LM， WD[i]-ED[i]+WL[i]-EL[i]+(WU[i]+P[i]-EU[i]-R[i]-UM)-LM>DM时，
WU[i+1]=UM;
WL[i+1]=LM;
WD[i+1]=DM;
W[i+1]=WU[i+1]+WL[i+1]+WD[i+1]；
式中：i表示第i天。
3.2汇流计算
根据流域净雨和流域径流单位线，采用卷积的差分形式算出流域出口的流量过程。直接径流汇流可根据该流域的时段单位线推求，地下径流汇流由线性水库演算法推求。
1）直接径流汇流计算公式：
QS(i)=RS(i)×UH；
式中：UH为该流域的单位线
2）地下径流汇流计算公式：
QG(i)=CG×QG(i-1)+(1-CG)×RG(i)×U
Q(i)=QS(i)+QG(i)
式中：U为单位转换系数，
3.3模型参数
1）Kc：蒸散发能力折算系数，它主要反映流域平均高程与蒸发站高程之间差别的影响和蒸发皿散发与路面蒸散发间差别的影响；
2）WM:流域平均张力水容量，它表示流域蓄满的标准；
3）WUM:上层张力水蓄水容量，它包括了植物截留量；
4）WLM：下层张力水蓄水容量；
5）b：流域包气带蓄水容量分布的不均匀指数，在一般情况下其取值与单元流域面积有关；
6）C:深层蒸散发扩散系数，它主要取决于流域内深根植物的覆盖范围。
7）IM:不透水面积占全流域面积的比例，它的值可由大比例尺的地形图，通过地理信息系统现代技术量测出来，也可用历史上干旱期小洪水资料来分析。

第四章 模型检验、结果评定及分析
水文预报是一项直接服务于国家安全和国民经济建设的不可或缺的重要基本工作，是帮助人类有效地预防洪水、减少洪灾损失，有效利用水资源的非工程措施之一。随着经济、社会发展及其全球化进程的需要，水文预报的服务面进一步拓展，对水文预报提出了更高的要求。
水文预报结果的准确率与可信程度是衡量服务质量的前提，为了更好地为国家安全和国民经济建设服务，必须对水文预报结果的可靠性和有效性进行评定和检验。
4.1产流模式的检验
定性分析
该流域集雨面积856km2。流域地处南方，海洋性气候显著，气候温和，雨量丰沛。暴雨成因主要是锋面雨和台风雨，常受热带风暴影响。降雨年际间变化大，年内分配不均，多年平均降雨量为1800mm，实测年最大降雨量为3417mm，汛期4—9月降雨量占年降雨量的81％左右：径流系数0.5~0.7。流域内地势平缓，土壤主要有黄壤和砂壤，具有明显的腐殖层，淀积层和母质土等层次结构，透水性好。台地、丘陵多生长松、杉、樟等高大乔木；平原则以种植农作物和经济作物为主，植被良好。
流域上游有一水文站，控制流域面积553km2，占流域集雨面积的64.6%。该水文站以上有4个雨量站。雨量站分布较均匀，有10年逐日降水资料和时段降水资料；该水文站具有10年以上水位、流量资料；流域属山区性小流域且受到地形、地貌等下垫面条件影响，洪水陡涨缓落，汇流时间一般2—3h，有时更短；一次洪水总历时2~5d。
由流域概况可以看出，该地区属于湿润多雨地区，雨量集中，地势平坦，土壤层容易蓄满，而且有场系列的降雨资料和水位流量资料，综合从气候条件、土壤状况、植被组成以及洪水机制看，新安江模型适用于该地区的水文预报。
4.2率定期优选蒸发折算系数Kc
4.2.1计算年径流与实测年径流的比较结果及误差分析
采用试算法，设置Kc取值在0.9-4.9之间，步长为0.001，分别用1987-1992年的资料数据进行计算，得到各年全年径流量理论计算值和实测值的相对误差值随着Kc的变化过程如下图1：
图1 各年全年径流量理论计算值和实测值的相对误差值随着Kc的变化

自左至右分别为1991、1988、1992、1989、1987、1990年，各年相对误差最小时，Kc取值如下表：
表1 各年Kc最佳取值
年份
Kc最佳取值
相对误差
1987
1.483
0.00013336
1988
1.227
0.000295935
1989
1.423
5.32E-05
1990
2.565
0.000136318
1991
1.082
0.000106609
1992
1.418
0.000149347
可以看出，1990年的资料不可取，因此舍弃不用。
4.2.2 蒸发折算系数Kc值的优选方法
在流域整体的计算径流量与实测径流量的相对误差满足5%以内的前提下，尽可能使1987-1989年连续三年的相对误差均在10%以内，并使得各年的相对误差差别尽可能小。根据相对误差规律得以下三种方法，同时说明：
①通过人为观察、比较，人工选取所给Kc的最优值；
②计算1987-1989年三年的相对误差累计值，相对误差累计值在一定程度上反映了Kc对计算径流量造成的影响，相对误差累计值越小越好；
③计算1987-1989年三年相对误差的均方差，它在一定程度上说明了各相对误差偏离平均相对误差的程度，也反映了Kc取值造成误差的稳定程度，相对误差的均方差越小越好。
4.2.3确定Kc的取值：
缩小Kc的取值范围，取Kc=0.9-1.3，步长为0.001进行计算，得到各年全年径流量理论计算值和实测值的相对误差值随着Kc的变化过程如下图2.
图2 各年全年径流量理论计算值和实测值的相对误差值随着Kc的变化

表2 各年Kc最佳取值

年份
Kc最佳取值
相对误差
1987
1.3
0.074130117
1988
1.227
0.000295935
1989
1.3
0.068792805
1991
1.082
0.000106609
1992
1.3
0.049825263

因为1990年数据舍弃，原定4年率定，改为采用1987-1989年资料进行率定，得到KC取1.263-1.430之间数据时，误差均在5%以内，其中，Kc取值为1.3450时，误差最小。
4.3 模型检验及评价
4.3.1 模型检验结果
通过对Kc取值在1.263-1.430之间进行变化绘制日径流量理论计算值和实测值的对比，得到Kc取值为1.300时，拟合较好，此时各年相对误差如下表3.
表3 Kc=1.300时的各年份相对误差值
年份
87
88
89
90
91
92
相对误差
0.0741
-0.0169
0.0983
0.5440
-0.0754
0.0715

4.3.2 1987-1992年各年计算径流与实测径流的拟合结果
Kc取1.300时，各年逐日径流理论计算值和实测值对比图见下图3-8.

图3 1987年逐日径流理论计算值和实测值对比图

图4 1988年逐日径流理论计算值和实测值对比图

图5 1989年逐日径流理论计算值和实测值对比图

图6 1990年逐日径流理论计算值和实测值对比图

图7 1991年逐日径流理论计算值和实测值对比图

图8 1992年逐日径流理论计算值和实测值对比图

以上各图是在日径流理论计算时未计算基流的情况下得到的，可见其整个起伏趋势还是相当契合的，但局部仍存在差异，基本可以满足常次预报的要求。因此确定Kc取1.300。
4.4误差来源
设计的蓄满产流模型结构与流域的实际产流过程和规律不完全相符， 出现的问题以及可能误差影响因素包括：
（1）1990年数据计算出的结果误差过大，可能是因为1990年数据存在问题。
（2）各年先对误差均已经很小，但不能全部达到5%的范围之内，可能是因为调试不够，或者是因为模型假设与市级的情况不尽相同。
（3）逐日径流计算理论值和实测值之间存在较大的相对误差，可能是因为未对基流进行计算，模型本身精度也有限，导致这种误差.。
（4）由于流域地理、气候、气象、水文条件上与模型假设条件存在一定程度上的差异，导致计算结果存在误差。
4.5模型的应用-暴雨预报
对2004年暴雨过程进行洪水预报。
运用单位线法和出流系数法分别计算直接径流出流量和地下径流出流量，两者之和即为总的流量，其中Kc取1.300。具体结果见下表1。
表4 2004年暴雨过程
时间

月
日
时
R
Rs
Rg
地下径流Qg
直接径流Qs
总径流Q
9
23
12
3.70
2.73
0.97
55.30
0.00
55.30

15
4.30
3.18
1.12
55.35
10.92
66.27

18
3.64
2.35
1.28
55.57
34.58
90.15

21
8.44
6.87
1.57
56.11
70.38
126.49

24
24.95
22.21
2.74
57.96
115.01
172.98

24
3
34.24
31.24
3.00
60.07
228.10
288.17

6
36.53
33.53
3.00
62.13
454.13
516.26

9
6.04
3.04
3.00
64.14
781.11
845.25

12
28.75
25.75
3.00
66.11
984.05
1050.16

15
40.88
37.88
3.00
68.03
1095.58
1163.62

18
92.05
89.05
3.00
69.91
1106.62
1176.53

21
82.33
79.33
3.00
71.76
1495.12
1566.88

24
46.05
43.05
3.00
73.56
2053.13
2126.69

25
3
54.26
51.26
3.00
75.32
2673.36
2748.68

6
48.20
45.20
3.00
77.04
2953.15
3030.19

9
14.31
11.31
3.00
78.72
2909.37
2988.10

12
3.69
0.69
3.00
80.37
2668.85
2749.23

15
6.16
3.16
3.00
81.98
2147.52
2229.50

18
0.45
0.00
0.45
80.69
1500.58
1581.27

21
0.00
0.00
0.00
78.91
964.96
1043.87

24
0.00
0.00
0.00
77.17
540.63
617.81

26
3
0.00
0.00
0.00
75.48
254.57
330.05

6
0.00
0.00
0.00
73.82
122.02
195.84

9
0.00
0.00
0.00
72.19
50.44
122.64

12
0.00
0.00
0.00
70.60
12.66
83.26

15
0.00
0.00
0.00
69.05
3.50
72.55

18
0.00
0.00
0.00
67.53
1.58
69.11

21
0.00
0.00
0.00
66.05
0.00
66.05

24

64.59
0.00
64.59

27
3

63.17
0.00
63.17

6

61.78
0.00
61.78

9

60.42
0.00
60.42

12

59.09
0.00
59.09

15

57.79
0.00
57.79

18

56.52
0.00
56.52

21

55.28
0.00
55.28

24

54.06
0.00
54.06

28
3

52.87
0.00
52.87

将直接径流、地下径流出流过程和总流量过程绘制出洪水流量过程线，见下图9.
图9 洪水流量过程线

第五章 总结和心得
此次课程设计，做了很久，期间碰到诸多棘手的问题，终于一一解决，得以完成，虽然还有些不尽人意，但总体上还是较好得完成了此次课程设计的各项要求。通过此次课程设计，收获良多。首先是对于水文预报这门课程，为了完成这次课设，对课本进行了深入系统的复习，尤其是与新安江模型相关的内容，使得对这门课程的掌握更加扎实牢固，理解也更加深入；其二，本次课设采用的是matlab软件进行数据处理，这也是基于数据较多，计算分析比较复杂的缘故，利用原有的一点儿编程基础，在这期间进一步不断地学习，对该软件的使用能力有了很大程度的增强，这对以后的学习工作都是大有裨益的。这次课设的意义就是在于不断逼迫自己去学习更多的新东西，并把以前学过的东西进行整合贯通，达到提升自身水平的效果；其三，在做课程设计的过程中，老师对我们进行了多次耐心认真的辅导，同学们之间也不断进行互相的交流，也得到了很多的乐趣。
总的来说，这次课程设计受益匪浅。这是毕业设计之前的最后一次课设，也是由书本上知识运用到实践中的一次尝试，让我们意识到自身知识的匮乏和有限，长叹“书到用时方恨少”，同时这也是一个有趣的过程。试想，通过自己的不懈努力终于把它攻克，把成果展现在自己眼前时的喜悦，又有什么能比得上？感谢老师，感谢此次课设。

附件：
1.1987-1992年连续六年，前三年率定期和后两年检验期的计算径流量与实测径流量绝对误差和相对误差计算程序：
%**************************************************************%
% 1987-1989年计算结果 %
% %
%**************************************************************%
clc
clear
ZL=load('1987到1989资料.txt');%ZL表示导入的数据
WM=140;UM=20;LM=60;DM=60;B=0.3;C=0.16; IM=0.002;
%流域平均张力水容量WM（mm),上层张力水容量UM（mm),下土层张力水容量LM（mm),深层张力水容量DM(mm),张力水蓄水容量曲线方次B，深层蒸散发折算系数C，
%不透水面积占全流域面积的比例IM
Q=ZL(:,3); %日径流量实测值Q（m3/s）
E0=ZL(:,4); %日蒸发量蒸发皿实测值（mm）
P1=ZL(:,5);
P2=ZL(:,6);
P3=ZL(:,7);
P4=ZL(:,8); %P1,P2,P3,P4分别表示流域四个地区的日降雨量实测值（mm）
sumQ=0; %sumQ表示年径流量实测值（mm）
sumR=zeros(4000，1);

h=length(ZL);%求出数据的天数
R=zeros(h,1) ; %R表示日径流量的理论计算值
EU=zeros(1,h);
for i = 1:h %流域平均降雨量计算
P(i) = 0.33 * P1(i) + 0.14 * P2(i) + 0.33 * P3(i) + 0.2 * P4(i); %逐日降雨深的计算
sumQ = sumQ + Q(i) * 24 * 3.6 / 553; %实际测得年径流量的计算
end
%求逐日降雨深和实测全年径流量

W(1) = 110;
WU(1) = 10;
WL(1) = 40;
WD(1) = 60; %初始的土壤总，上，下，深层土壤含水量
%流域三层蒸发计算
WMM = WM * (1 + B); %初始土壤含水量的确定
a(1) = WMM * (1 - (1 - (W(1) / WM)) ^ (1 / (1 + B)));
for j=1:4000
Kc(j)=0.9+0.001*j; %设定Kc值在0.9到4.9的范围内变动，步长取0.001
for i = 1:h
EP(i) = E0(i) * Kc(j); %流域逐日蒸散发深

if WU(i) + P(i) >= EP(i)
EU(i) = EP(i);
EL(i) = 0;
ED(i) = 0;
end
if WU(i) + P(i) < EP(i)
if WL(i) >= C * LM
EU(i) = WU(i) + P(i);
EL(i) = (EP(i) - EU(i)) * WL(i) / LM;
ED(i) = 0;
elseif WL(i) < C * LM & WL(i) >= C * (EP(i) - EU(i))
EU(i) = WU(i) + P(i);
EL(i) = (EP(i) - EU(i)) * C;
ED(i) = 0;
elseif WL(i) < C * (EP(i) - EU(i))
EU(i) = WU(i) + P(i);
EL(i) = WL(i);
ED(i) = (EP(i) - EU(i)) * C - EL(i);
end
end

E(i) = EU(i) + EL(i) + ED(i);%求得各层的蒸发量
PE(i) = P(i) - E(i); %流域产流计算净降雨量PE

%三层蒸散发计算，求总蒸发量和流域净降水量

if PE(i) > 0 %当产流时
if PE(i) + a(i) < WMM
R(i) = PE(i) + W(i) - WM + WM * (1 - (PE(i) + a(i)) / WMM) ^ (B + 1);
W(i + 1) = W(i) + PE(i) - R(i);
a(i + 1) = PE(i) + a(i);
elseif PE(i) + a(i) >= WMM
R(i) = PE(i) + W(i) - WM;
W(i + 1) = WM;
a(i + 1) = WMM;
end
end
%产流计算完毕

if WU(i) + P(i) - EU(i) - R(i) <= UM
WU(i + 1) = WU(i) + P(i) - EU(i) - R(i);
WL(i + 1) = WL(i) - EL(i);
WD(i + 1) = WD(i) - ED(i);
else
WU(i + 1) = UM;
if WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) <= LM
WL(i + 1) = WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM);
WD(i + 1) = WD(i) - ED(i);
else
WL(i + 1) = LM;
if WD(i) - ED(i) + WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) - LM <= DM
WD(i + 1) = WD(i) - ED(i) + WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) - LM;
else
WD(i + 1) = DM;
end
end
end
%此处计算的是降雨补充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，补充来源为降雨减去蒸散发减去径流量，顺序为上、下、深依此补充。

if PE(i) <= 0 %当不产流时
R(i) = 0;
W(i + 1) = W(i) + PE(i);
a(i + 1) = WMM * (1 - (1 - W(i + 1) / WM) ^ (1 / (1 + B)));
end

sumR(j) = sumR(j) + R(i);
end
end

%选出误差最小的年份和差额

cha=zeros(4000,1);
for j=1:4000
kc(j)=0.9+0.001*j;
cha(j)=abs((sumR(j)-sumQ)/sumQ); %cha表示年径流量理论值和实测值的相对误差
end

plot(kc,cha); %计算1989年和1990年相对误差随着Kc值变化图
xlabel('Kc值');ylabel('相对误差值'); %坐标轴表示对象标签

%1987-1989
cha=zeros(4000,1); %直接筛选出1987年误差最小的年份和差额
i=1;
for j=1:4000
cha(j)=abs((sumR(j)-sumQ)/sumQ);
if cha(j)<=0.05
KC(i)=kc(j);
i=i+1 ;
end
end %由此得到KC取1.263-1.430之间数据时，误差均在5%以内

[y,l]=min(cha); %y是误差值，l是误差最小时j的取值
Kc8=0.9+l*0.001; %Kc的最优取值为1.3450

2.1987-1992年六年日模产流量计算程序：
此处仅以1987年为例
%**************************************************************%
% 1987年日径流量理论值和实测值的检验 %
% %
%**