水文预报课程设计报告-江西滁州流域不同参数组合下产流量模  拟效果对比分析

水文预报课程设计报告心得体会1绪论1.1设计任务根据已给的降水、蒸发资料，选择不同的产流参数范围，编写程序，将滁州流城作为整体进行产流量计算；分析不同参数组合对产流量计算结果的影响程度。1.2设计目的1、通过课程设计熟悉降水径流预报的内容，熟悉新安江模型的原理和结构，掌握降水径流预报方案编制的基本步骤；2、掌握参数的物理意义，了解产流参数大小设置对结果的影响。3、加深对水文预报方法的掌握，学会各种方法的综合运用；4、锻炼从事专业工作的基本能力，培养分析问题和解决问题的能力。5、掌握计算机编程的基本方法和基本操作；熟练OFFICE文档操作和研究报告的规范编制。1.3技术要求1、采用泰森单元法进行面平均降水量计算；2、产流采用三江源新安江模型；3、每人计算3种组合情景下的产流量，并进行对比。分析参数改变对计算径流量的影响。1.4工作步骤1、产流量计算；2、划分水源，得出三种水源径流量(径流深)；3、根据参数设置，分析参数改变对径流深的影响。1.5设计要求根据已给资料、参数，自己独立编写程序，将流域作为整体进行产流量计算、水源划分;绘制径流过程；给出参数组合情景下的产流量对比结果，分析参数设置对产流量的影响。设计时间为2周，要求独立完成，遇到问题可以相互讨论，人工计算与计算机计算相结合。所有数据资料共享，成果表达方式要求多样化(包含文字、图、表等)。报告撰写要认真，不得抄袭。可以先人工计算一定时间长度的降水径流过程，以便于计算机编程并作程序调试。1.6有关说明1、程序代码编制中建议变量名：PJ表平均雨量，Ep表示蒸散发能力；其他参数名称设置参考教材。2材料与方法2.1流域基本情况2.1.1自然地理概况遂川江流域主要位于江西省遂川县境内，是赣江一级支流，是吉安市五大河流之一，发源于湖南省桂东县北部的龙潭脑，河源位于东经113º56′，北纬26º11′。遂川江流域以山地和丘陵为主，流域内植被较好，森林覆盖率达70%。流域上游为山区，河床多为砾石、岩石，河面较窄，比降较大，水浅流速大，中下游为低丘区，河床多为砾沙、卵石、细沙。遂川流域属中亚热带湿润季风气候区。遂川江由左、右溪组成，从西南向东北纵贯全线，总长249.5km，流域面积2824km2，流域平均高程511m，流域平均坡度0.31m/km2，流域长度100km，流域形状系数0.29，主河道长度176.0km，河道纵比降2.36‰。2.1.2水文特征流域多年平均降水量为1637mm。降水的年际变化较大，降水主要集中在汛期4～9月，特别其中在4、5、6月，即主汛期，约占全年降水的50%。年降水天数为138～183天，降雨量总的分布是东部和西部大于中部，北部大于南部。多年平均蒸发量为700～1100mm，蒸发量的年变化与气温的年变化趋势大体一致，一般是12～1月份的蒸发量为最小，7～9月份的蒸发量为最大。多年平均径流量1416×108m3，多年平均径流深847mm。径流年内分配是：从1月的2～3%开始逐月上升，5月或6月达到18～25%，从6月或7月开始逐月下降，12月约占2～3%，径流年际变化较大。2.1.3暴雨洪水特征暴雨是江西省遂川县主要灾害之一，以4～8月出现较多。暴雨过程一般持续1～2d，个别年份暴雨带在流域内南北摆动，持续10d以上，例如1962年6月，1998年7月。暴雨强度一般是日降水量为50～100mm。遂川县各河流主汛期大致在4～7月，各站最大洪峰流量和最大峰现时间多发生在5月和6月，干流各站一次洪水过程多表现为峰高量大，历时一般为15～20d，单峰上游出现较多，历时一般为1～3d，由于暴雨带的南北摆动，暴雨范围大、强度大、历时长，往往遂川江的洪水为几条河流的洪水同时组成。2.1.4水文资料遂川江流域面积288.24km2，包括小夏、桐古、营盘墟、沙湖里、阡陌、淋洋、七岭、滁州8个水文站的集水面积。本次课程设计选取滁州水文站（东经114º08′，北纬26º21′）以上集水区间。站点分布情况参见图1和图2。站点权重见表1。模型参数取值范围见表2。图1遂川江流域图图2遂川江流域滁州区间泰森多边形表1遂川江流域滁州区间泰森多边形面积权重站号站名面积（km2）面积权重1小夏44.9115.58%2桐古14.715.10%3营盘墟42.2814.67%4沙湖里37.7213.09%5阡陌43.2114.99%6淋洋53.7418.65%7七岭27.519.54%8滁州24.168.38%合计288.24100%表2模型参数取值范围KCWMWUMWLMWDMBCIMPSMEX0.901302540650.20.120.0013510.951402550650.250.140.002401.21.001202560350.30.160.003451.3初始张力水含量初始SWWUWLWDS110104060352.2新安江流域水文模型2.2.1模型介绍1973年，河海大学赵人俊教授领导的研究组在编制新安江洪水预报方案时，汇集了当时在产汇流理论方面的研究成果，并结合大流域洪水预报的特点，设计了国内第一个完整的流域水文模型——新安江流域水文模型，以下简称新安江模型[]。最初研制的是二水源新安江模型，20世纪80年代中期，借鉴山坡水文学的概念和国内外产汇流理论的研究成果，提出了三水源新安江模型。三水源新安江模型蒸散发计算采用三层模型；产流计算采用蓄满产流模型；用自由水蓄水库结构将总径流划分为地表径流、壤中流和地下径流3种；流域汇流计算采用线性水库；河道汇流采用马斯京根分段连续演算或滞后演算法[]。2.2.2模型结构为了考虑降水和流域下垫面分布不均匀的影响，新安江模型的结构设计为分散性的，分为蒸散发计算、产流计算、分水源计算和汇流计算四个层次结构。每块单元流域的计算流程如图3所示。图3三水源新安江模型流程图图3中方框外为参数，方框内为状态变量。模型各层次结构的功能、计算采用的方法和相应参数见表3。表3新安江模型各层次结构功能、计算采用的方法和相应参数表层次（第一层次）（第二层次）（第三层次）（第四层次）功能蒸散发计算产流计算水源划分汇流计算二水源三水源坡面汇流河道汇流方法三层模型蓄满产流稳定下渗率自由水蓄水库单位线或线性水库或滞后演算法马斯京根或滞后演算法参数KC、UM、LM、CWM、B、IMFCSM、EX、KG、KIUH或CS、CI、CGKE、XE或L2.2.3计算方法2.2.3.1降雨量计算为了考虑降雨分布不均匀和下垫面分布的不均匀性，采用泰森多边形法将遂川江流域滁州区间划分为八个单元流域，在每块单元流域内各有一个雨量站，以每个多边形内雨量站的雨量代表该多边形面积上的降雨量，最后按面积加权推求流域平均降雨量。计算公式如下：P=式中：P为流域平均降雨量；Pi为流域内第i个雨量站同时段降雨量；f为第i个雨量站所在多边形的面积；F2.2.3.2蒸发量计算在新安江模型中，流域蒸散发计算没有考虑流域内土壤含水量在面上分布的不均匀性，而是按土壤垂向分布的不均匀性将土层分为三层，用三层蒸散发模型计算蒸散发量。参数有流域平均张力水容量WM（mm），上层张力水容量UM（mm），下土层张力水容量LM（mm），深层张力水容量DM（mm），蒸散发折算系数KC和深层蒸散发扩散系数C，计算公式如下：WM=W=WUE=EP式中：W为总的张力水蓄量，mm；WU为上层张力水蓄量，mm；WL为下层张力水蓄量，mm；WD为深层张力水蓄量，mm；E为总的蒸散发量，mm；EU为上层蒸散发量，mm；EL为下层蒸散发量，mm；ED为深层蒸散发量，mm；EP为蒸散发能力，mm；蒸散发折算系数三层蒸发模式按照先上层后下层的次序，具有分以下四种情况计算：（1）当WU+EU（2）当WU+EU（3）当WU+EU（4）当WU+EU2.2.3.3产流计算产流计算中采用蓄满产流机制。蓄满是指包气带的土壤含水量达到田间持水量。蓄满产流是指：降水在满足田间持水量以前不产流，所有的降水都被土壤所吸收；降水在满足田间持水量以后，所有的降水（扣除同期蒸发量）都产流。其概念就是设想流域具有一定的蓄水能力，当这种蓄水能力满足以后，全部降水变为径流，产流表现为蓄量控制的特点。湿润地区产流的蓄量控制特点，解决了产流计算在这些地区处理雨强和入渗动态过程的问题;而降雨径流理论关系的建立，解决了考虑流域降雨不均匀的分布式产流计算问题。按照蓄满产流的概念，采用蓄水容量-面积分配曲线来考虑土壤缺水量分布不均匀的问题。所谓蓄水容量—面积分配曲线是指：部分产流面积随蓄水容量而变化的累计频率曲线。应用蓄容量-面积分配曲线可以确定降雨空间分布不均匀情况下蓄满产流的总径流量。对于闭合流域，流域蓄水容量-面积分配曲线采用抛物线形为宜。流域蓄水容量-面积分配曲线与降雨径流相互转换关系如图4所示。(a)(b)图4流域蓄水容量-面积分配曲线与降雨径流间关系图(a)流域蓄水容量-面积分配曲线(b)流域蓄水容量-面积分配曲线与降雨径流关系主要计算公式如下：WM=a=若PE+R=若PE+R=2.2.3.4分水源计算按蓄满产流模型计算出的总径流量R中包括了各种径流成分，由于各种水源的汇流规律和汇流速度不同，相应采用的计算方法也不同。因此，必须进行水源划分。本次使用的新安江模型采用三水源划分。三水源的水源划分结构借鉴了山坡水文学的概念，去掉了FC，用自由水蓄水库结构解决水源划分问题。自由水蓄水库结构如图5所示。自由水蓄水库结构考虑了包气带的垂向调蓄作用。按蓄满产流模型计算出总径流量R，先进入自由水蓄水库调蓄，再划分水源。从图可见，产流面积上自由水蓄水库设置了两个出口，一个为旁侧出口，形成壤中流RI；另一个为向下出口，形成地下径流RG。根据蓄满产流的概念，只有在产流面积上才可能产生径流，因为产流面积是变化图5自由水蓄水库结构图的，所以，自由水蓄水库的底宽也是变化的。在图中还设置了一个壤中流水库，该水库用于壤中流受调蓄作用大的流域，即将划分出来的壤中流再进行一次调蓄计算。该水库一般是不需要的，故在图中用虚线表示。由于饱和坡面流的产流面积是不断变化的，所以在产流面积FR上自由水蓄水容量分布是不均匀的。三水源水源划分结构是采用类似于流域蓄水容量—面积分配曲线的流域自由水蓄水容量—面积分配曲线来考虑流域内自由水蓄水容量分布不均匀的问题。主要计算公式如下：流域平均自由水蓄水容量为：Sm时段初自由水蓄量S1S1相应的纵坐标AU为：AU=当PE+AU<RS=当PE+AU≥RS=时段末的自由水蓄量S2S2相应的壤中流RI为：RI=地下径流RG为：RG=3结果与讨论3.1产流计算根据VBA编程计算得出的数据，应用oringe软件绘制降雨及径流的相关图形如下。(a)81-83年径流深过程图(b)径流深与降雨对比图图6组合一(a)81-83年径流深过程图(b)径流深与降雨对比图图7组合二(a)81-83年径流深过程图(b)径流深与降雨对比图图8组合三图981-83年三种组合径流深过程对比图由图6-9可知，三种参数组合情况下的径流过程基本一致，变化较小，且径流过程与降雨过程有明显的对应关系，径流受降雨影响较大。3.2水源划分利用VBA编程将径流划分为地表径流、壤中流与地下径流三种成分，应用oringe软件绘制的三种水源过程图如图10-12所示。图1081-83年三水源过程图(组合一)图1181-83年三水源过程图(组合二)图1281-83年三水源过程图(组合三)由图10-12可知，三种参数组合情况下的地表径流、壤中流及地下径流过程基本一致。每种组合的三种水源成分占径流的比例也基本相同。其中，地表径流在集中降雨大时出现较高值，原因是流域已经蓄满从而产生大量地表径流。壤中流与地下径流受蓄水容量限制变化较为平缓，几乎相等。3.3模型参数新安江模型参数可分蒸散发计算、产流计算、分水源计算和汇流计算四个层次，本次课设只进行到第三层次，相应的参数如表4所示。表4新安江模型各层次参数表层次参数符号参数意义敏感程度参数值第一层次蒸散发计算KC流域蒸散发折算系数敏感0.90/0.95/1.00UM上层张力水容量（mm）敏感20/25/30LM下层张力水容量（mm）敏感40/50/60C深层蒸散发折算系数不敏感0.12/0.14/0.16第二层次产流计算WM流域平均张力水容量（mm）不敏感120/130/140B张力水蓄水容量曲线方次不敏感0.20/0.25/0.30IM不透水面积占全流域面积的比例不敏感第三层次水源划分SM表层自由水蓄水容量（mm）敏感35/40/45EX表层自由水蓄水容量曲线方次不敏感1.0/1.2/1.3KG表层自由蓄水水库对地下水的日出流系数敏感0.35KI表层自由蓄水水库对壤中流的日出流系数敏感0.35分别进行三个层次的参数敏感性分析：1、蒸散发计算表5参数KC的敏感性分析KC0.900.951.00年径流R/mm4432.74340.64250.1WM=130，WUM=25，WLM=40，WDM=65，WMM=156，B=0.2，C=0.12KC是影响产流量计算最为重要和敏感的参数，产流计算中KC控制着水量平衡，因此，对水量计算是最重要的。它主要反映流域平均高程与蒸发站高程之间差别的影响和蒸发皿蒸散发与陆地蒸散发间差别的影响。该参数值的改变，影响初始土壤含水量和雨期蒸发，导致计算产流量的改变。如表5所示，KC值增大，蒸发量增大，产流量减少；KC减小，蒸发量降低，产流量增大。表6参数WUM的敏感性分析WUM202530年径流R/mm4341.84334.74321.6KC=0.95，WM=140，WLM=50，B=0.25，C=0.14WUM与作物根系土层厚度有关，其取值范围为10-50mm。由表6可知，产流量R随着上层张力水容量WUM值的增加而减小，且变化幅度明显，较为敏感。表7参数WLM的敏感性分析WLM405060年径流R/mm4432.74433.04443.2KC=0.90，WM=130，WUM=25，B=0.2，C=0.12WLM与包气带的土层结构和物理特性有关，其值约在60-90mm（湿润地区）之间。产流量R随着下层张力水容量WLM值的增加而变大，其变化幅度较WUM更明显，敏感性更强。表8参数C的敏感性分析C0.120.140.16年径流R/mm4334.74334.74334.7KC=0.95，WM=140，WUM=25，WLM=50，B=0.25，C主要取决于流域内深根植物的覆盖面积，是反映深根植物作用的参数，决定深层散发。由表8可知，该值的变化对年径流的影响很小，参数不敏感。2、产流计算表9参数WM的敏感性分析WM120130140年径流R/mm4367.14341.64340.9KC=0.95，WUM=30，WLM=50，B=0.25，C=0.14流域平均张力水容量WM表示流域干旱程度和影响土壤水分变化的土层深度。WM在模型中不是很灵敏，如WM分别取值130和140时，年径流变化不大，但WM值不能太小，可能会出现负值。表10参数B的敏感性分析B0.200.250.30年径流R/mm4432.74435.44437.9KC=0.90，WM=130，WUM=25，WLM=40，C=0.12B值反映划分单元流域张力蓄水分布的不均匀程度。B值越大表示流域越不均匀。其值的改变对年径流的影响较小，参数敏感性不强。3、水源划分表11参数SM的敏感性分析SM354045RS1111.3997.0902.8RI1520.51573.61617.4RG1520.51573.61617.4WM=130，EX=1.0，B=0.2，KI=0.35，KG=0.35SM反映表土蓄水能力，其值受降雨资料时段均化的影响明显。当以日为计算时段长时，一般流域的SM值约为10～50mm，当所选取的计算时段长较小时，SM要增大，这个参数对地面径流的多少起着决定性作用。如表11所示，SM增大，RS减小，RI与RG均增大，三者变化显著。表12参数EX的敏感性分析EX1.01.21.3RS1111.31151.61168.7RI1520.51502.01494.2RG1520.51502.01494.2WM=130，SM=35，B=0.2，KI=0.35，KG=0.35EX反映流域自由水蓄水分布的不均匀程度，该值范围一般为1.0～1.5，由表12可知，随着EX的增大，RS逐渐增大，RI与RG均减小。3.4误差分析影响流域降雨径流过程的因素很多，三水源新安江模型的结构与参数能够反映湿润地区降雨径流过程的主要规律与特点，因而能获得较好的精度。但是，模型本身以及模型计算中有许多概化，会造成误差。造成遂川江流域年径流量模拟误差来源主要有以下几方面：1、雨量站代表性的影响遂川江流域面积288.24km2，包括小夏、桐古、营盘墟、沙湖里、阡陌、淋洋、七岭、滁州8个水文站的集水面积。对于多年平均降雨来讲，基本上能控制降雨的空间分布。但对于不同年份不同时期，不同类型洪水，其差别较大。对于次洪水来说，在不同的气象条件下计算出的面平均雨量也会有较大的误差，面平均雨量的计量误差是产流量误差的主要影响因素。因此，在计算流域径流过程中，要充分考虑流域的气象条件和地质的时空不均匀性。2、人类活动的影响随着社会经济的快速发展，人类活动的影响加剧，流域内先后修建了一些中小型水库，这些中小型水库大多数没有固定调度原则和调洪方式，干旱季节或初汛的时候，流域内的大小水利工程均需蓄水，会使实测的径流量偏小。人类活动的影响也会造成流域下垫面的变化，下垫面的改变与不均匀分布将会产生一定的误差。3、模型参数的影响模型参数是根据输入，通过模型计算输出，再将输出过程与实测过程进行比较，用系统识别的方法作优化调试的，因为需要率定的参数较多，有一部分参数是直接给定的，因此上述所率定的模型参数可能不是最优。此外，用来率定参数的资料序列太短，资料的代表性不强，并且概念性水文模型本身精度也有限。3.5问题讨论3.5.1流域划分根据遂川江滁州水文站（东经114º08′，北纬26º21′）以上控制流域水系及水文站雨量站网的分布特点，天然流域划分为小夏、营盘墟、桐古、沙湖里、淋洋、阡陌、七岭、滁州8块单元面积，然后再在每一块面积内用泰森多边形细分。这样可以更符合降雨空间分布实际情况，解决面雨量计算误差问题。3.5.2实时校正模型计算值与实测值（流量或水位）之间总是存在一定的误差。造成两者间误差的因素很多，若针对每一个单一因素，它们是难于描述和预见的，一般采用实时校正模型来解决。实时校正模型的种类有很多，主要与所选择的预报模型有关。可以充分利用流域内水文（或水位）站的信息，对模型计算值进行实时校正。4参考文献[1]刘新仁.流域水文模型的研究途径[J].全国水文学术讨论会论文集，2004，189-192.[2]赵人俊.流域水文模型-新安江模型与陕北模型[M].北京：水利电力出版社，1984.[3]赵人俊.新安江模型参数的分析[J].水文，1998（6）：2-9.[4]包为民.水文预报（第4版）[M].北京.中国水利水电出版社，2016.[5]《水文预报课程设计任务书指示书》.扬州大学水利与能源动力工程学院.5附录：程序代码5.1蓄满产流计算Subtest1()DimPJAsDouble'流域日平均降雨量，mmDimE0AsDouble'流域水面蒸发，mmDimEPAsDouble'流域蒸发能力，mmDimEUAsDouble'土壤上层蒸发量，mmDimELAsDouble'土壤的下层蒸发量，mmDimEDAsDouble'土壤的深层蒸发量，mmDimEAsDouble'土壤总蒸发量，mmDimPEAsDouble'扣除雨期蒸发后的降雨量，mmDimWUAsDouble'上层土壤含水量，mmDimWLAsDouble'下层土壤含水量，mmDimWDAsDouble'深层土壤含水量，mmDimWAsDouble'总含水量，mmDimaAsDouble'流域初始含水量最大值，mmDimRAsDouble'流域蓄满产流量，mmDimKCAsDouble'蒸发折算系数DimWMAsDouble'流域平均蓄水容量，mmDimWUMAsDouble'上层土壤含水容量，mmDimWLMAsDouble'下层土壤含水容量，mmDimWDMAsDouble'深层土壤含水容量，mmDimWMMAsDouble'流域蓄水容量最大值，mmDimBAsDouble'蓄水容量曲线指数，mmDimCAsDouble'蒸发扩散系数DimxAsIntegerForx=3To1097Step1PJ=Range("b"&x)E0=Range("c"&x)EP=Range("d"&x)EU=Range("e"&x)EL=Range("f"&x)ED=Range("g"&x)E=Range("h"&x)PE=Range("i"&x)WU=Range("j"&x)WL=Range("k"&x)WD=Range("l"&x)W=Range("m"&x)a=Range("n"&x)R=Range("o"&x)KC=0.9WM=130WUM=25WLM=40WDM=65WMM=156B=0.2C=0.12'三层蒸发模式计算IfRange("j"&x)+Range("b"&x)>=Range("d"&x)Then'第一种情况Range("e"&x)=Range("d"&x)Range("f"&x)=0Range("g"&x)=0ElseIfRange("j"&x)+Range("b"&x)<Range("d"&x)AndRange("k"&x)>=C*WLMThen'第二种情况Range("e"&x)=Range("j"&x)+Range("b"&x)Range("f"&x)=(Range("d"&x)-Range("e"&x))*Range("k"&x)/WLMRange("g"&x)=0ElseIfRange("j"&x)+Range("b"&x)<Range("d"&x)AndC*(Range("d"&x)-Range("e"&x))<=Range("k"&x)AndRange("k"&x)<C*WLMThen'第三种情况Range("e"&x)=Range("j"&x)+Range("b"&x)Range("f"&x)=C*(Range("d"&x)-Range("e"&x))Range("g"&x)=0ElseIfRange("j"&x)+Range("b"&x)<Range("d"&x)AndRange("k"&x)<C*(Range("d"&x)-Range("e"&x))Then'第四种情况Range("e"&x)=Range("j"&x)+Range("b"&x)Range("f"&x)=Range("k"&x)Range("g"&x)=C*(Range("d"&x)-Range("e"&x))-Range("f"&x)EndIfRange("h"&x)=Range("e"&x)+Range("f"&x)+Range("g"&x)'土壤总蒸发量计算Range("i"&x)=Range("b"&x)-Range("h"&x)'扣除雨期蒸发后的降雨量计算Range("m"&x)=Range("j"&x)+Range("k"&x)+Range("l"&x)'总含水量计算Range("n"&x)=WMM*(1-(1-Range("m"&x)/WM)^(1/(1+B)))'流域初始含水量最大值计算'流域蓄满产流计算IfRange("i"&x)>0Then'产流SelectCaseRange("n"&x)+Range("i"&x).ValueCaseIs<=WMM'全流域蓄满前产流量计算Range("o"&x)=Range("i"&x)+Range("m"&x)-WM+WM*(1-(Range("i"&x)+Range("n"&x))/WMM)^(B+1)CaseElse'全流域蓄满后产流量计算Range("o"&x)=Range("i"&x)-(WM-Range("m"&x))EndSelectIfRange("j"&x)+Range("i"&x)-Range("o"&x)<=WUMThen'上层土壤未蓄满时Range("j"&x+1)=Range("j"&x)+Range("i"&x)-Range("o"&x)Range("k"&x+1)=Range("k"&x)Range("l"&x+1)=Range("l"&x)ElseIfRange("j"&x)+Range("i"&x)-Range("o"&x)>WUMAndRange("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WUM<=WLMThen'上层土壤蓄满、下层土壤未蓄满时Range("j"&x+1)=WUMRange("k"&x+1)=Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WUMRange("l"&x+1)=Range("l"&x)ElseIfRange("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WUM>=WLMAndRange("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WDM<WDMThen'下层土壤蓄满、深层土壤未蓄满时Range("j"&x+1)=WUMRange("k"&x+1)=WLMRange("l"&x+1)=Range("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WDMElseIfRange("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)-Range("o"&x)-WDM>=WDMThen'深层土壤蓄满时Range("j"&x+1)=WUMRange("k"&x+1)=WLMRange("l"&x+1)=WDMEndIfElseIfRange("i"&x)<=0Then'不产流IfRange("j"&x)+Range("i"&x)>0Then'上层土壤含水量未蒸发完时Range("j"&x+1)=Range("j"&x)+Range("i"&x)Range("k"&x+1)=Range("k"&x)Range("l"&x+1)=Range("l"&x)ElseIfRange("k"&x)+Range("j"&x)+Range("i"&x)>0Then'上层土壤含水量蒸发完、下层未蒸发完时Range("j"&x+1)=0Range("k"&x+1)=Range("k"&x)+Range("j"&x)+Range("i"&x)Range("l"&x+1)=Range("l"&x)ElseIfRange("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)>0Then'下层土壤含水量蒸发完、深层未蒸发完时Range("j"&x+1)=0Range("k"&x+1)=0Range("l"&x+1)=Range("l"&x)+Range("k"&x)+Range("j"&x)+Range("i"&x)ElseIfRange("l"&x)<=0Then'深层土壤含水量蒸发完时Range("j"&x+1)=0Range("k"&x+1)=0Range("l"&x+1)=0EndIfEndIfNextxEndSub5.2水源划分Subtest2()DimPEAsDouble'扣除雨期蒸发后的降雨量，mmDimRAsDouble'流域蓄满产流量，mmDimFRAsDouble'产流面积DimAUAsDouble'时段初始自由水蓄量相应的纵坐标DimSAsDouble'自由水蓄量，mmDimRSAsDouble'地面径流，mmDimRIAsDouble'壤中流，mmDimRGAsDouble'地下径流，mmDimWAsLong'总含水量，mmDimWMAsDouble'流域平均蓄水容量，mmDimSMAsDouble'深度，mmDimEXAsDouble'自由水蓄量分布曲线指数DimBAsDouble'产流面积计算指数DimSmmAsDouble'流域最大自由蓄水量，mmDimKIAsDouble'壤中流的出流系数DimKGAsDouble'地下径流的出流系数DimxAsIntegerForx=4To1098PE=Range("b"&x)R=Range("c"&x)FR=Range("d"&x)AU=Range("e"&x)S=Range("f"&x)RS=Range("g"&x)RI=Range("h"&x)RG=Range("i"&x)W=Range("j"&x)WM=130SM=35EX=1B=0.2Smm=70KI=0.35KG=0.35'三水源划分SelectCaseRange("b"&x).ValueCaseIs<=0Range("d"&x)=1-(1-Range("j"&x)/WM)^(B/(1+B))Range("g"&x)=0Range("h"&x)=Range("f"&x)*KI*Range("d"&x)Range("i"&x)=Range("f"&x)*KG*Range("d"&x)Range("f"&x+1)=Range("f"&x)*(1-KI-KG)CaseElseRange("d"&x)=Range("c"&x)/Range("b"&x)Range("e"&x)=Smm*(1-(1-Range("f"&x)*Range("d"&x-1)/Range("d"&x)/SM)^(1/(1+EX)))SelectCaseRange("b"&x).Value+Range("e"&x).ValueCaseIs<SmmRange("g"&x)=Range("d"&x)*(Range("b"&x)+(Range("f"&x)*Range("d"&x-1)/Range("d"&x))-SM+SM*(1-((Range("b"&x)+Range("e"&x))/Smm)^(1+EX)))Range("h"&x)=(SM-SM*(1-(Range("b"&x)+Range("e"&x))/Smm)^(1+EX))*KI*Range("d"&x)Range("i"&x)=(SM-SM*(1-(Range("b"&x)+Range("e"&x))/Smm)^(1+EX))*KG*Range("d"&x)Range("f"&x+1)=(1-KI-KG)*(SM-SM*