第2章 分布式水循环模拟模型EasyDHM

2.1 水文模拟

2.1.1 产流理论

本书提出的EasyDHM产流模型是在Wetspa模型（Water and Energy Transfer Between Soil，Plants and Atmosphere）、SWAT模型（Soil and Water Assessment Tool）、新安江模型等产流模型的基础上，通过集成创新而提出的一种产流模型，其在不同地区、不同水文地质条件下均具有通用性，其产流过程概化见图2-1。EasyDHM模型在垂向上划分了4层，即植被冠层、地表层、土壤层和地下水含水层。首先降雨进入植被冠层，发生冠层截留，穿过冠层的水分则会进入地表层，一部分超渗的水量会发生地表填洼，一部分形成地表径流，其他水分则通过入渗进入到土壤，土壤水侧向流出形成壤中流，土壤水下渗进入地下水含水层，形成地下径流。地表径流、壤中流和地下径流之和即为总径流。植被冠层、地表填洼、土壤水和地下水均会产生蒸散发。

EasyDHM产流模型除了支持植被冠层、地表过程、土壤水过程及地下水过程模拟外，还可实现寒区水文模拟，如积雪/融雪、冻土/冻融模拟等。

2.1.1.1 冠层过程模拟

一般，冠层截留损失在降水初期相对较高，随后会逐渐趋于零。在EasyDHM产流模型中，认为冠层截留达到最大存储量时，截留率才会下降。如果在某一时段内，总降水量比冠层截留蓄水量大，则降水先会充满所有冠层截留容量，剩余降水则会继续下降；否则全部降水均认为被冠层截留，剩余的冠层截留量则在之后的时段里从降水中除去，其表达式为：

式中：Ii（t）为t时段内i单元的冠层截留损失，mm；Ii，0为单元截留蓄水能力，mm；SIi（t-1）为t-1时间段的单元截留蓄水量，mm；Pi（t）为单元降水量，mm。

在某一单元中，截留蓄水量平衡的计算公式为：

式中：SIi（t-1）和SIi（t）分别为t-1时间段和t时间段i单元的截留蓄水量，mm；EIi（t）为i单元截留量中的蒸发部分，mm。

当截留量为0或是在降水过程中，EIi（t）=0；当Pi（t）=0且EP&gt；SIi（t-1）&gt；0时，EIi（t）=SIi（t-1），其中EP为潜在蒸发量，mm；其他条件下，EIi（t）=EP。

截留蓄水能力是有关叶面积指数和植被种类的函数。由于截留蓄水能力随时间不断变化，为方便模型计算，建立了一个简单的正弦变化曲线。基于长期测量值的统计分析，其经验方程为：

式中：Ii，min为单元网格i中的最小截留蓄水能力，mm；Ii，max为最大截留蓄水能力，mm；d为日序数；指数b控制正弦变化曲线的形状，可根据局部条件进行调整。

EasyDHM模型中假定一天中每小时截留蓄水能力恒定，因此截留蓄水能力为日函数。

2.1.1.2 地表过程模拟

扣除截留后的降雨到达地表后，部分降雨将发生地表入渗、填洼、地表径流及填洼蒸发等过程，这些统称为地表过程。

在EasyDHM中，计算地表入渗量和扣除入渗量后的剩余降雨量的方法有两种：一种是与土壤含水量有关的修改系数法；另一种是SWAT模型中采用的SCS曲线法。修改系数法的思想是采用超渗产流，即在特定降雨过程中，当降雨强度超过了土壤下渗能力时产生的地面径流，用该方法计算剩余降雨量和地表入渗量公式为：

式中：PEi（t）为一定时间间隔中i网格上超渗降水量，mm；Pi（t）为t时段的降水量，mm；Ii（t）为植被冠层截留损失，mm；θi（t）为网格内t时间的土壤含水量，m3/m3；θs，i为土壤孔隙度，即饱和含水量；a为与降水强度有关的指数；Ci为网格潜在降水超渗系数或潜在径流系数；Fi（t）为单元网格的入渗量，mm。

SCS曲线法又称曲线数值法（curve number method），是美国农业部水土保持局（Soil Conservation Service，SCS）于20世纪50年代年开发研制的，其几乎可以在所有的土地利用类型和土壤类型的区域使用。该方法计算剩余降雨量和地表入渗量的公式为：

式中：St（t）为滞留参数。

St（t）的计算公式为：

式中：CNi（t）为SCS曲线参数。

地表下渗后，剩余降雨量部分进行填洼，部分成为地表径流参与汇流。某一时刻地表填洼量与地表填洼能力、地表填洼水储量及累积剩余降雨量有关，其计算经验公式为：

式中：SDi（t）为地表填洼水储量，mm；SDi，0为地表填洼能力，mm；PCi（t）为累积剩余降雨量，mm，即t时刻时计算单元i的累积剩余降雨量。EasyDHM产流模型仅根据式（2-8）模拟初始时刻的地表填洼水储量，其后按更新填洼需水量公式（2-9）计算。t时刻的地表填洼量为：

进行地表填洼后的累积剩余降雨量则为：

地表径流RSi（t）为计算单元i在计算时刻t的剩余降雨量减去地表填洼量，其计算公式为：

根据EasyDHM产流模型，认为填洼水量均通过蒸发逐渐消逝。填洼水的蒸发量计算式为：

式中：PE为潜在蒸发量，mm；EIi（t）为截留蒸发量，mm。

当发生填洼蒸发后，即需要更新填洼水储量，其计算公式为：

2.1.1.3 土壤水过程模拟

EasyDHM产流模型中，土壤水的模拟采用分层土壤模型，逐层计算土壤层的下渗和壤中流。分层土壤模型将土壤分为若干层（层数可自行设定），从上到下，假设土壤水在每一层均匀分布，当某一层的土壤含水量θi达到田间持水量θfc，i时形成重力水，而其下层土壤未饱和，该层土壤水往下层移动，当最下层土壤水蓄满后即向土壤层下的包气带渗漏，继而对地下含水层进行补给。

分层土壤模型将每个计算单元进行了垂向分层，对于任一单元的第i层土壤，其每天可下渗水量SW0，i为：

式中：θi为第i层土壤的水含量，mm；θfc，i为第i层土壤的田间持水率，mm。

另外，若第i层土壤温度较低，低于某一临界温度时，则视为冻土，冻土层土壤水不发生下渗。实际上，并不是所有可下渗水量均进入下一层土壤，第i层土壤的实际下渗量Wp，i为：

式中：Δt为计算时间步长；TTi为下渗持续时间。

TTi计算式为：

式中：θs，i为第i层土壤的饱和含水率，mm；Ki为该层土壤渗透系数，mm/s。

由于对土壤进行了分层处理，则可对土壤垂向土层分布多样化的情况进行分析，如某层土壤为不透水层，则水体下渗到此层即停止。

当某层水蓄满且该层存在一定的坡度时，则会发生横向流动，即形成壤中流，否则其上层土壤水发生横向流动，如此逐层上推，直至最上层土壤发生蓄满产流。土壤水横向流动Qlat，i计算式为：

式中：H0，i为土壤层等效水头，mm；φd，i为土壤层的排水孔隙度，即土壤层孔隙度θs，i与田间持水量θfc，i的差值；Li为土壤层横向等效长度，m；vlat，i为横向流动速度，是渗透系数的分量；αi为土壤层坡度，其值一般较小，可认为sinαi=tanαi，即sinαi为等效水头与横向距离的比值。若土壤层的一段裸露于地表，则横向流动量将贡献地表径流。

计算土壤水蒸发量时，首先应区分出不同深度土壤层需要的蒸发量。土壤深度层次的划分决定土壤允许的最大蒸发量，其计算公式为：

式中：Esoil，z为z深度处蒸发需要的水量，mm；z为地表以下土壤层中心的深度，mm；Es为最大可能土壤水分蒸发量，mm；式中的系数项是为了保证50%的蒸发所需水分来自土壤表层0～10mm，95%的蒸发所需水分来自0～100mm土壤深度范围内。

2.1.1.4 地下水径流过程模拟

储存在地下饱和区的水量即为地下水，其按含水层是否具有自由表面分为潜水含水层和承压含水层。由于承压含水层一般不发生水平方向的水量交换，并且没有蒸发，不会对地表径流产生直接影响。因此，模型仅对潜水含水层进行模拟。

当对河床状况了解较少时，流域的地下水模拟可使用线性水库或非线性水库的概念，即认为整个地下含水层是一个地下水库。地下水出流计算公式为：

式中：QGi为计算单元i的出口平均地下水出流量，m3/s；SGi为t时刻计算单元i的地下水储水量，mm；m为指数，m=1为线性水库，m=2为非线性水库；cg为考虑计算单元面积后的地下水回归系数。对每一个计算单元i，其地下水水量平衡表达式如下：

式中：SGi（t）和SGi（t-1）分别为t时刻和t-1时刻的子流域地下水储水量，mm；Fi为计算单元的面积，m2；EGi（t）为计算单元i的地下水储水量的计算平均蒸散发量，mm；QGi（t）为计算单元i的地下水出流量，m3/s；Qi，rch（t）为单元i的补给量，mm。Qi，rch（t）由上覆土壤层渗漏量Qi，srch（t）和上一时段的补给量Qi，rch（t-1）加权平均求得：

式中：δgw为有效排水时间，由地质条件决定。

EasyDHM模型中，地下水蒸发量公式为：

式中：EGi（t）为计算单元i的地下水储水量的平均蒸散发量，mm；EP为潜在蒸散发量，mm；cv为作物系数；EIi（t）为计算单元i的截留水蒸发量，mm；EDi（t）为计算单元i的填洼水蒸发量，mm；ESi（t）为计算单元i的土壤水蒸发量，mm；cd是一个可变参数。cd的计算公式为：

式中：SGi（t）为t时刻计算单元i的地下水储水量，mm；SGs，0为计算单元i的地下水储水量的最大值，mm。

2.1.1.5 潜在蒸散发计算

EasyDHM模型提供了Penman—Monteith、Priestley—Taylor和Hargreaves三种计算潜在蒸散发能力的方法，另外还可以读入实测的水面蒸发能力等直接得到的逐时段的潜在蒸散发数据［145］。本研究采用世界气象组织推荐的无论在干旱还是湿润地区计算精度均较高的Penman—Montieth公式法计算潜在蒸发量。

Penman-Montieth公式需要输入的资料为辐射、气温、风速和相对湿度，其公式为：

式中：ET0为潜在蒸散发能力，mm；Δ为饱和水汽压斜率，kPa/℃；Rn为净辐射量，MJ/m2；G为土壤热通量，MJ/m2；ρ为空气密度，g/m3；D为饱和水汽压差，kPa；ra为边界层阻力，s/m；L为汽化潜热，MJ/kg；γ为湿度计常数。

空气密度ρ计算公式为：

式中：T为气温，℃；PB为大气压力，kPa。PB的计算公式为：

式中：Elev为计算点的高程，m。

边界层阻力ra的公式：

式中：zd为零平面位移高度，m；zd=0.702H0.979，其中H为冠层高度，m；zo为蒸散面粗糙长度，m，zo=0.131H0.997；V为日均风速，m/s。

2.1.2 汇流理论

本次研究采用马斯京根法对河道汇流进行演算，该法模拟了沿渠道长度柱蓄（Prism Storage）和楔蓄（Wedge Storage）组成的蓄水容量。河段槽柱蓄与楔蓄如图2-2所示。

洪水波行进到某个河段槽，当入流量大于出流量便形成了楔形蓄水体；当洪水波退去，在河段槽便出现了出流量大于入流量的负楔蓄。另外，对于楔蓄水体，河段槽内始终包含一个体积为流域长度上横截面不变的柱蓄水体。

总的蓄水容量为：

式中：Vstored为蓄水容量；qin为入流量；qout为出流量；K为稳定流情况下的河段传播时间；X为流量比重因素。该公式可重新整理为：

流量比重因素X的下限为0.0，上限为0.5。这个因子是楔蓄量的函数。对于水库式蓄水，没有楔蓄，X=0.0；而对于一个完全的楔蓄，X=0.5；对于河流，X落在0.0～0.3之间，其平均值接近0.2。

对于蓄水容量的定义可以加入连续公式并简化为：

式中：qin，1为该时间段开始时的入流量；qin，2为该时间段结束时的入流量；qout，1为该时间段开始时的出流量；qout，2为该时间段结束时的出流量。

其中

C1+C2+C3=1

为用体积单位表示所有值，式两边都要乘以该时间段：

为了保持数值稳定和避免负出流量的计算，必须满足以下条件：

流量比重因素X的值由使用者输入，蓄水时间常数的值估计如下：

式中：K为稳定流情况下的河段传播时间，s；coef1和coef2为由使用者输入的权重系数；Kbnkfull为稳定流情况下渠道蓄满水的河段传播时间；K0.1bnkfull为渠道蓄满1/10水量时河段传播时间。Cunge（1969）提出了一个计算Kbnkfull和K0.1bnkfull的公式：

式中：K为稳定流情况下的河段传播时间，s；Lch为渠道长度，km；ck为指定深度处的波速。波速ck计算公式为：

式中：流速qch由曼宁公式定义。

2.1.3 水库调度

模型将分布式水文模型与水库调度模型相结合，实现上游来水预报与水库调度演算的耦合。因此，研究将分布式水文模型EasyDHM与水库调度模型充分结合，开发水库群联合调度模型的同时也对分布式水文模型EasyDHM进行了改进。考虑水库的兴利、防洪调度作用后，水库节点的水文过程能得到更合理的重现，并能对下游河道的径流过程产生影响。此外，由于分布式水文模型存有水系、水文站、水库的空间拓扑关系，能很方便地描述水库群上下游间的影响关系，这也是基于分布式水文模型的水库群联合调度模型的优点所在。分布式水文模型与水库群调度模型耦合过程如图2-3所示，共有三种耦合方式，依据水库群调度模型中是否包含产流计算和汇流计算分为三类，即产流—汇流—水库调度的紧密耦合，汇流—水库调度的简化，以及两模型分开模拟的离线耦合方式。

分布式水文模型紧密耦合水库群调度模块。这种方式是直接在分布式水文模型中开发一个水库调度模块，实现了水文模拟、水库调度的紧密耦合。使用时，各子流域首先进行子流域内部计算单元的产流计算；然后进行子流域内部的汇流计算，汇流至子流域出口点后判断子流域出口点是否存在水库，如果存在水库则进行水库调度计算，得到水库出流作为子流域的出流汇入到下游子流域，否则直接汇入下游子流域，水库调度模型与水文模型所用时间步长是完全一致的。其优点是使用方便，在计算水文模拟的同时计算完水库调度，同步输出各模型结果，方便进行比较。缺点是模型计算时占据电脑内存较大，计算速度相对慢；同时，水文模型输出的中间结果数据量较大，从水库调度角度来看并不需要，增加了额外存储空间。

基于分布式水文模型区间产流的汇流/水库调度简化耦合模型。和上述方式不同，它是一个耦合河道汇流模型和水库调度模型的水库群调度模型，需事先输入水文模型的产流模拟结果。这种方式中，汇流模型和水库调度模型时间步长必须一致，但可以与水文模型（即模拟产流过程）的时间步长不一致。如果时间步长不一致，则需在使用水文模型的产流结果时预先进行时间转化。这种方式既能节省模型计算时间也能准确考虑河道汇流过程。因此，这种方式的水库群调度模型可以直接利用分布式水文模型率定好的汇流参数，避免单独率定水库群之间河道参数的困难，可以方便地应用到水库群实时调度中去。

分布式水文模型与水库群调度模型的离线耦合方式。这种方式的水库群调度模型和传统的水库群调度模型一样，但区间来水（区间入流）是从分布式水文模型中得出的。首先，借助分布式水文模型EasyDHM计算出各水库的天然径流过程，依据各方案水库间的拓扑关系以及各水库的天然径流过程即可得出各水库的区间入流量，即区间入流量为区间产流量减去河道槽蓄量及损失量。由于水库群调度模型不包含产流模块和汇流模块，在研究中进行相应简化，给出以下假定：尽管考虑上游水库的调蓄作用后，上游分区汇入的某水库分区的流量会发生改变，但认为这种流量改变对该分区的河道槽蓄作用（包含蒸发、渗漏过程）影响不大，即河道的槽蓄过程不随河道入口点输入流量的改变而发生变化。这种方式是一种水库调度模型与分布式水文模型间离线耦合的方式，水库调度模型与分布式模型的时间步长同样可以不同。

总的来说，以上三种方式中，第一种方式计算最为耗时，但使用最为方便，对整个流域水循环模拟过程考虑得更为详尽，计算结果也更为精确；第二种方式是为节省第一种方式计算时间而改进的，因为在同等的气象及参数条件下，区间产流结果是不会发生变化的，改变水库调度规则或者设置不同水库情景改变的只是河道汇流过程；第三种方式则为纯粹的水库群调度模型，通过输入各水库区间入流量进行自上而下的水库群调度模拟，水库间的河道汇流过程暗含在区间入流中。