第3章 基于 SWAT的海伦市地表水径流模拟

        SWAT (Soil and Water assessment Tool) 模型是反映流域下垫面空间分布特征和地表径流模拟的半分式水文模型。

通过将 ArcSWAT、AvSWAT 模块嵌入 ArcGIS 软件，实现了地表径流模型的可视化效果。相比传统水文模型，其空间数据管理、数据处理和展示功能大大增强 [75]。模型建立所需的基础资料包括土壤特征、流域的地形、气象资料、土地利用、水文资料等 [76]。

SWAT 模型与 ArcGIS 紧密结合，在建模前需对研究的基础数据进行整理与分析，在此基础上构建模型需要的用户数据库，例如土壤、土地利用类型、天气发生器等。

本文以通肯河流域上游为研究区，联合利用 ArcGIS 和 SWAT 软件，建立海伦市基础数据库和 SWAT 数据库。

3.1模拟理论与运算过程

SWAT 模型由许多变量（蒸散发量、河道汇流量、地表径流量、土壤水分运移、地下水补给）组成，实现对流域内的水文循环的计算与模拟 [78]。

受流域下垫面性质和降雨过程的影响，将 SWAT 模型分为子流域计算模块和河道产汇流模块 [79]。

子流域计算模块对流域河道内水的输入量进行模拟，其受气候、地表产流和地表覆盖等因素的影响；河道汇流模块计算水、泥沙、污染物等物质从个子流域向流域总出口的输送情况[80]。

水量均衡是模型计算和统计流域水分随时间变化的根本依据，SWAT 模型内置的水量平衡均衡模块，运算式如下[81-84]：

其中：SW0为初始土壤水量，单位mm；SWt为t时间节点土壤含水量，单位mm；t为时间，单位d；Rday为第i天的降水量，单位mm；Qsurf为地表汇水量，单位mm；E为蒸发量，单位mm；Wseep为土壤层的旁侧流量，单位mm；Qgw为地下水的回归流量，单位mm。

3.1.1地表径流

在SWAT模型中，通过将所有HRU的计算结果汇聚到子流域，得到河道出水口的径流量。地表径流的多少受土壤水分的影响较大。如果土壤处于干燥状态，此时下渗率最大，地表径流不会形成；随着土壤逐渐湿润且水分的增多，下渗率逐渐减小，逐渐达到填洼临界点，即地表径流产生[85]。对于SWAT模型，模型自带地表径流估算方法：径流曲线法(SCS)和格林-安姆普特(GREEN-AMPT)。本次研究获取到日尺度数据，因此使用SCS径流曲线法来进行地表径流计算。计算公式如下[86-90]：

 Q_{surf} = \frac{(R_{day}-I_a)^2}{(R_{day}-I_a+S)}

S = 25.4 \times (1000 / CN - 10)

其中: Qsurf—地表汇水量,单位mm;            Rday—当天降水量,单位mm;

Ia—降雨初损,单位mm; S—潜在最大滞留量,mm;            CN—曲线系数,流域降雨前期特征的参数,无量纲。 

3.1.2蒸散发量

        SWAT模型中的蒸散发包括:植物蒸腾与升华、树冠截留的水分蒸发、土壤水的蒸发。流域水量输出的主要途径是蒸散发。蒸散发计算:Priestley-Taylor、Hargreaves、Penman-Monteith方法。本次研究选用Penman-Monteith方法[91-93],详细计算公式如下: 

\lambda E_t = \frac{\Delta(H_{net} - G) + \gamma k_1 (0.622 \frac{\lambda \rho_{air}}{P}) e^0_s - e_s/\gamma_a}{\Delta +\lambda(1+\gamma_c/\gamma_a)}

其中:SW0—热通量密度,MJ/m2·d;

E—潜发深率,单位mm/d;

Δ—曲线斜率,kpa/℃;

Hnet—净辐射,单位MJ/m2·d;

G—热通量密度,MJ/m2·d;

ρair—空气密度,单位kg/m3;

Cp—热量,MJ/kg·℃;

P—大气压,单位KPa;

k1—8.64×104;

es—大气水汽压,单位KPa;

es0—饱和水蒸汽压;

γc—植物叶孔阻抗;

γa—空气动力学阻抗。 

3.1.3 土壤水分运移

        入渗到土壤中的水存在多种运动形式，有一部分通过蒸散发损失掉，其中一部分通过入渗补给进入地下水含水层，余下为壤中流。SWAT利用动力贮水方法计算 [94-96]：

Q_{lat} = 0.024 \times \left( \frac{2 \times SW_{ly,excess} \times K_{sat} \times SLP}{\Phi_d \times L_{hill}} \right)

其中:

Qlat—初始土壤水量,单位mm;

SWly,excess—流出的水,mm;

Rday—第i天的降水量,单位mm;

Ksat—水力传导率,mm/h;

SLP—坡度,无量纲;

φd—土壤层的孔隙度,无量纲;

Lhill—坡长,单位m。

3.1.4 地下水

        对于地下水，SWAT中分为浅层和深层。浅层是指产生的基流供产流；深层指流域以外的水 [97-98]。根据水量平衡方程，浅层地下水的计算如下：

aq_{sh,i} = aq_{sh,i-1} + w_{rchrg} - Q_{gw} - w_{revap} - w_{pumpsh}

其中:aqsh,i—i天蓄水量,mm; aqsh,i-1—i-1天蓄水量,mm; wrchrg—浅层地下水的补给量,单位mm; Qgw—基流量,单位mm; wrevap—第i天土壤带的水量,单位mm; Wpump,sh—提取的水量,mm。

 深层地下水的水量平衡公式为:

Q_{gw} = \frac{8000 \times K_{sat} \times h_{wtbl}}{(L_{gw})^2}

其中:hwtbl—浅层地下水的水位,单位m; Ksat—地下水的水力传导率,单位mm; Lgw—至主河道距离,m。 

3.1.5 河道汇流

        根据产汇流及水流方向，将子流域可细化为支流与干流。支流的水主要来源是地表汇流，汇流时间取决于河道长度和地形坡度。干流水量来源于支流和地下水的补给。为了进行河道水流计算，SWAT模型建立变动存储系数模型，计算公式如下[99-100]：

q_{\text{out},2} = SC \times q_{\text{in},\text{ave}} + (1 - SC) \times q_{\text{out},1}

其中:qout,1—初始时间的水流速率,单位m3/s; qout,2—结束时间的水流速率,单位m3/s; qin,ave—平均引流速率; SC—储存系数,单位mm。 

3.2 数据库构建

        SWAT水文模拟需要构建空间数据和属性数据库，其中主要包含土地利用、数字高程系统DEM、土壤属性数据、水文数据、土壤空间数据、气象数据等，具体数据来源、格式与精度见表3.1。由于SWAT模型需要对不同空间数据叠加，所以要求所有空间数据拥有统一坐标系。投影坐标系选定为GCS_Xian_1980。

表3.1 空间数据和属性数据情况表 

数  据	分辨率或比例尺	格  式	来  源
DEM	90m×90m	GRID	地理空间数据云
土地利用	1:25万	mapgis	OLI和GF-1遥感数据
土壤类型	23m×23m	ArcGIS	世界土壤数据库
气象数据	2008年-2016年	DBF	寒区旱区科学数据中心
水文数据	月尺度	XLSX	水文测站

3.2.1 DEM高程数据

        DEM全称是数字高程模型,用地形高程来模拟不规则的地形曲面,从而建立数字化的地形模型[101]。在SWAT模型中,利用DEM数据来实现流域空间离散,随后对子流域进行划分、生成水文响应单元(HRU),从而得到流域的边界、子流域、河网分布图。此外,通过DEM还可以得到河道长度、坡度、坡向、流域面积、平均高程、最小高程、最大高程及剖面曲率等地形参数[102]。在本次研宄中,利用ArcGIS进行对DEM数据进行拼接,获得研究区完整的DEM图。之后,对DEM进行投影变换,即从WGS84转为GCS_Xian_1980。最后,利用ArcGIS工具箱中的图形裁剪工具,得到海伦市的数字高程图,如图3.2所示。

3.2.2土地利用类型数据

        根据解译结果,2019年海伦市土地类型包括:水田、林地、旱田、草地、湖泊、河流、沼泽、水库、沙地等。土地利用分类则以旱地为主,占地面积2968.95km²,占总面积的63.52%;其次为水田,面积为783.07km²,占总面积的16.74%,分布在扎音河、海伦河、通肯河漫滩及阶地地区;林地,占地面积460.10km²,占总面积的9.85%,分布海伦东北部的丘陵地区;建设用地,占地面积191.56km²,占总面积的4.10%,分布海伦市区及各乡镇;草地,占地面积181.94km²,占总面积的3.89%,分布东方红水库下游、通肯河及海伦河下游地区;水库,占地面积为44.70km²,占总面积的0.94%;河流及沼泽,占地面积为32.60km²,占总面积的0.68%;沙地占地面积最少,仅占地面积4.47km²,占总面积的0.09%(图3.3)。 

SWAT模型中的土地的利用类型,由美国地质调查局(USGS)定义,需对2019年遥感解译结果按照USGS标准进行重分类,具体分类标准及SWAT查找表见表。将重分类结果输入SWAT模型,得到研究区SWAT模型中土地利用类型分类[103],图3.4所示。 

3.2.3土壤类型数据

土壤数据是SWAT模型最基本的参数,土壤参数的质量影响模拟的精度。土壤参数:物理性质、空间分布、索引关系表[104]。本次研究使用的联合国粮农组织与维也纳国际应用系统研究所联合建立,名为世界和谐土壤数据库(HWSD)。土壤空间分辨率1km,分类系统采用FAO-90。  (1)土壤属性库 土壤属性库包括:土壤物理属性和化学属性。 土壤的物理属性反映土壤侵蚀的程度及规律,其主要包括土壤名称、饱和导水率、土壤分层数、土壤水文学分组、有效持水量、湿密度、土壤孔隙度、层结构和土壤有机碳含量等参数。化学属性反映氮、磷含量的初始浓度,是非点源污染研究中重要参数之一。由于本研究模拟的是通肯河流域(海伦)的径流过程,因此只考虑土壤的物理属性,不考虑其化学属性。模型需要输入的参数如表3.3所示。 

(2)土壤水文分组 根据土壤的渗透特征,美国自然资源保护局(NRCS)将HYDGRP分成四类(A、B、C、D)[105]。如表223.4所示,土壤的水文性质如下。 

(3)SPAW参数 土层有效含水量(mm)、土壤湿容重(g/m3)、饱和水力传导系数(mm/hr)和土层结构参数都是通过SPAW软件计算。如图3.5所示,在输入Sand、Clay、Organic Matter、Salinity、Gravel、Compaction等参数值后,就可输出计算结果。

(4)土壤侵蚀力因子 土壤数据库中的土壤侵蚀力因子[106],可使用公式(9)-(13)计算:

 K_{USLE} = f_{csand} \cdot f_{cl-si} \cdot f_{orgc} \cdot f_{hisand} 

 f_{\text{csand}} = 0.2 + 0.3 \times e^{-0.25 \times sd \times (1 - \frac{st}{100})}

f_{\text{csand}} = 0.2 + 0.3 \times e^{-0.25 \times sd \times (1 - \frac{st}{100})}

f_{\text{orgc}} = 1 - \frac{0.25 \cdot c}{c + e^{(3.72 - 2.95 \cdot c)}} 

其中:fcsand—粗砂土的侵蚀参数,无量纲; fcl-si—黏壤土的侵蚀参数,无量纲; forgc—土壤有机质参数,无量纲; fhisand—高沙土的侵蚀参数,无量纲; sd—沙粒的百分比,无量纲; si-粉粒的百分比,无量纲;cl—黏粒的百分比,无量纲; c—有机碳的百分比,无量纲。

通过计算,流域土壤的物理属性参数均已获取。将整理好的土壤物理属性数据库与SWAT模型中的user soil子模块参数对应,构建研究区土壤数据库,如图3.6。 

3.2.4气象资料

气象数据输入:降水量(逐日)、气温(逐日、范围)、观测站逐日太阳辐射、相对湿度(观测站)和子流域平均风速,以上参数均以天尺度的实测值 [107]。太阳辐射数据不是直接获取数据,只能利用天气发生器模拟产生,其余数据来源为中国大气同化数据集(China Meteorological Assimilation Driving Datasets for the SWAT model,CMADS Version 1.0)所提供数据,观测站基本情况见表。选取的时间尺度为2008年-2016年。数据格式为数据集dbf和文本txt。时间分辨率:日(day);空间分辨率:1/3°。

3.3模型建立与运行

土壤数据库构建起来之后,SWAT模型还需进行子流域的划分和流域HRU生成。 

3.3.1子流域

本次研究,子流域的划分主要基于DEM数字高程数据,运用SWAT模块中的 Watershed Delineat子程序自动划分。在划分过程中,首先,导入DEM和已知河网的矢量(.shp)文件;其次,调整子流域面积的阈值,确保DEM生成的河网与导入的河网基本一致;然后,根据研究区实际情况,人工添加子流域的出口(42#、37#),水文测站(35#)及水库信息(10#、33#);最后,点击参数计算,生成流域分区图。 流域划分情况如图3.7所示,子流域总数为42个。其中,海伦河流域有29#、30#、31#、32#、33#、34#、36#共7个子流域;扎音河流域有9#、10#、13#、14#、17#、18#、19#、22#共8个子流域;克音河流域有15#、16#、23#、24#、37#共5个子流域;三道乌龙沟流域仅有41#子流域,其余均为21个均为通肯河干流及其子流域,各子流域参数见表3.6。 

3.3.2水文响应单元

水文响应单元(HRU)在SWAT模型中,其意义是参与水文计算与模拟的最小单元,同时也是子流域的重要组成部分。每一个HRU计算单元包含了单元内的坡度、地形、土壤类型、土地利用类型等影响水量计算结果的重要信息[108]。利用Land/ Soils/Slope Definition子模块,在输入土壤分类、土地利用类型、地形坡度等参数初始值后,整个流域的子流域的HRU就被成功定义,详见图3.8。

3.3.3模型运行

在将气象数据和水文响应单元基本数据导入模型后,再选择径流计算方法、模拟方法、潜在蒸散发估算等方法的选择与输入。径流值以实测月平均值为单位,输出结果亦选择月尺度;选择模型默认的偏正态分布方法对降雨量进行模拟与计算;选择Muskingun方法对河道径流进行分析与计算[109];利用SCS径流曲线数法对河道径流进行计算;Penman-Monteith方法用以估算与模拟潜在蒸散发;SWAT模型运行结果选择输出:子流域(.sub)、水文响应单元(.hru)、河道(.rch)。

3.3.4敏感性分析

选取与径流密切相关的25个参数进行分析[110],这些参数分为三类:土壤特性参数、气温参数、基流特性参数和河道参数。 土壤特性4个参数:有效水容量(SOL_Z)、生物混合效率(SOL_AWC)、饱和水力传导度(SOL_K)、土壤湿润反照率(SOL_ALB)。SOL-K代表土壤的入渗速率,随着该值的增大,土壤和含水层得到入渗补给逐渐增多,不利于产流;土壤的蓄水能力由SOL-AWC表征,其值越大,表明蓄水能力越强,而产流会越小。 气温参数:融雪基温(SMTMP)、积雪温度滞后(TIMP)、融雪因子(SMFMX:6.21)、降雪温度(SFTMP)和融雪因子(SMFMN:12.21)。融雪基温、降雪温度、融雪因子和积雪温度滞后系数都会影响径流的大小和径流滞后降雨的时间。 基流特征参数:“基流”阈值(GWQMN)、延迟时间(GW_DELAY)、再蒸发系数(GW_REVAP)、alpha因子(ALPHA_BF)、再蒸发系数(REVAPMN)。在浅层含水层中,蓄水量超过GWQMN阈值,则地下水向子流域河道进行排泄;浅层含水层的临界持水量用ALPHA-BF表示,当含水层中蓄水量超过ALPHA-BF阈值时,地下水再蒸发发生;再蒸发过程中,GW-REVAP决定含水层损失的水量;在浅层含水层,水分进入有一定的延迟,该时间GW-DELAY表征。 河道参数:水力传导度(CH_K2)、SCS曲线(BIOMIX)、径流曲线系数(CN2)、曼宁系数(CH_N2)。CN2代表下垫面的产流能力,该值大小影响径流量的多少。对于降雨量一定, CN2值越大,地表径流增大,产流也会增加[111]。 在SWAT-CUP软件中,利用SUF-2算法,对以上25个参数进行参数率定。为了寻求最优的参数值,选代次数为2000次,最优解的范围不断缩小,最终得到最优解,率定结果见表3.7。

对25个率定期的参数敏感性分析,见表3.7,ALPHA-BF、CH_K2、CANMX、REVAPMN、SOL-K和TIMP是影响水文过程中的重要因子,因此这几个参数对流域径流敏感性较强。 3.3.5模型适用性分析 通常情况下,纳什系数(ENS)和相关系数(R2)评判所构建的模型的适用性[112]。 (1)相关系数(R2) R2是实测径流值和模拟值之间的相关性系数,利用Excel按线性回归法求得。R2→1时,越接近1,预测值与实测值越吻合,模型适用性越强;R2<1时,其值越大,两者的相似度就越高。 

第4章 地表水-地下水耦合模型 地球上水的循环过程,是多环节同时参与的过程,涉及到热量传递与转换、水量(地表水、地下水)蒸发、水汽输送、地表水循环、地下水循环以及其他形式的能量、物质的存储与转换。其中,水循环过程中降水、径流、蒸发是最重要的环节,这三个环节决定了流域内的水资源总量,是整个流域实现水量平衡的重要途径。流域内水循环包括:降水→地表水、降水→地下水、地表水与地下水相互转化、地表水→排泄、地下水→排泄等过程。根据水循环的机理,为了使计算结果更为符合流域实际情况,应将地表水径流和地下水径流过程有机结合[113]。本章根据ArcSWAT计算结果,调用各子流域的输出项,作为Visual MODFLOW Flex对应单元的输入项,进行流域地下水的数值模拟,从而实现了海伦地区地表水和地下水的耦合模型。

（3）SWAT水文过程数值模拟系统收集土地利用资料、地形等高线和高程点、气象资料、土壤类型资料、渠系和灌区分布、水系和水文站点的径流资料、蒸散发量等数据个图件，并利用遥感技术对相关数据进行补充。利用ArcGIS软，基于高程点和等高线等高程数据，结合DEM数据集，构建研究区数字高程模型，提取子流域所需的地形参数，生成研究区内河网，划分子流域（水位响应单元）。基于研究区倍土壤利用类型实际调查数据，结合土地利用类型现状图（中科院资源环境科学与环境数据中心），并与SWAT土地利用数据库中的分类进行对照，划分土地利用类型。利用区内土壤调查及其物理参数，结合HWSD数据库，结合颗粒分析和非饱和水力参数实验结果，对土壤水分特性进行分组，形成ArcSWAT模型中的土壤数据库。利用ArcSWAT软件构建研究区水文模型，利用实测水文数据（径流量、地下水位波动等）对模型参数进行敏感性分析和参数识别；并与均衡法计算的入渗补给、蒸散发排泄、地下水-地表水相互作用进行验证对比，修正模型参数，更为准确地确定地下水系统的入渗补给和排泄。

在SWAT模型中，利用DEM数据来实现流域空间离散，随后对子流域进行划分、生成水文响应单元(HRU)，从而得到流域的边界、子流域、河网分布图。