2.2　分布式水文模型

为了研究人工取用水变化和下垫面变化等区域人类活动对水资源的影响，在变化环境下流域水资源演变的归因方法中选用WEP-L（Water and Energy transfer Processes in Large river basins）（贾仰文，2005）模型来模拟变化环境下的水资源情况。

WEP（Water and Energy transfer Process）模型的开发始于1995年，1999—2002年，经过改进和完善，逐渐趋于成熟，在日本、韩国的多个流域得到了验证和应用（Jia 等，2001、2005），并于2002年10月获日本国著作权登录。2003年，针对我国内陆河流域的特点，对WEP模型进行了改进，形成IWHR-WEP模型。2003—2004年，在网格单元型WEP模型的基础上，开发了耦合模拟天然水循环过程与人工侧支水循环过程的大尺度流域分布式水文模型WEP-L，经过模型验证后，应用于黄河流域水资源评价和人类活动影响下的水资源演变规律研究（Jia等，2006）。

WEP-L模型在IWHR-WEP模型的基础上又做了如下改进：

（1）为了适应超大规模流域水文模拟的需要，改网格单元为子流域套等高带单元，既加快了计算速度，又避免了汇流失真。

（2）针对各水循环要素过程的特点，采用变时间步长进行模拟计算（强降雨时期入渗产流过程采用1h、坡地与河道汇流过程采用6h，其余的过程采用1d），既保证了水循环动力学机制的合理表述，又提高了计算效率。

（3）将灌区引水口、水库落实到具体的计算河段，将灌溉用水落实到田间，用水类别分农业、工业、生活三大类，进一步划分为七小类，按地表水、地下供水分别计算，加强了人工用水系统与自然水循环系统的耦合模拟。

目前，WEP-L模型已在海河流域、黄河流域、松辽流域、汉江流域等进行了广泛应用。

2.2.1　模型结构

WEP-L模型各基本计算单元内的垂向结构如图2.1（a）所示。从上到下依次为植被或建筑物截留层、地表洼地储留层、土壤表层、过渡带层、浅层地下水层和深层地下水层（承压层）等。状态变量包括植被截留量、洼地储留量、土壤含水率、地表温度、过渡带层储水量、地下水位及河道水位等。主要参数包括植被最大截留深、土壤渗透系数、土壤水分吸力特征曲线参数、地下水透水系数和产水系数、河床透水系数和坡面、河道糙率等。为了考虑计算单元内土地利用的不均匀性，采用了“马赛克”法，即把计算单元内的土地归成若干类，分别计算各类土地类型的地表面水热通量，再取其面积平均值作为计算单元的地表面水热通量。土地利用首先分为裸地-植被域、非灌溉农田、灌溉农田、水域和不透水域5个大类。裸地-植被域又分为裸地、草地和林地3类，不透水域分为城市地面与城市建筑物2类。另外，为了反映表层土壤的含水率随深度的变化，便于描述土壤蒸发、草或作物根系吸水和树木根系吸水，将透水区域的表层土壤分成3层。

WEP-L模型的平面结构为子流域套等高带，如图2.1（b）所示，其中①～⑨为子流域编号，Q1～Q9为各子流域的坡面汇流量。前述WEP-L模型下垫面划分及垂向结构主要用于产流计算，产流计算完成后按照平面结构进行汇流过程演算。子流域内各等高带之间进行坡面汇流演算，从最高等高带开始顺坡而下，最终汇入河道。各子流域之间进行河道汇流演算，从流域最上游开始依次向下直到流域出口。坡面汇流计算根据各等高带的高程、坡度与Manning糙率系数，采用一维运动波法将坡面径流由流域的最上游端追迹计算至最下游端。各条河道的汇流计算，根据有无下游边界条件采用一维运动波法或动力波法由上游端至下游端追迹计算。地下水流动分山丘区和平原区分别进行数值解析，并考虑其与地表水、土壤水及河道水的水量交换。

除了自然水循环过程之外，WEP-L模型还进行人工水循环过程的模拟。每个计算单元内使用经验统计方法进行人工水循环模拟，需要通过外部输入各计算单元内相关社会用水量（包括农业灌溉用水量以及工业、生活用水量），进行“自然-社会”二元水循环过程模拟，以反映人类活动影响下的流域水循环过程。其中，社会用水量可以是历史数据，也可以是规划预测数据。输入的社会用水量都必须展布到各计算单元上，以反映社会用水的空间变化。

2.2.2　自然水循环过程模拟

在流域内生成“子流域套等高带”的计算单元后，根据改进的Pfafstetter规则对每个计算单元建立包含拓扑关系的编码。根据此编码，在流域内从上游的子流域到下游子流域依次计算；而在子流域内按等高带由高至低顺序计算。自然水循环过程模拟主要包括蒸发蒸腾、入渗、地表径流、壤中流、地下水运动、地下水出流和地下水溢出、坡面和河道汇流以及积雪融雪过程。而降水等气象要素作为水文模型的主要驱动因素之一，将其展布到分布式水文模型的计算单元上，是进行自然水循环过程模拟的前提。

2.2.2.1　气象要素数据展布

1. 气象数据的空间展布

WEP-L模型是日尺度模型，需要逐日气象要素数据作为输入。气象数据一般来自流域内的气象站、雨量站，属于点上的数据，而水文模型的计算单元是等高带，属于面上的数据，因此需要通过空间展布手段将点上数据转化为面上数据。

空间展布也可以理解为空间插值，在选择不同的方法时，精度是首先需要考虑的因素，然而计算效率也不能忽视。由于所研究的流域空间尺度一般较大，一般水文模型计算的时间系列又较长，所以在选用插值方法时，应该尽量避免采用计算量特别大的方法。WEP-L模型中采用距离平方反比法和泰森多边形法相结合的方法，以距离平方反比法为主，少量地区采用泰森多边形法。

（1）泰森多边形法。泰森多边形法是一种广泛使用的空间插值方法，该方法一般需要绘制泰森多边形，以每个多边形所包含的雨量站点的降水量值作为该区域内各点的降水量值，其实质是平面上每个点取距离最近的站点的实测值。如果用计算机实现，可以先把研究区域划分为网格，计算每个站点与该网格中心的距离，取最近的雨量站点的观测值作为该网格的降水量。

（2）距离平方反比法（Reversed Distance Squared，RDS）。 假设待估点的降水可以用它周围的一些雨量站插值得到，同时假设待估点的降水量和雨量站点的降水量大小成正比，与雨量站点的距离成反比。根据此假设，待估点的雨量估计值可以表示为

式中：P为待估点的降水量；Pi为雨量站的降水量；di为雨量站到待估点的距离；IR为用来插值的雨量站点集合；n为权重系数，一般取2时是距离平方反比法，取0时是算术平均法。

由泰森多边形法和距离平方反比法的计算公式可见，两类方法概念清楚，计算简便，计算速度快，且距离平方反比法能够考虑降水空间趋势性，两种方法结合，适合于大尺度流域空间插值计算。

泰森多边形法的参证站为最近的一个站点，然而距离平方反比法需要选择一定数量的站点作为参证站。目前选择站点的方法一般有两类：一是固定参证站点个数m，即选择离待估点最近的m个站点进行插值；二是固定距离D，即选择离待估站点距离小于D的站点作为参证站。

由于本次收集到的海河流域降水气象站空间分布不均，如果采用第一类选择站点的方法，对于雨量站点分布密集的地方，选择的站点离待估点很近，插值效果很好；而对于雨量站点分布稀疏的地方，可能选择的参证站离待估点很远，两处在气象上关系并不密切，插值效果很难保证。如果采用第二类选择站点的方法，对于雨量站点分布密集的地方，可能有很多站点都会入选，而对于雨量站点分布稀疏的地方，可能得到的站点很少，甚至一个站点也没有。由于站点的分布不均匀性，采用固定个数或距离的方式选取参证站有时是不合理的，需要采取一种比较灵活的、有弹性的方法。

考虑到不管采用什么插值方法，研究者都希望参证站的降水量和待估点的降水量相关性比较好，二者的相关系数可以用来判断空间各点之间雨量相关性的好坏。因此，本研究引入相关系数作为选取参证站点的指标。但是对于站点比较稀疏或者影响降水气象的因素特别复杂的地方，可能某个站点和其他所有站点的相关系数都小于该阈值，这样就没有相关站点了。对于这种情况，采用任何方法效果都不会好，为便于计算，则采用最简单的泰森多边形法。这就是本次研究所采用的考虑相关系数的综合插值方法（ARDS，Amended RDS）。

采用ARDS进行降水、气象要素插值的具体计算步骤为：①两两计算所有站点之间的相关系数；②确定一个相关系数的阈值，针对每个站点，根据阈值判断与之相关的所有站点，并将最远一个相关站点与之相间的距离作为最大相关距离；③对每个子流域，计算该子流域的形心与每个降水气象站点之间的距离，如果小于某个站点的最大相关距离，该站点即为该子流域的参证站点；④如果某一子流域存在相关站点，采用距离平方反比法进行降水气象要素进行插值，如果不存在相关站点，采用泰森多边形法进行插值。

有关前述插值方法的详细介绍请参考（周祖昊等，2006）。

2. 降雨数据的时间降尺度展布

WEP-L模型进行产流计算时除采用日尺度外，还釆用小时尺度（主要在暴雨期）。一般提供的降雨资料都是逐日数据，需要降尺度到小时数据。当然，也可以直接使用小时数据进行计算，但受降水资料精度限制以及数据系列长度限制，某些年份根本无法获取小时数据，此时需要进行降尺度处理。

WEP-L模型采用分区建立日雨量向下尺度化模型，将插值所得日降雨量进行向下尺度化。由于大强度降雨的日内分布对流域产流影响较大，所以仅考虑大于10mm的日降雨进行降尺度处理。实际上，在WEP-L中以10mm为界分为暴雨期和非暴雨期，且不同时期釆用不同的入渗计算方式。暴雨期采用Green-Ampt模型模拟，而非暴雨期则根据水量平衡按饱和导水率计算，详细见2.2.2.3节的入渗过程模拟部分。

向下尺度化模型假设每个大强度降雨日内只有1次，具体计算公式为

式中：i为历时t内最大降水平均雨强；S为暴雨参数（或称之为雨力），等于单位时间内最大平均雨强；t为时段；n为暴雨衰减系数，与气候区域有关，可用实测资料率定获取；P为日降雨量；T为日降雨总历时；a、b为参数。

使用日雨量向下尺度化时需要先对研究区域使用实测小时降雨量率定参数a、b和n。日内分配首先根据式（2.14）由日降雨计算雨力参数S，然后根据式（2.13）计算当日降雨总历时T，最后由式（2.12）计算降雨历时内每个小时的降雨量。

2.2.2.2　蒸发蒸腾模拟

蒸发蒸腾是水循环过程中的一个重要环节，主要通过改变土壤的前期含水率来影响降水入渗，进而影响产流，是生态用水和农业节水等应用研究的对象，因此准确计算蒸发蒸腾具有重要的意义。WEP-L模型采用了“马赛克”结构，考虑计算单元内的土地利用变异性，每个计算单元内的蒸发蒸腾可能包括植被截留蒸发、土壤蒸发、水面蒸发和植被蒸腾等多项，按照土壤-植被-大气通量交换方法（SVATS）、采用Noilhan-Planton模型、Penman公式和Penman-Monteith公式等进行详细计算。蒸发蒸腾过程往往伴随着能量交换过程，计算蒸发蒸腾，必须计算地表附近的辐射、潜热、显热和热传导，而这些热通量又都是地表温度的函数。为减轻计算负担，热传导及地表温度的计算采用强制复原法，详见2.2.3节的能量循环过程模拟部分。

1.阻抗参数计算

（1）空气动力学阻抗。地表面附近大气中的水蒸气及热的输送遵循紊流扩散原理。近似中立大气的空气动力学阻抗（ra）的计算公式（Jia，1997）为

式中：z为风速、湿度或温度的观测点离地面的高度；к为 Karman常数即0.41；U为风速；d为置换高度；zom、zox分别为风速、比湿的粗糙度。

根据Monteith理论，若植被高度为hc，则zom= 0.123hc，zox= 0.1zom，d=0.67hc。

大气安定或不安定时，运动量、水蒸气及热输送还受浮力的影响，空气动力学阻抗需根据Monin-Obukhov相似理论计算。

（2）植被群落阻抗。植被群落阻抗是各个叶片的气孔阻抗的总和，Dickinson等（1991）提出了计算公式，即

式中：LAIi为n层植被的第i层的叶面指数；rsi为第i层的叶气孔阻抗；（rs）为群落的气孔阻抗的平均值；rsmin为最小气孔抵抗；f1为温度的影响函数；f2为大气水蒸气压的饱和差[饱和水蒸气压与大气的水蒸气压的差VPD（Vapor Pressure Deficit）]的影响函数；f3为光合作用有效放射PAR（Photosynthetically Active Radiation flux）的影响函数；f4为土壤含水率的影响函数。

若忽视LAI对叶气孔阻抗（rs）的影响，则可得到以下公式（Dickinson，1991）：

式中：Ta为气温，℃；VPDc为叶气孔闭合时的VPD值（约为4kPa）；PARc为PAR的临界值（森林为30W/m2 ；谷物为100W/m2）；rsmax为最大气孔阻抗，为5000s/m；θ为根系层的土壤含水率；θw为植被凋萎时的土壤含水率（凋萎系数）；θc为无蒸发限制时的土壤含水率（临界含水率）。

2.蒸发蒸腾量计算

计算单元内的蒸发蒸腾包括来自植被湿润叶面（植被截留水）、水域、土壤、城市地表面、城市建筑物等的蒸发，以及来自植被干燥叶面的蒸腾。计算单元的平均蒸发蒸腾量由下式通过面积加权平均算出：

式中：Fw、Fu、Fsv、Fir、Fni分别为计算单元内水域、不透水域、裸地-植被域、灌溉农田及非灌溉农田的面积率；Ew、Esv、Eu、Eir、Eni分别为计算单元内水域、不透水域、裸地-植被域、灌溉农田及非灌溉农田的蒸发量或蒸发蒸腾量。

水域的蒸发量（wE）由Penman公式（Penman，1948）算出，同时也用于计算区域蒸发能力：

式中：RN为净放射量；G为传入水中的热通量；Δ为饱和水蒸气压对温度的导数；eΔ为水蒸气压与饱和水蒸气压的差；ra为蒸发表面的空气动力学阻抗；ρa为空气的密度；cp为空气的定压比热；λ为水的气化潜热；。

裸地-植被域的蒸发蒸腾量（Esv）计算公式为

式中：Ei为植被截留蒸发（来自湿润叶面）；Etr为植被蒸腾（来自干燥叶面）；Es为裸地土壤蒸发；下标1为高植被（森林、城市树木）；下标2为低植被（草）。

植被截留蒸发（Ei）使用Noilhan-Planton模型（Noilhan和Planton，1989）计算：

式中：Veg为裸地-植被域的植被面积率；δ为湿润叶面的面积率；Ep为可能蒸发量（由Penman方程式计算）；Wr为植被截留水量；P为降雨量；Rr为植被流出水量；Wrmax为最大植被截留水量；LAI为叶面积指数。

植被蒸腾由Penman-Monteith公式（Monteith，1973）计算：

式中：RN为净放射量；G为传入植被体内的热通量；rc为植被群落阻抗（canopy resistance）。

蒸腾属于土壤、植被、大气连续体（Soil-Plant-Atmosphere Continuum，SPAC）水循环过程的一部分，受光合作用、大气湿度、土壤水分等的制约。这些影响通过式（2.32）中的植被群落阻抗（rc）来考虑。

植被蒸腾是通过根系吸水由土壤层供给。根系吸水模型参见雷志栋（1988）的《土壤水动力学》。假定根系吸水率随深度线性递减、根系层上半部的吸水量占根系总吸水量的70％，则

式中：Etr为蒸腾；ℓr为根系层的厚度；z为离地表面的深度；Sr （z）为深度z处的根系吸水强度；Etr （z）为从地表面到深度z处的根系吸水量。

灌溉农田和非灌溉农田的作物蒸腾计算与裸地-植被域的计算类似。根据以上公式，只要给出植被根系层厚度，即可算出其从土壤层各层的吸水量（蒸腾量）。本研究认为草与农作物等低植被的根系分布于土壤层的一层、二层，而树木等高植被的根系分布于土壤层的所有三层。结合土壤各层的水分移动模型，即可算出各层的蒸腾量。

裸地土壤蒸发由修正Penman公式（Jia，1997）计算，即

式中：β为土壤湿润函数或蒸发效率；θ为表层（一层）土壤的体积含水率；θfc为表层土壤的田间持水率；θm为单分子吸力（pF6.0～7.0）对应的土壤体积含水率（Nagaegawa，1996）。

不透水域的蒸发和地表径流用计算公式为

式中：P为降雨；Hu为洼地储留；Eu为蒸发；Ru为表面径流；Humax为最大洼地储留深；Eumax为潜在蒸发（由Penman公式计算）；c为城市建筑物在不透水域的面积率；下标1为城市建筑物；下标2为城市地表面。

总体而言，以上所计算的蒸发蒸腾量均为理论值，实际蒸发蒸腾量需要根据下垫面实际水量而定，如果水量充足则使用计算值，否则按实际水量计算。

2.2.2.3　入渗过程模拟

降雨时的地表入渗过程受雨强和非饱和土壤层水分运动所控制。非饱和土壤层水分运动的数值计算既费时又不稳定，许多研究也表明，除坡度很大的山坡以外，降雨过程中土壤水分运动以垂直入渗占主导作用，降雨之后沿坡向的土壤水分运动才逐渐变得重要。因此，为提高模型模拟速率，WEP-L模型将入渗过程划分为两种情况进行模拟：暴雨期和非暴雨期。划分标准是日降雨量是否超过10mm，其中暴雨期采用Green-Ampt入渗模型按小时进行模拟计算，只考虑土壤水的垂向运动，小时降雨量使用日降水通过降尺度获得；而非暴雨期，由于雨量相对较小，使用水量平衡原理进行日尺度模拟，考虑土壤水分垂向和水平向运动，土壤入渗能力按饱和导水系数计算。

Green-Ampt模型假定在入渗过程中存在一个湿润锋将土壤层划分为上部饱和部分和下部非饱和部分，应用达西定律和水量平衡原理进行计算。入渗模型物理概念明确，所用参数可由土壤物理特性导出，已经得到大量应用验证。Mein-Larson（1973）将Green-Ampt入渗模型应用于均质土壤降雨时的入渗计算，Moor-Eigel（1981）将Green-Ampt入渗模型扩展到稳定降雨条件下的二层土壤的入渗计算。考虑到由自然作用和人类活动（如农业耕作）等引起的土壤分层问题，Jia 和Tamai于1997年提出了实际非稳定降雨条件下的多层土壤入渗Green-Ampt模型，以下称通用Green-Ampt模型。WEP-L采用该模型进行计算，如图2.2所示。

当入渗湿润锋到达第m土壤层时入渗能力计算公式为

式中：f为入渗能力，mm/h；km为m土壤层的导水系数，mm/h；Am-1为上面m-1层土壤层总共可容水量，mm；Bm-1为上面m-1层土壤层因各层土壤含水率不同而引起的误差，mm；F为累积入渗量，mm，其计算方法根据地表面有无积水而不同。

如果自入渗湿润锋进入第1m-土壤层时起地表面就持续积水，那么累积入渗量由式（2.45）计算；如果前一时段tn-1地表面无积水，而现时段tn地表面开始积水，那么由式（2.46）计算。

其中

Δθ=θs-θ0

式中：Fm-1为第m-1层累积入渗量，mm；Fp为地面开始积水时的累积入渗量，mm；m为目标入渗土壤层；ki为第i层土壤层导水系数，mm/d；SW为入渗湿润锋处的毛管吸引压所引起的入渗量，mm；θs为土壤层的含水率；θ0为土壤层的初期含水率；t为时刻；tp为积水开始时刻，s；Ip为积水开始时的降雨强度，mm；tm-1为入渗湿润锋到达第m层与第m-1层交界面的时刻，s；L为入渗湿润锋离地表面的深度，mm；Li为第i层的土壤厚度，mm。

2.2.2.4　地表径流

影响流域径流产生的主要因素包括降雨、地表覆盖、地形特征、土壤类型以及人类活动等。根据土壤径流来源的不同，可以将流域径流划分成地表径流、壤中流以及地下径流三种类型。其中，地表径流为经地表直接进入河道的部分，壤中流是地下水位以上包气带中重力自由水横向移动产生的部分，地下径流则是地下水补给产生的部分。地表径流和壤中流随降雨量呈波动变化，而地下径流相对比较稳定。

从径流产生的方式来看，又可分为超渗产流（霍顿坡面径流）、蓄满产流（饱和坡面径流）、回归流以及基流。超渗产流是指当降雨强度超过了土壤入渗能力时无法下渗的部分，多发生在渗透能力较低的土壤上（如干旱半干旱区域、城镇不透水域等）。蓄满产流是指在降雨将土壤水蓄满后导致后续降水无法入渗而形成的径流，多发生在植被覆盖较好的湿润地区。回归流则是溢出地面且从地表汇入河道的壤中流。基流则是经土壤直接进入河道的地下径流。详细如图2.3所示（刘佳嘉，2013）。区分超渗产流和蓄满产流的关键在于判别降雨强度同土壤入渗能力之间的大小，前者产生过程中土壤水由上而下逐渐饱和，而后者由下而上逐渐饱和。蓄满产流和回归流的区别在于，蓄满产流是因土壤水饱和而无法下渗的降水，而回归流则是因土壤饱和而溢出的土壤水。

由于山坡地形及土壤特性的差异性，同一个山坡上不会只有一种产流方式。一般而言，山坡上部以超渗产流为主，山坡底部则以蓄满产流、回归流为主。且不同产流机理之间可以相互转化，例如在超渗产流区域，当土壤水饱和后则自动转换成蓄满产流。因此，对流域产流模拟过程中需要综合考虑各种产流机制（即混合产流）。

在WEP-L模型中，水域的地表径流等于降雨量减去蒸发量，不透水域的地表径流按上述式（2.37）～式（2.43）计算，裸地-植被域（透水域）的地表径流则根据降雨强度是否超过土壤的入渗能力分以下两种情况计算。

1.霍顿坡面径流（Hortonian overland flow）

当降雨强度超过土壤的入渗能力时将产生这类地表径流1ieR，即超渗产流，其计算公式为

式中：P为降水量，mm；Hsv为裸地-植被域的洼地储留，mm；Hsvmax为最大洼地储留深，mm；Esv为蒸散发，mm；fsv为由通用Green-Ampt模型算出的累积入渗量，mm。

2.饱和坡面径流（Saturation overland flow）

对于河道两岸及低洼的地方，由于地形的作用，土壤水及浅层地下水逐渐汇集到这些地方，土壤饱和或接近饱和状态后遇到降雨便形成饱和坡面径流，即蓄满产流。此时，Green-Ampt模型已不再适用，需根据非饱和土壤水运动方程来求解。为减轻计算负担，将表层土壤分成若干层，利用非饱和状态的达西定律和连续方程进行计算。

地表洼地储留层：

土壤表层：

土壤中层：

土壤底层：

中间参数：

式中：Hs为洼地储留，mm；Hsmax为最大洼地储留，mm；Veg1、Veg2为裸地-植被域的高植被和低植被的面积率；Rr1、Rr2为从高植被和低植被的叶面流向地表面的水量，mm；Q为重力排水，mm；QDj，j +1为吸引压引起的j层与j+1层土壤间的水分扩散，mm；E0为洼地储留蒸发，mm；Es为表层土壤蒸发，mm；Etr 为植被蒸散（第一下标中的1表示高植被、2表示低植被；第一下标表示土壤层号）；R2为壤中流，mm；k（θ）为体积含水率θ对应的土壤导水系数，mm/d；Ψ（θ）为体积含水率θ对应的土壤吸引压，kPa；d为土壤层厚度；W为土壤的蓄水量，mm；W10为表层土壤的初期蓄水量，mm。另外，下标0、1、2、3分别为洼地储留层、表层土壤层、第2土壤层和第3土壤层。

3. 壤中流

在山地丘陵等地形起伏地区，同时考虑坡向壤中流及土壤渗透系数的各向变异性。壤中流包括从山坡斜面饱和土壤层中流入溪流的壤中流，以及从山间河谷平原不饱和土壤层流入河道的壤中流两部分。第一部分的计算类似地下水出流计算，而从山间河谷平原不饱和土壤层流入河道的壤中流计算公式为

式中：k（θ）为体积含水率θ对应的沿山坡方向的土壤导水系数，mm/d；slope为地表面坡度；L为计算单元内的河道长度，m；d为不饱和土壤层的厚度，m。

4. 地下水运动、地下水出流和地下水溢出

地下水运动按多层模型考虑。将非饱和土壤层的补给、地下水取水及地下水流出（或来自河流的补给）作为源项，按照Bousinessq方程进行浅层地下水二维数值计算。 在河流下部及周围，河流水和地下水的相互补给量根据其水位差与河床材料的特性等按达西定律计算。另外，为了考虑包气带层过厚可能造成的地下水补给滞后问题，在表层土壤与浅层地下水之间设一过渡层，用储流函数法处理。

浅层（无压层）地下水运动方程：

承压层地下水运动方程：

式中：h为地下水位（无压层）或水头（承压层），m；C为储留系数；k为导水系数，m/d；z为含水层底部标高，m；D为含水层厚度，m；Q3为来自不饱和土壤层的渗透量，m；WUL为管道渗漏，m；RG为地下水出流，m；E为蒸发蒸腾，m；Per为深层渗漏，m；GWP为地下水开采量，m；下标u和1分别表示潜水层和承压层。

地下水出流的计算根据地下水位（hu）和河川水位（Hr）的高低关系，地下水流出或河水渗漏计算公式为

式中：kb为河床土壤的导水系数，mm/d；Ab为计算单元内河床处的浸润面积，m2；db为河床土壤的厚度，m。

在低洼地，地下水上升后有可能直接溢出地表。出现这种情况时，则令地下水位等于地表标高，多余地下水蓄变量计为地下水溢出。

5. 坡面汇流和河道汇流

对于坡面汇流的计算，由于超渗产流的存在，加上沟壑溪流的汇流也可以用等价坡面汇流近似，WEP-L模型采用基于数字高程模型（DEM）的运动波模型来计算。利用DEM和GIS工具，按最大坡度方向定出各计算单元的坡面汇流方向，并定出其在河道上的入流位置。

对于河道汇流的计算，根据DEM并利用GIS工具，生成数字河道网，根据流域地图对主要河流进行修正。搜集河道纵横断面及河道控制工程数据，根据具体情况按运动波模型或动力波模型进行一维数值计算。

运动波方程：

式中：A为流水断面面积，m2；Q为断面流量，m3/s；qL为计算单元或河道的单宽流入量（包含计算单元内的有效降雨量、来自周边计算单元及支流的水量），m3/（s·m）；n为Manning糙率系数；R为水力半径，m；S0为计算单元地表面坡降或河道的纵向坡降；Sf为摩擦坡降。

动力波方程（Saint Venant 方程）：

式中：V为断面流速，m/s；Vx为单宽流入量的流速在x方向的分量，m/s；其余符号意义同前。

6. 积雪融雪过程

尽管“能量平衡法”对积雪融雪过程的描述提供了很好的物理基础，但由于求解能量平衡方程所需参数和数据过多，因此在实践中常用简单的“温度指标法”（Temperature-index approach）或“度日因子法”来模拟积雪融雪的日或月变化过程。WEP-L模型目前采用“温度指标法”计算积雪融雪的日变化过程，积雪融雪计算在产流过程之前进行，通过气温、初始积雪量决定是否进行该过程，计算公式为

式中：SM为融雪量，mm/d；Mf为融化系数或称“度日因子”，mm/（℃·d）；Ta为气温指标，℃；T0为融化临界温度，℃；S为积雪水当量，mm；SW为降雪水当量，mm；E为积雪升华量，mm。

“度日因子”既随海拔高度和季节变化，又随下垫面条件变化，常作为模型调试参数。一般情况下为1～7mm/（℃·d），且裸地高于草地、草地高于森林。气温指标通常取日平均气温。融化临界温度通常为-3～0℃。另外，为将降雪与降雨分离，还需要雨雪临界温度参数（通常为0～3℃）。

2.2.3　能量循环过程模拟

蒸发蒸腾与能量循环过程密切相关，WEP-L模型对地表面-大气间的能量循环过程进行了比较详细的模拟。地表面的能量平衡方程为

式中：RN为净放射量，MJ/m2；Ae为人工热排出量，MJ/m2；LE为潜热通量，MJ/m2；H为显热通量，MJ/m2；G为地中热通量，MJ/m2。

净放射量（RN）是短波净放射量（RSN）和长波净放射量（RLN）相加求得：

WEP-L模型包括日以内的能量循环过程模拟与日间平均能量循环过程模拟两个模块，本书只对日平均能量循环过程模拟部分做简要介绍，包括短波放射、长波放射、潜热通量、地中热通量、显热通量和人工热排出量。

1.短波放射

在没有短波放射观测数据的情况下，通常利用日照时间观测数据推算日短波放射量。推算公式（Jia，1997）为

式中：RS为到达地表面的短波放射量，MJ/（m2·d）；α为短波反射率；RS0为太阳的地球大气层外短波放射量，MJ/（m2·d）；as为扩散短波放射量常数（在平均气候条件下为0.25）；bs为直达短波放射量常数（在平均气候条件下为0.5）；n为日照小时数；N为可能日照小时数、dr为地球与太阳之间的相对距离；ωs为日落时的太阳时角；φ为观测点纬度（北半球为正，南半球为负）；δ为太阳倾角；J为Julian日数（1月1日起算）。

2.长波放射

在没有长波放射观测数据的情况下，日长波放射量的推算公式（Jia，1997）为

式中：RLN为长波净放射量，MJ/（m2·d）；RLD为向下（从大气到地表面）长波放射量，MJ/（m2·d）；RLU为向上（从地表面到大气）长波放射量，MJ/（m2·d）；f为云的影响因子；ε为大气与地表面之间的净放射率；σ为Stefan-Boltzmann常数，即4.903×10-9 MJ /（m2 ·K4 ·d）；Ta为日平均气温，℃；aL为扩散短波放射量常数（在平均气候条件下为0.25）；bL为直达短波放射量常数（在平均气候条件下为0.5）；n为日照小时数，h；N为可能日照小时数，h。

3. 潜热通量

潜热通量的计算公式为

式中：ℓ为水的潜热，MJ/kg；Ts为地表温度；E为蒸散发（根据Penman-Monteith公式计算）。

4. 地中热通量

地中热通量的计算公式为

式中：cs为土壤热容量，MJ/（m3·℃）；ds为影响土层厚度，m；T1为时段初的地表面温度，℃；T2为时段末的地表面温度，℃；Δt为时段，d。

5. 显热通量

显热通量可根据空气动力学原理计算，其公式为

式中：ρa为空气的密度；CP为空气的定压比热；Ts为地表面温度；Ta为气温；ra为空气动力力学的抵抗。

虽然上式可与式（2.90）及能量平衡方程联合，用迭代法求解H、G和Ts，但计算量大并且结果不稳定。考虑到日平均热通量和其他热通量相比很小，此处用日气温变化近似日地表面温度变化，在求出地中热通量后，根据能量平衡方程求解显热通量：

6. 人工热排出量

在城市地区，工业及生活人工热消耗的排出量（Ae）对地表面能量平衡有一定影响，根据城市土地利用与能量消耗的统计数据加以考虑。

2.2.4　社会水循环过程模拟

WEP-L模型中的社会水循环模块受外部输入的社会取用水量驱动，需要知道每个计算单元每天的社会取用水量。然而，通常可获取的社会用水数据往往是某个省市的年统计数据。因此，社会用水数据的时间、空间展布是进行社会水循环过程模拟的基础。

2.2.4.1　社会用水数据展布

社会用水按用途可分为农业灌溉用水、居民生活用水、工业生产用水以及生态用水，其中农业灌溉用水和工业生产用水占据比重较大。按来源可分为地表水、地下水、跨流域调水。按时间可分为历史用水和预测用水，历史数据来自统计年鉴，预测数据通过其他方法估算。社会用水根据行政区域数据展布到计算单元上，供社会水循环模拟使用。

农业灌溉用水历史数据主要来自统计年鉴，一般为各省或各地市年数据，需要降尺度展布到计算单元上。根据统计年鉴，灌概用水类型主要有水田、水浇地、林果地、草场、菜田以及鱼塘补水。WEP-L模型将6类用水重新归类成4类：水田用水、水浇地用水（含菜田）、林草用水和鱼塘补水，各类农业用水又包含地下、地表两个部分，即农业用水数据共8个类型。其中水田用水、水浇地用水属于灌溉农田域，林草用水属于植被裸地域，鱼塘补水属于水域。水田、水浇地用水以灌概面积作为依据进行展布：首先，根据年灌溉面积以及土地利用分布图估算各计算单元灌溉面积；其次，根据年降水量估算单位面积灌概定额以及整个计算单元灌溉定额；第三，将水田、水浇地地表用水按灌溉定额加权平均分配到计算单元内。其他类型农业用水则按计算单元面积加权平均分配。时间尺度上，根据当地作物种植结构及灌概制度，将年农业用水分配到旬，从而反映年内灌溉不均匀性，而旬内则按日进行平均分配。

工业生活用水历史数据也来自统计年鉴。空间尺度上，农业生活用水基于农村面积进行面积加权平均分配；工业及城市生活则基于城市面积进行加权平均。时间尺度上，则年内平均分配。

2.2.4.2　社会水循环过程模拟

社会水循环是以“取水-输水-用水-耗水-污水回用-排水”为环节的独立于自然水循环以外的侧支循环过程，水循环通量是受社会经济因素驱动影响的。社会水循环过程是从自然水循环过程中独立出来的不同的循环过程，分析各环节同自然水循环的关系发现，只有“取水、排水”两个环节直接同自然水循环相关，其他环节属于社会水循环内部环节，同自然水体之间无交互关系。因此，在模型拟合“自然-社会”二元水循环过程中，将社会水循环划分成3个模块分别进行模拟：取水模块、用水模块以及排水模块。取水模块主要计算自然水循环不同水体在取水过程中的减少量；用水模块包括输水、耗水以及污水回用等环节；排水模块则将排泄的污水作为一个输入项添加到自然水循环过程中，参与后期自然水循环。

（1）取水过程模拟。取水模块中认为地下水用水直接从计算单元获取，直接引发的效果是降低地下水位；而地表水用水则来自河道或水库，主要导致河道、水库水量的减少。由于WEP-L模型以等高带为计算单元进行水量平衡计算，还需考虑本地河道取水和异地河道取水，即是否从计算单元所在子流域河道取水，模型构建过程中，需要对灌区及计算单元分别进行地表灌溉水来源设置，即指定取水子流域编号及水库编号，优先从水库取水，如果没有水库（或水库没有修建）则从其他子流域取水，如果前两者都没有则从当地子流域取水。其中对灌区的设置主要用于农业灌溉地表用水过程计算；而对计算单元的设置主要用于工业、生活地表用水过程计算。

（2）用水过程模拟。用水模块使用损失消耗系数进行社会水循环用水分量的估算，包括管道输水损失（含管道蒸发及下渗）、工业、生活耗水率（主要指用水过程中蒸发掉的部分，含产品带走的部分；分农村生活、城镇生活、工业生产3个耗水系数）。首先，使用管道损失系数，将社会用水量划分成损失量和净用水量（或农业净灌溉量），其中损失量又进一步划分为下渗量和蒸发量。其次，农业净灌溉量作为额外用水加入到计算单元的水循环过程；工业生活净用水量则使用耗水率划分成蒸发量和污水排放量。

（3）排水过程模拟。排水模块主要功能是按污水来源将用水模块计算所得污水排放量返还到自然水循环过程。农村排放污水直接参与当地坡面汇流过程，即认为农村排水是计算单元坡面汇流的来源之一。由于城市污水通过管网输送到污水处理厂处理后再进行排放，模拟过程中将工业、城镇生活污水使用污水回用率划分成污水回用部分和污水排泄部分（指处理后的）。由于城市污水主要经城市管网进行排泄，模拟过程中认为污水排泄部分直接进入河道。污水回用部分主要作为工业用水、生态用水的一个来源，以模拟城市化过程中的污水回用。

2.2.5　WEP-L模型的构建和模拟步骤

2.2.5.1　模型构建

WEP-L模型采用子流域套等高带作为计算单元，在模型运行前需要获取所有计算单元相关参数数值。基于水文气象、自然地理及社会经济等各类基础数据的收集整理，可以构建研究区域的WEP-L模型，模型具体的构建步骤如图2.4所示。

（1）模拟河网提取。通过DEM提取流域模拟河网，以便后期子流域划分及相关地形信息统计。由于提取的模拟河网往往同实际河网相差很大，一般都需要使用实际河网对DEM进行修正，增加实际河网所在栅格的汇水能力，从而提取出同实际河网相似度高的模拟河网。对DEM进行修正后，使用坡面流累积算法（O'Callaghan和Mark，1984）提取模拟河网，使用D8算法计算栅格流向，使用流向计算栅格汇流累积数，设置河网阈值提取河网。

（2）计算单元划分。使用已提取的模拟河网进行流域计算单元的划分，包括子流域划分、子流域编码和等高带划分。其中，子流域划分是根据模拟河网将流域划分成一系列相互独立的子流域以反映区域的空间异质性；子流域编码则将划分的子流域进行编码标识以反映子流域之间的空间拓扑关系；等高带划分则将各子流域内山区按高程分为一系列等高带以反映高程的影响。

（3）空间属性统计。“子流域属性”统计过程主要计算子流域级别相关信息，如子流域面积、子流域拓扑关系、河道坡降、河道曼宁系数等参数。“等高带属性”主要计算等高带级别相关信息（即计算单元相关信息），如等高带面积、坡度、所属行政区等。“土地利用属性”主要指统计计算单元内不同类型土地利用面积百分比，用以划分5大类下垫面面积，反映下垫面的空间变异性。“土壤类型属性”主要指统计计算单元内土壤类型，该过程采用计算单元内最大面积的土壤类型作为整个计算单元土壤类型。具体土壤类型对应参数则从土壤基础参数文件中读取。“气象要素展布”过程主要使用ARDS算法将气象站数据展布到各子流域形心，作为整个计算单元气象数据输入。“农业用水展布”过程和“工业生活用水展布”过程根据具体展布规则将社会用水展布到计算单元。“社会用水参数设定”过程则主要设定各计算单元相关社会水循环参数，如工业生活取水子流域、输水损失系数、耗水系数等。

（4）基本参数设置。该过程主要设置相关非空间类参数。“土壤基础参数”表示土壤相关的基本参数信息，如饱和导水率、田间持水率、饱和持水率等。WEP-L模型中将土壤类型概化成4类：砂土类、壤土类、黏壤土类以及黏土类。“输入输出控制文件”用于记录模型所使用的相关输入输出文件信息。“基本控制参数”主要指模拟时间、是否考虑用水、是否考虑水库等控制参数。“情景设置参数”主要指用以描述特定情景的参数。“自动调参参数设置”主要用以程序自动调参相关功能。“水库信息”指水库月蓄变量资料以及水库基本信息（如兴利库容、修建年份等），主要用于水库调节及水库取用水。“灌区信息”指灌区相关取水子流域和取水水库编号信息。

经过上述步骤可以完成区域WEP-L模型构建，之后根据具体问题进行模型模拟和结果分析。

2.2.5.2　模型模拟

WEP-L模型可以单次顺序模拟多个情景，每个情景按逐年、逐月、逐日循环模拟，进而实现长系列模拟，其循环模拟步骤如图2.5所示（刘佳嘉，2013）。

各情景模拟开始时读取相关输入数据并进行参数初始化，如子流域信息、等高带信息、水库信息等，这些参数信息最主要的特征是在年月日循环中不会被改变。可通过设置自动调参开关，决定是否进行自动调参模拟。模型支持GLUE调参方式。设置好相关自动调试参数后，进入模型年月日循环体进行具体模拟。如果不自动调试，则直接进入年月日循环体。

WEP-L模型以日为基本单元进行模拟，分年、月、日三个级别进行嵌套循环，每个级别循环开始都进行相关初始化，主要工作是相关功能参数初始化以及和年月日相关的数据读取，如土地利用类型数据的设置以及社会年用水的读取均在年循环初始化中实现；而计算单元日气象数据则在日循环初始化中完成。后续扩展中，只需知道参数变量的时间尺度，只要在对应尺度循环初始化模块进行设置即可。各尺度循环模拟结束后，按要求输出逐日、逐月、逐年相关参数，用以结果展示。

逐日循环中的“循环主体”模块是WEP-L模型核心所在，用于实现“自然-社会”二元水循环的模拟。概括而言，循环主体可以分为4个子模块：社会水循环模块、产流模块、土壤水-地下水模拟模块、汇流模块。

（1）社会水循环模块。社会水循环模块包括两个部分：社会用水量的计算和河道、水库取水量计算。其中社会用水量的计算指使用输入数据进行相关分量的累加，获取水域、植被-裸地域、灌溉农田域相关地表水、地下水使用量。如果存在跨流域调水，则相应减少受水区域地表、地下水开采量。根据统计所得地表水使用量分别从本地河道、异地河道以及水库中减去相应数值，模拟社会水循环的地表取水过程。地下水取水量则在地下水模块使用。

（2）产流模块。模型产流模块分暴雨期和非暴雨期分别进行，其中暴雨期采用逐小时模拟，而非暴雨期采用逐日模拟。分5类下垫面分别进行产流模拟，通过面积加权平均累加到整个计算单元。

（3）土壤水-地下水模拟模块。模型土壤水、地下水模块主要进行土壤水、地下水相关垂直运动、水平运动模拟以及地下水位变化情况。该模块模拟地下水补给、土壤水动态变化、壤中流出流、地下径流等一系列土壤中相关水分运移，同时也模拟土壤蒸散发过程。对地下水而言，通过输入输出水量平衡计算日尺度地下水位变化。

（4）汇流模块。模型汇流过程主要包括坡面汇流过程和河道汇流过程。汇流过程以6h为一个时段进行模拟，采用“牛顿下山法”进行运动波方程求解。河道汇流计算输入项包括坡面产流和上游来水，输出项主要是水域蒸发（其中人工河道取水量在模拟开始社会水循环模块进行了消减）。根据方程计算出河道流量后，还需根据是否有水库进行水库调蓄计算。水库调蓄使用月蓄变资料进行调蓄，根据逐日计算流量累加计算，使得全月蓄变量满足历史资料。如果没有月蓄变量资料，则使用最小下泄流量进行调蓄，即如果计算流量小于最小下泄量，则实际计算流量等于最小下泄量，否则使用计算值。