新安江流域水文模型

1、第二章 新安江流域水文模型60年代初，河海大学(原华东水利学院)水文系赵人授等开始研究蓄满产流模型，配合一定的汇流计算，将模型应用于水文预报和水文设计。1973年，他们在对新安江水库做人库流量预报的工作中，把他们的经验归纳成一个完整的降雨径流流域模型新安江模型。模型可用于湿润地区和半湿润地区的湿润季节径流模拟和计算。最初的新安江模型为两水源模型，只能模拟地表径流和地下径流。80年代初期，模型研制者将萨克拉门托模型与水箱模型中，用线性水库函数划分水源的概念引入新安江模型，提出了三水源新安江模型，模型可以模拟地面径流、壤中流、地下径流。1984至1986年，又提出了四水源新安江模型，可以模拟地面径

2、流、壤中流、快速地下径流和慢速地下径流。三水源新安江模型一般应用效果较好，但模拟地下水丰富地区的日径流过程精度不够理想。在新安江三模型中增加慢速地下水结构就成为四水源新安江模型。 当流域面积较小时，新安江模型采用集总模型，当面积较大时，采用分块模型。分块模型把流域分成许多块单元流域，对每个单元流域做产、汇计算，得到单元流域的出口流量过程。再进行出口以下的河道洪水演算，求得流域出口的流量过程。把每个单元流域的出流过程相加，就求得了流域出口的总出流过程。划分单元流域的主要目的是处理降雨分布的不均匀性，因此单元流域应当大小适当，使得每块面积上的降雨分布比较均匀并有一定数目的雨量站。其次尽可能使单元流

3、域与自然流域相一致，以便于分析与处理问题，并便于利用已有的小流域水文资料。如果流域内有大中型水库，则水库以上的集水面积即应作为一个单元流域。因为各单元流域的产汇、流计算方法基本相同，以下只讨论一个单元流域的情况。2.1新安江两水源模型 1模型结构和参数 新安江两水源模型的产流子模型采用蓄满产流模型，蒸发计算采用三层蒸发计算模型。利用稳定下渗率将径流划分为地面径流和地下径流两种水源。地面径流采用单位线汇流，地下径流采用一次线性水库汇流。模型把流域面积划分为透水面积和不透水面积两部分，不透水面积上的降水在满足蒸发后将直接转化为地面径流。透水面积上将发生下渗，下渗的水量一部分存储于土壤层，后期耗于蒸

4、发；满足了流域土壤蓄水容量后的下渗水量才能转化为径流。不透水面积用参数表示，它是用流域内不透水面积占全流域面积的百分比表示的。新安江模型的输出是流域出流过程和流域蒸散发过程，输入则为时段降雨量、蒸发皿观测蒸发量。 新安江两水源模型共有9个参数，一条单位线。 流域蒸发折算系数，是流域蒸散发能力与蒸发皿蒸发量之比； C深层蒸散发系数； 不透水面积占全流域面积的百分比，； 流域平均蓄水容量(指张力水)，； 流域平均上层蓄水容量, ； 流域平均下层蓄水容量, ；蓄水容量曲线指数；稳定下渗率；地下水消退系数；单位线。模型结构如下图所示：图0 新安江两水源模型结构示意图图中方框内标注为状态变量，方框外标注

5、的是模型参数。2模型参数的调试当模型初定后就可上机调试，在调试的过程中，应注意各参数的径流响应特征，以及参数之间相依性的影响，以便调试能有目的地顺利进行。以下简介的调试技术可供参考。流域蒸发折算系数的调试，对一个具体流域来说，此参数完全靠优选。调试此参数主要考察年径流模拟误差，当年径流模拟误差达最小时，此值最优，但个别参数的最优并不表示模型总体最优。根据现有经验，年径流模拟误差控制在58左右就可以了。流域平均蓄水容量的调试，一般根据实测的降雨径流资料分析或经验选取后，适当微调即可。经验选定后，此值调试范围不大，比较容易确定。通常在南方湿润地区可选80150左右，江淮一带约在110左右,燕山东北

6、东部地区，选150左右。流域平均上层蓄水容量，一般取1020，流域平均下层蓄水容量一般取6090。蓄水容量曲线指数的调试，表示流域蓄水容量分布的不均匀性，当全流域蓄满后，这个参数就不起作用了，因此，应该选取流域没有达到蓄满的那些洪水点据作为调试的依据。此值一般在0.20.5之间取值。注意值与值之间有相依性，二者对径流模拟的结果相互有影响。稳定下渗率的调试，对水源划分起决定性作用，但至今对的研究还不充分，目前可对作一些简单处理，如果取为一个常数使模拟效果不好时，可考虑作为变动参数。考察是否合适主要观察地下径流的模拟精度。地下水消退系数的调试，主要观察洪水退水段的拟合精度。不透水面积参数对模型拟合

7、精度影响很小，一般取00.05。2.2新安江三水源模型一、新安江三水源模型结构 新安江三水源模型包括4个计算环节：流域产流计算；径流的划分；蒸散发计算；汇流计算。流域的产流计算和蒸散发计算与新安江二水源模型相同，水源划分则完全不同，因此汇流计算也不相同。新安江三水源模型结构示意图如下 二、新安江三水源模型流域产流计算 新安江模型产流部分的计算是蓄满产流模式， 蓄满产流指在流域包气带土湿满足田间持水量以前不产流，所有的降雨都被土壤吸收；而在土湿达到田间持水量之后，所有的降雨（除去同期的蒸散发）都产流。在产流后，流域包气带土壤的下渗能力为稳定下渗率，下渗的水分成为地下径流和壤中流，超蓄的部分成为地

8、面径流。考虑到流域内各点的蓄水容量并不相同，实际产流时常常是在部分面积上产流，新安江模型引入流域蓄水容量曲线来刻划流域内各点蓄水容量的不均匀性，把流域内各点的蓄水容量概化成如图2.3所示的一条抛物线（也可概化成其它函数形式），其方程为： （2-1）式中： 为产流面积；是点蓄水容量,；是蓄水容量小于等于蓄水量的全部点的面积,；是单元流域面积, ；是蓄水容量曲线指数。是流域内点最大蓄水容量,；其与流域平均蓄水容量的关系为: （2-2）引入蓄水容量曲线后，在降雨过程中，只有在满足了包气带蓄水容量的面积上才可能产生径流，其余面积上不产生径流。扣除蒸发后的有效降雨量中，未转化为净雨产流量的水分通过下渗进

9、入流域包气带土层，补充土层蓄水量。引入流域蓄水容量曲线后的产流计算式为：当时为部分面积产流 （2-3）当时为全流域面积产流 （2-4） 式中： 为有效降雨量， ；为与流域土壤初始蓄水量相应的前期影响雨量，由下式计算 （2-5）三、新安江三水源模型水源划分1不考虑自由水蓄水容量分布不均匀的划分方法求得的产流量包括地面径流 、壤中流和地下径流 三部分。新安江三水源模型用一个自由水蓄水库解决水源划分，自由水蓄水库有两个出流孔，底孔为地下径流出流孔，边孔为壤中流出流孔。新安江模型考虑了产流面积（）的变化，自由水蓄水库实际只发生在产流面积上，其底宽为产流面积，显然它是随时间变化的。产流量进入水库即在产流

10、面积上产生的径流深，也就是自由水蓄水库所增加的蓄水深，当自由水蓄水深超过其最大值时，超过部分成为地面径流。壤中流和地下径流按线性水库出流，其出流系数分别为和。底孔出流量和边孔出流量分别进入各自的水库，并按线性水库的退水规律流出，分别成为地下水出流和壤中流出流。模型认为蒸散发在张力水中消耗，自由水蓄水库的水量全部为径流。水源划分示意图见图3。水源划分计算式为：地面径流计算式：当 （2-6）当 （2-7）壤中流计算式： （2-8）地下径流计算式： （2-9）注意每步计算时应将自由水蓄水库的蓄量折算为当前时段产流面积上的深度。2.考虑自由水蓄水容量分布不均匀的划分方法流域自由水蓄水容量实际上也是不均

11、匀的，作为图3表示的水源划分方法的改进，人们引入流域自由水蓄水容量分布曲线来刻化自由水蓄水量的不均匀性，其线型与流域蓄水容量曲线类似，曲线见图4，曲线方程如下式。 （2-10）式中： 为产流面积；是点自由水蓄水容量；是自由水蓄水容量小于等于蓄水量的全部点的面积；是单元流域面积；是自由水蓄水容量曲线指数。是流域内点最大自由水蓄水容量；与流域平均自由水蓄水容量的关系为： （2-11）由于认为在产流面积上才有自由水，因此，产流面积上的点最大自由水蓄水容量是随产流面积变化而变化的，见图4。对当前时段，已知，则将其代人（2-10）得 （2-12）1.0从（2-12）解出得： （2-13）同样，在产流面积

12、上的自由水平均蓄水容量与该面积上的点最大自由水蓄水容量的关系为： （2-14）在划分水源时，产流面积上的只有在满足了蓄水容量的水库面积上才可能产生地面径流；其余面积上不产生地面径流。进入水库的水分扣除地面径流量后，剩余的补充水库蓄水量S，这时的水源划分计算式为：地面径流出流当时为部分面积产流 （2-15）当时为全流域面积产流 （2-16）式中：为计算时段初的自由水蓄水库的蓄水量（用产流面积上的平均深度表示）；为计算相应于的自由水蓄水库前期影响入流量，由下式计算； （2-17）进入水库的水分扣除地面径流量后，剩余的补充水库蓄水量S： （2-18）自由水蓄水库的时段入流量为时段产流量，但由（2-3

13、）或（2-4）计算的产流量是用流域平均水深表示的水量，应将折算为产流面积上表示的水深。可以证明折算为产流面积上表示的水深值恰为流域有效降雨值，证明如下：为折算为产流面积上表示的水深，根据降雨径流关系的性质，水文学原理中已经证明相对产流面积与产流量之间有关系： （2-19）从此式可看出壤中流出流： （2-20）地下径流出流： （2-21）由于，都是产流面积上的深度表示的水量，计算所得、应折算回全流域深度表示的水量，即：、在当前时段计算结束后应乘。蒸散发计算四、流域蒸散发计算蒸散发模型不考虑面上分布的不均匀性，但可考虑土湿垂向分布的不均匀性，根据需要采用两层或三层蒸散发模型计算流域蒸发量。两层蒸散

14、发模型将土层分为上、下两层，各层蓄水容量分别为、（）。流域土层下渗蓄水和蒸散发计算过程按下述原则进行：降雨先补充上层，上层蓄满后再补充下层；蒸发时先蒸发上层，上层蓄水量蒸发殆尽后再蒸发下层。1.两层蒸散发模型的计算式如下：， ： ，式中：是流域降雨量，；是流域蒸散发能力，；是上层土壤蓄水量，；是下层土壤蓄水量，； 是下层土壤蓄水容量，；、分别为上层、下层蒸发量，。2.三层蒸散发模型的计算式如下：当：， ，时：，时：，时：，式中：是深层蒸散发系数； 五、流域汇流计算汇流计算可以采用两种方法：1.传统流域汇流这种汇流方法是：地面径流采用单位线汇流，壤中流和地下径流采用线性水库汇流至出口断面后与单位

15、线汇流结果叠加，最后得到流域出口断面的模拟出流。即图3中下面的两个水库的出流就是壤中流和地下径流在流域出口断面的出流流量，该两个水库分别称为壤中流调蓄水库和地下径流调蓄水库。壤中流汇流计算地下径流汇流计算式中： t时段壤中流出流流量，；壤中流消退系数，可假定为；t时段地下径流出流流量，；地下径流消退系数，可假定为；是径流深转换为流量时的单位转换系数，； 上式中假定了调蓄水库与自由水水库的出流系数相同。在实践中，为延缓壤中流和地下径流出流，有人把下面的两个水库的出流看作是坡地调蓄水库，而在它们后面再分别串联一个水库反映壤中流和地下径流受河槽调蓄的作用。2.河网汇流单位线汇流河网汇流单位线也称时变

16、河网汇流单位线，这种汇流方法忽略地面径流的坡面汇流阶段，直接将地面径流产流量转换为流量作为河网地面径流入流（看作集中在计算时段末入流）。地下径流和壤中流汇流均采用线性水库作坡地调蓄后得到相应河网入流，即：把图3中的转换为流量作为河网地面径流入流，下面的两个水库的出流看作壤中流和地下径流经坡地调蓄后注入河网的入流量。三种水源的河网入流流量相加就是所谓的河网总入流。河网总入流经流域河网调蓄后得到流域出口断面的流量过程。河网汇流单位线的计算式如下： （2-22）式中：t时段河网的总入流，；河网消退系数；河网汇流时间，；t时刻的河网出流（即t时刻的流域出口断面流量），；t时刻的河网总入流经河网调蓄后，

17、延迟时间的出流，。是随时间而变化的时变参数，它取决于t时刻的河网的总入流量和流域的河槽特性。一般用下式计算： （2-23）式中：反映流域河槽特性的系数。2.3线性水库不同单位时段的出流系数转换 从(2-20)和(2-21)知道新安江模型的自由水水库出流是按线性水库出流处理的。在客观世界中，流域的降雨、蒸散发、径流等水文物理量都是连续量，流域内的蓄水量也是连续量，观测和计算时都只能处理为一个一个的时段平均量；从数学观点来看，这实际上是将连续量处理为离散量，将表达物理过程的连续方程处理为表达离散量的差分方程，这样就必然带来差分误差，特别是水文模型计算通常采用的时段较大，如果在模型的算法设计上不注意

18、，差分误差可能较大。为了减小计算的差分误差，本可以缩小计算时段，但是，在水文上常常受到限制，这是由于原始水文观测数据时段本身较大，如果在所有的计算环节上都缩小计算时段并不能提高成果的精度。尽管如此，在模型的算法上，至少应要求算法不会带来更大的差分误差。仔细考察新安江流域水文模型计算的各个环节，就会发现(2-20)和(2-21)的自由水水库出流计算可能会出现较大的差分误差，为此，有的模型程序设计者在算法上作了一定处理，以设法减小差分误差。基本的想法是在这个计算环节上减小计算步长，即将模型总体的计算时段在这个环节上进一步细分；为此，需要将水库长时段的出流系数转换为短时段的出流系数。设水库在一个计算

19、时段内的入流量为，为减小差分误差，将等分为份入流，这实际上是将原计算时段细分为个更小的计算时段，每个时段的水库入流为。可以证明水库计算时段为的出流系数和与计算时段为的出流系数和之间存如下关系： ； （2-24）证明如下：1. 只有一个出流孔的线性水库不同时段出流系数的转换 先对只有一个出流孔的线性水库进行证明。设线性水库的出流方程和退水时的水量平衡方程如下： （2-25） （2-25）式中： 为线性水库的瞬时出流系数；为线性水库时刻的蓄量；为线性水库时刻的出流流量。合并（2-25）和（2-26）得到 （2-27）在上积分上式得 （2-28）此式称为线性水库退水方程，称为线性水库消退系数（又称退

20、水系数）。此式的另一个表达形式是： （2-29）这是（2-27）在区间上积分的结果。令，即取为一个单位计算时段，则上式记为 （2-30）在区间上积分得到水库在时段内流出的水量 （2-31）令，即取为一个单位计算时段，则上式记为 （2-32）（2-32）中的就是线性水库单位时段为的出流系数。比较（2-25）和（2-32）可知和的意义是完全不同的，总是与一个单位时段相联系的，而且单位时段不同，其值也不同。令，记单位时段为的水库消退系数为。由（2-30）的退水方程，可推导得 （2-33）此式中的1代表一个时段，表示个时段。显然。对单位时段的退水方程，有； 或 （2-34）比较（2-33）最后一个式子

21、与（2-34）式可知，线性水库单位计算时段为的出流系数与单位计算时段为的出流系数之间有如下关系： （2-35）已知对单位计算时段，有；记单位计算时段的出流系数为，则可推得 （2-36）此式就是所需要得结果。2有二个出流孔的线性水库不同时段出流系数的转换记单位计算时段为的出流系数分别为和，单位计算时段为的出流系数分别为和，不同计算时段的两孔总出流系数分别为和，代人（2-36）式得到： （2-37）要求解和还须增加一个条件，可假定不同计算时段的两孔出流系数分配比例相同，即成立： （2-38）利用（2-38）消去得到： （2-39）利用（2-38）消去得到： （2-40）利用（2-39）与（2-40

22、）就可进行线性水库不同时段的出流系数转换。2.4新安江三水源模型的参数调试 采用河网汇流单位线的新安江三水源模型共有14个参数，根据参数的功能可分为4类 与蒸散发计算有关的参数：蒸散发折算系数，上层土壤蓄水容量下层土壤蓄水容量，深层蒸散发系数。与产流计算有关的参数：流域蓄水容量，流域蓄水容量曲线指数，流域不透水面积比值。与水源划分有关的参数：流域自由水蓄水容量，流域自由水蓄水容量曲线指数，自由水蓄水库壤中流出流系数与自由水蓄水库地下径流出流系数。与汇流有关的参数：壤中流消退系数和地下径流消退系数；河网汇流单位线参数和。壤中流汇流调蓄水库和地下径流汇流调蓄水库与自由水水库是不同的水库，但如果假定

23、二调蓄水库的出流系数与自由水水库对应的径流的出流系数相同，即假定，则参数可减为12个。4类参数相对独立，同类参数有一定相依性。1参数初定根据流域自然地理资料、水文气象资料和前人以积累的经验初定参数。蒸散发折算系数：可采用汛期连续大雨，流域蓄水量处于蓄满状态的累计降雨、累计径流、累计蒸发观测资料，用下式估计: （2-41）流域蓄水容量、上层土壤蓄水容量、下层土壤蓄水容量：值一般根据实测的降雨径流资料分析或经验选取后，适当微调即可。经验选定后，此值调试范围不大，比较容易确定。通常在南方湿润地区可选80150左右，北方半湿润地区可达170，江淮一带约在110左右,燕山东北东部地区，选150左右。 中

24、包含有植物截留，少林地510，多林地1020左右。在6090范围。 流域蓄水容量曲线指数：此值一般在0.20.5之间取值，通常流域越小，该值越小，特小流域可能小于0.2。1000平方公里以上流域一般大于0.3。 不透水面积比值：值对模型拟合精度影响很小，一般取00.05。深层蒸散发系数： 南方多林地区左右 ；北方半干旱地区左右。 流域自由水蓄水容量： 一般由经验初定，土层浅薄的山区，土厚多林地区，甚至更大。 流域自由水蓄水容量曲线指数在； 根据经验，常在11.5左右。 自由水蓄水库壤中流出流系数与自由水蓄水库地下径流出流系数：一般来说，约为0.7左右，二者分配比与具体流域有关，只能调试决定。壤

25、中流消退系数和地下径流消退系数：一般趋近0.9；约为0.980.998，具体取值与具体流域有关，可进一步调试决定。河网汇流单位线参数和在0.02左右，与具体流域的河网特性有关，可进一步调试决定。与具体流域的河网特性有关，为计算方便，通常取计算单位时段的整倍数。2参数率定首先率定第一类参数，主要观察率定期的总水量是否与实测水量相近。其次率定第二类参数，主要观察次洪水的水量是否与实测水量相近。第三类参数的率定，主要观察次洪水的退水段过程的水量是否与实测水量相近。第四类参数的率定，主要观察次洪水的过程线是否拟合良好。参数率定过程中要反复对参数进行多遍调试，直至满足精度要求。2.5新安江三水源模型的V

26、B程序介绍新安江三水源，两层蒸发模型（时变单位线）VB源程序Function XingAnJang4() As Integer Dim Fno As Integer, I As Integer，Nian(2) As Integer，AX(20) As Single, Qo As Single, Dim Wm As Single, Wum As Single, Wlm As Single, k As Single, Sm As Single, b As Single, Bx As Single Dim Kss As Single, Kg As Single, CRss As Single, CR

27、g As Single, Immp As Single, F As Single, U As Single Dim P As Single, E As Single, Pe As Single, Eu As Single, El As Single, EE As Single, FR0 As Single Dim R As Single, Rs As Single, Rss As Single, Rg As Single, Q As Single, X As Single, XX As Single Dim W As Single, Wu As Single, Wl As Single, S

28、As Single, Wmm As Single, Smm As Single Dim QRss0 As Single, QRg0 As Single, Cr As Single, Cs As Single, Cx As Single, Qz As Single, Qa() As Single Qo = 5 Qo是基流。 Open e:test流域模型enshi9699parm3.txt For Input As #1 Open e:test流域模型enshi9699ensh9699.txt For Input As #5 Open e:test流域模型enshi9699end96994.tx

29、t For Output As #6 Fno = 1 ReadRecord Fno:TurnRecord Nian, AX, 11 + 1, 1: Close (Fno) Wm = AX(1): Wum = AX(2): k = AX(3): b = AX(4): Sm = AX(5): Bx = AX(6): Kss = AX(7): Kg = AX(8) CRss = AX(9): CRg = AX(10): Immp = AX(11) W = 85: Wu = 1: S = 0: FR0 = 0: Dt = 3: F = 2928: U = F / (3.6 * Dt): QRss0 =

30、 0: QRg0 = 10.1 Wmm = (1 + b) * Wm: Smm = (1 + Bx) * Sm: Wlm = Wm - Wum: Wl = W - Wu: X = 1 / (1 + b): XX = 1 / (1 + Bx) Cr = 0.025: Cx = 6 Fno = 5: n = 0: Do Until EOF(Fno): Line Input #Fno, TextLien: n = n + 1: Loop: Seek #Fno, 1 ReDim Qa(n + 1000) For I = 1 To n ReadRecord Fno: TurnRecord Nian, A

31、X, 4 + 1, 1 n1 = Nian(1): m1 = Int(AX(1): P = AX(2): E = AX(3): qb = AX(4) Pe = P - k * E If Pe 0 Then Eu = k * E: El = 0 XuManChanLiu1 Pe, R, Wm, Wmm, Wum, Wlm, W, Wu, Wl, b, X Else R = 0 EuEl P, k, E, Wlm, Wu, Wl, Eu, El End If RsRssRg1 Pe, R, S, Sm, Smm, Rs, Rss, Rg, Bx, XX, FR0, Kss, Kg, Immp Rs

32、RssRg2 Pe, R, S, Sm, Smm, Rs, Rss, Rg, Bx, XX, FR0, Kss, Kg, Immp EE = Eu + El: W = Wu + Wl HeWangZongRuLiu Rs, Rss, Rg, QRss0, QRg0, Qz, CRg, CRss, Qo, U Tao = 1 Cs = 1 - Cr * Qz 0.4 Tao = Int(Cx / (1 - Cs) / 3 Qa(I + Tao) = Cs * Qa(I) + (1 - Cs) * Qz Q = Qa(I + Tao) Write #6 I, P, EE, R, Rs, Rss,

33、Rg, W, Q, qb End If Next I Close (6) Form1.Text1.Text = CStr(666)End FunctionPublic Sp(20) As VariantFunction ReadRecord(Fno) As Integer Dim TextLien As String, Stext As String, sss As String Dim I As Integer, j As Integer, LineLen As Integer j = 1 Line Input #Fno, TextLien LineLen = Len(TextLien) F

34、or I = 1 To LineLen Stext = Mid(TextLien, I, 1) Select Case Asc(Stext) Case 45, 46, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57 sss = sss + Stext If I = LineLen Then Sp(j) = sss End If Case Else If sss Then Sp(j) = sss: j = j + 1: sss = End If End Select Next I ReadRecord = jEnd Function蓄满产流函数一Function X

35、uManChanLiu1(Pe, R, Wm, Wmm, Wum, Wlm, W, Wu, Wl, b, X) As Integer Dim Dr As Single, A As Single A = Wmm * (1 - (1 - W / Wm) X) If Pe + A Wmm Then R = Pe - Wm + W + Wm * (1 - (Pe + A) / Wmm) (1 + b) Else R = Pe + W - Wm End If If Wu + Pe - R Wlm Then Dr = Wl - Wlm Wl = Wlm R = R + Dr End If End IfEnd Function两层蒸发函数Function EuEl(P, k, E, Wlm, Wu, Wl, Eu, El) As Integer If Wu + P k * E Then Eu = k * E Wu = Wu + P - Eu El = 0 Else Eu = Wu + P El = (k * E - Eu) * Wl / Wlm Wu = 0 Wl = Wl - El End IfEnd Function三水源划分函数一Function RsRssRg1(Pe, R, S, Sm, Smm, Rs, Rss, Rg, Bx, XX, FR0, Kss, Kg, Im