水文模拟一直是水科学研究的热点和重点, 也是水资源研究的基础^[1], 通过水文模拟可以了解各地区水资源状况, 揭示水文规律, 为水资源开发利用和进一步的水文研究奠定基础。国内外大量的研究在流域尺度上对水文过程进行了模拟和水资源状况进行了分析^[2-4]。中国的各大流域的水文过程得到广泛模拟, 其中SWAT水文模型也成为水文模拟的重要工具, 在黄河、长江、淮河、渭河等流域取得了较好的模拟效果^[5-9]。

水是西北干旱区生命和经济活动之源, 有水就是绿洲, 无水便成沙漠。水资源短缺是制约中国干旱地区社会经济与生态环境可持续发展的主要因素^[10]。不少研究者对中国干旱地区的流域进行了水文模拟^[11-13], 例如宋一凡等^[11]将SWAT模型应用到艾布盖河流域进行水文过程的模拟, 丁程峰等^[12]利用SWAT模型模拟了乌鲁木齐河流域对森林变化的水文响应过程。这些研究对掌握西北干旱区对应流域水文水资源的形成、演化和分布规律, 理解气候变化和人类活动对流域水文过程的影响做出了重要贡献。然而干旱区由于水资源更加短缺以及对经济社会发展的制约作用, 国家对水资源的管理和配置不断加强, 例如黑河流域水资源和水权的分配。但是这种管理和配置是基于县域尺度, 以往流域尺度的水文模拟无法满足实践的需求, 必须在流域尺度上进一步细化。

总之, 虽然流域尺度的水文模拟相关研究相当多, 取得了重要进展, 但是流域尺度上的水文模拟需要进一步细化到县域尺度上才能更好地满足水资源管理和分配的需求。因为流域水资源管理和分配往往是以行政单位进行, 决策者迫切需要掌握县域范围内的水文状况。然而县域尺度的水文模拟研究鲜有报道, 县域水文过程又是怎样及其对水资源产生何种影响？这些问题并不清楚。要回答这些问题需要在县域尺度上进行水文模拟, 这对流域水资源管理和配置具有重要的科学指导意义。

基于此, 本文以西北干旱区的临泽县为典型研究区, 在县域尺度上利用SWAT模型进行水文过程的模拟, 揭示1980—2012年临泽县水文要素变化规律并结合Mann-Kendall方法分析其年际变化, 进而估算临泽县1980—2012地表水资源和生态用水的数量并阐明其动态变化特征。

1 数据和方法 1.1 研究区概况

临泽县(38°57′—39°42′N, 99°51′—100°30′E)位于甘肃省河西走廊中段, 地处巴丹吉林沙漠南缘, 隶属于张掖市, 是河西走廊上的一个重要农业县。县域面积2727 km², 东西宽约49.7 km, 南北长约77 km；地势南北高, 中间低, 由东南向西北逐渐倾斜, 海拔1380—2278 m(图 1)。研究区年均气温7.6℃；年均降水量117 mm, 多集中在7—9月, 约占全年降水量的65%, 年均日照时数可达3045 h；土地利用类型有耕地、林地、草地、水域、沙漠、戈壁、盐碱地、沼泽等^[14]。地表水资源主要来自黑河、梨园河。全县为灌溉农业, 有平川、板桥、鸭暖、蓼泉、小屯、梨园河(新华、沙河、倪家营)6个灌区。农作物有小麦、制种玉米、番茄、棉花、温室蔬菜等, 沙地植被以灌木为主, 有沙拐枣(Calligonum arborescens Litv)、梭梭(Chenopodiaceae)、柽柳(Tamarix chinensis Lour)、花棒(Hedysarum scoparium)、柠条(Caragana intermedia intermedia)等, 乔木有二白杨(Populus gansuensis C. Wang et H. L. Yang)等。土壤类型包括灰棕漠土、风沙土、灰钙土、灌耕土、潮土、盐土、草甸土和沼泽土8个类型。

1.2 SWAT模型

SWAT(Soil and Water Assent Tool)模型是美国农业部(USDA)农业研究中心(ARS)1993年为美国水文模型(HUMUS)项目开发的大、中尺度流域环境模型, 已经在美国18个主要流域以及其他国家得到广泛应用。其水量平衡公式如下^{[3, 4]}:

(1)

式中, ΔSW、P、Q、Ea、DP和GGW分别为土壤含水量、降水量、地表径流量、实际蒸散发量、深层下渗量和地下水补给量。DP可分为产流的下渗(PERC)和汇流的河道下渗(TLOSS), 浅层回归流包括壤中流(LATQ)和地下径流(GW), 地表径流采用改进的SCS模型的径流曲线数方法计算。

本研究利用决定系数(R²)和Nash\Sutcliffe效率系数(NS)来衡量模型校准期和验证期的模拟效果。决定系数衡量模拟值和实测值的拟合效果, 范围为0—1, 值越大表示模拟效果越好。一般来说R²大于0.5模拟效果就可以接受。纳西系数(NS)的取值范围为0—1, 当NS=1, 说明模拟值与实测值完美匹配。一般来说, 如果NS≥0.75, 模拟效果很好；如果0.36＜NS＜0.75, 模拟能够满足应用；NS≤0.36则模拟达不到要求^[3-7]。

(2)

(3)

式中, Q_sim为模拟值；meanQ_sim为模拟值的算术平均值；Q_abs为实测值；meanQ_abs为实测值的算术平均数；t指值的个数。

1.3 Mann-Kendall突变检验

序列观测值(Y₁, Y₂, …, Y_n), 用Y_q表示第q个样本Y_q＞Y_p的累积数(其中1≤p≤q)。

统计量：

(4)

式中, 如果Y_q＞Y_p, 则L_q的值为1；反之则为0。将统计量标准化为：

(5)

式中, E(F_j)为F_j的平均值, E(F_j)=j(j-1)/4；Var(F_j)为方差, Var(F_j)=[j(j-1)(2j+5)]/72。

在给定α显著水平情况下查正态分布表确定U_α/2, 如果UF_j＞U_α/2或者UF_j＜U_α/2, 则表明原时间序列具有很明显的增强或减弱趋势。将其所有值绘制一条曲线UF。将时间序列Y按逆序列排列。按公式4重新计算得到UB_j, 绘制第二条曲线UB, 并使两次计算结果满足

(6)

式中, 统计量序列曲线UF、UB和±U_a/2均绘在同一坐标平面内分析曲线走向。当曲线UF大于0则表明序列呈上升趋势, 反之则表明呈下降趋势。当超过信度线值则表明上升或下降趋势明显。如果两条曲线有交点且交点在信度线之间, 则表明这个点所对应的时刻便是突变点^[15]。

1.4 基础数据处理

从中国科学院环境数据中心下载1986、1995、2005的3期土地利用数据, 在ENVI软件平台中完成剪切、镶嵌和投影转换等数据预处理。1995年的土地利用数据用于模型的初始构建。模型需要的土地利用数据包括空间图和对应的索引表。本研究利用如下索引表将土地利用类型转换为SWAT模型能识别的4位代码(表 1)。

SWAT把土壤分为2层进行模拟, 土壤数据利用Harmonized World Soil Database(HWSD)的数据, HWSD数据土壤粒径分类方法和美国制的相差很小。土壤的物理属性可以通过HWSD的数据库中查到, 土层的湿容重、有效含水量和饱和渗透系数3个土壤属性数值直接通过SPAW软件计算获得。土壤侵蚀力因子(USLE_K)的计算可以通过美国通用方程(USLE)计算获得。

数字高程数据(DEM)从中国地理数据云(http://www.gscloud.cn/)下载, 分辨率为30 m。按照行列号共下载12景影像, 在ENVI软件中镶嵌后利用数字化边界进行裁剪, 投影转换为Alber投影系统。DEM栅格中存在“凹陷栅格”, 因此需要进行“填洼”的预处理, 研究中利用三次卷积算法对DEM进行“填洼”的预处理。

通过中国气象数据共享中心(http://cdc.cma.gov.cn/)下载7个气象站(札马什克、祁连、托勒、山丹、张掖、金塔、高高崖、临泽)1980—2012年日降雨量、日平均气压、日最高温、日最低温、日平均风速、日照时间长度和日相对湿度数据, 对数据进行整理与插补后构建SWAT模型气象数据库。

利用河道径流数据模型校准和验证, 河道径流数据来自于1980—2012年的甘肃省、青海省和内蒙古自治区的水文年鉴以及黑河流域管理局；选择札马什克、莺落峡和正义峡3个水文站作为验证站点。

2 结果与分析 2.1 模型校准和验证

SWAT模型首先根据DEM划分出子流域(Subbasin), 在此基础上根据土壤、坡度和土地利用数据进一步将子流域划分成水文响应单元(HRU)。研究中所有土地利用和土壤类别都参与HUR的划分。临泽县最终被划分成28个子流域, 98个HUR。河道实测径流数据用于模型的校准和验证, 校准期为1980—1997年, 验证期为1998—2012年。从表 2和图 2中可以看出模拟结果比较理想, 完全能够满足应用的需求。

2.2 年均水量平衡

计算得到临泽县1980—2012年的年均水量平衡。临泽县年均降雨量为117.99 mm, 年均实际蒸发量为130.53 mm, 年均下渗49.78 mm, 年均壤中流补给量(LATQ)为68.95 mm, 年均地下水补给量为48.00 mm, 年均产流量为53.47 mm。临泽县年均地表径流接近0, 年均产水系数为0.16。年均实际蒸发量比年均降水量大12.54 mm, 产流量主要由壤中流和地下径流补给形成, 二者的贡献比大约是3:4。

2.3 年水文分量的空间分布

从图 3可以看出, 临泽县产水系数范围为0—0.92, 南部区域数值比北部大, 最大值出现在北部区域。从整体上看, 临泽县产流系数有北低南高的趋势, 空间差异比较大。将灌区空间分布与产流系数空间分布叠置可以看出临泽县南部地区是灌区, 产流域系数以蓼泉灌区附近为极大值中心向周围递减。我们把这种以每一地区产流高值为中心向周围逐渐降低的现象称为“产流岛效应”。与非灌区相比, 灌区的产流系数要高。岛中心的多年平均降雨量、实际蒸散发、入渗量、地下水补给量(浅层地下水和土壤径流)为2.34×10⁷, 2.48×10⁷, 9.9×10⁷, 2.028×10⁸ m³, 产流量为1.0210⁸ m³。临泽县年均地下水补给量(图 3)的空间差异也比较大。变化范围为0—2.02×10⁷ m³。北部地区地下水补给量很少, 接近0。南部地区逐渐增加, 并在蓼泉灌区附近达极大值(2.02×10⁷ m³)；然后以此为中心向四周逐渐下降。地下水补给量空间分布与产水系数十分相似。这说明临泽县的产流在很大程度上是受地下水补给量的控制。临泽县属于内陆干旱区, 在北部区域, 实际蒸发已经大于降雨量, 地表基本不产流。而在灌溉农业区, 水分过程主要为降水、凝结水入渗、蒸发、人工取水-灌溉等, 以垂向水分运动为主。子流域之间的横向水力联系变弱, 横向水力联系主要在某位置地下水位变动时周边区域的侧向补充。河水、引灌河水、降水、凝结水入渗补给了地下水, 而地下水通过泉水溢出和人工开采等排泄方式补给地表径流, 而大量的地下水补给填补了蒸发量的亏损并产生河川径流^[14]。

临泽县年均实际蒸发量变化量变化范围为0.1×10⁵—565.42×10⁵ m³, 实际蒸发量都比较大。东部地区最大, 极大值出现在板桥灌区和邻近区域(图 3)。西南地区的梨园河灌区实际蒸发量次大。整体上而言, 临泽县的实际蒸发量异质性比较大, 但没有明显的南北分异。

2.4 水文分量的月变化

如图 4所示, 研究区月平均温度状况接近正态分布。最低气温出现在1月, 为-0.91℃, 随后气温持续性的上升, 在7月达到最大值, 为27.32℃；然后持续性的下降。临泽县有2次汛期, 夏讯出现在7月, 春汛出现在3月；7月以后降雨量开始下降, 10月后有所回升。降雨主要集中在5—9月, 占全年降水量的62.29%。临泽县实际蒸发量的月份变化与气温、降雨有大体相似的趋势, 但波动性更大。春季的3月开始实际蒸发量增加变缓, 对应的降雨量出现次高峰后开始下降。7月出现最大蒸发量, 大约7×10⁷ m³左右, 然后持续性的下降；实际蒸发量主要集中在5—9月, 占全年的65.08%。11月到次年3月之间, 降雨量要大于实际蒸发量；在4—10月, 实际蒸发量大于降雨量。壤中流的补给分布呈现倒“U”字形, 1月最低, 然后缓慢上升, 在5—10月基本保持, 11月后开始缓慢下降。研究表明实际蒸发(E_a)主要受气温和降水的影响, 但是降雨是临泽县E_a的限制性因素, 而不是气温所带来的蒸发能力, 因此其分布曲线与降雨曲线更为相似；但同时受到下垫面的干扰, 呈现出波动性。壤中流的补给主要集中在3—10月, 占全年的90%。产水量分布接近一条直线, 全年的变化很小, 每月大约1×10⁷ m³左右。这主要由于产水量主要来源于地下水补给, 比较稳定。

2.5 水量平衡年际变化

将临泽县1980—2012年年降雨量(P)、年地下水补给量(GGW=GW+LATQ)、年实际蒸散发量(E_a)、年下渗量(PECR)、年河道损失量(TLOSS)、年产水量(YILD)的时间序列利用Mann-Kendall法进行计算绘制成图 5。

由UF曲线和UB曲线可知, 1980—2012年间降雨量(图 5)的突变点较多, 分别为1983年, 2006年和2010年。1980—2005年, 降水保持着递增趋势(UF>0), 2006—2010年降水量下降(UF < 0), 21世纪初期开始上升。从图 5中可以看出, 实际蒸散发量(图 5)从1982年开始保持持续上升的趋势, 突变年份出现在1983年。年下渗量(图 5)在1980—1984年呈现出上升趋势, 20世纪80年代后期保持下降趋势, 21世纪以后持续下降, 突变点出现在2000年。河道损失量(图 5)的波动不大, 基本分为2个阶段, 2003年前是上升趋势, 2003年以后则下降。其中1980s和1990s上升趋势比较显著(UF>1.96)。地下水补给量(图 5)在1982年有稍微的下降, 在1980s和1990s保持上升趋势, 2001年以后则开始下降, 突变点出现在2002年。产水量(图 5)是各水文参数的综合反映, 其变化趋势与地下水补给量相似, 在1980—2003年持续上升, 2003年后则持续下降, 突变点出现在2003, 1990s增加趋势十分显著。

利用相同的方法对年均温(T_e)和年潜在蒸发量(PET)进行动态变化检测(图 6)。结果表明气温(图 6)1980—2012年尽管有所波动, 但整体上保持持续上升趋势。1998年后气温上升趋势十分显著, 突变点出现在1990年。潜在蒸发量(图 6)波动性相对气温波动性大一些, 突变点出现在2010年。1980—2007年保持着持续上升趋势, 其中1998—2004年上升显著。

IPCC第五次报告指出, 在过去的50年间全球气温以0.12℃/10a的速度上升^[16]。相关预测随着温度的上升将会导致陆面及近地面空气更加干燥, 进而造成蒸发的上升。然而同时期内多地的蒸发皿蒸散发量观测值呈现出下降趋势, 即所谓的“蒸发悖论”现象^[17]。临泽县2007—2012年以后伴随这气温的上升潜在蒸发量反而下降。左德鹏等^[18]研究表明西北地区蒸发皿蒸发量与PET相关系数高达99%, 因此, 临泽县2007—2012年出现“蒸发悖论”现象。蒸发互补理论假定, 在给定的辐射条件下, 当充分供水时实际蒸散发与潜在蒸散发相等；当下垫面供水量减少时, 实际蒸散发量会减少, 从而释放出更多的能量成为显热, 从而导致潜在蒸散发增加^[19-20]。王艳君等^[21]分析了实际蒸散发和潜在蒸发量的关系, 并利用干燥度指数R(潜在蒸发量与降雨量的比值)来判定研究区域的干湿条件。当R＜0.8时为湿润条件, 实际蒸发量与潜在蒸发量为正相关关系；而R＞1.0时为干旱环境, 实际蒸发量与潜在蒸发量为互补关系。黑河流域属于西北地区典型的干旱地带。参考王艳军等^[21]对区域干湿条件的划分标准, 可知临泽县属于干旱环境(R=5.73)。根据互补理论, 临泽县1980—2012年实际蒸发量持续增加, 显热减少导致潜在蒸散发量减少。

黑河流域是2000年开始实施分水方案, 对流域水资源进行统一管理和调配。各水量要素年际变化的KM分析突变点和转折点都在2000年之后的5年内, 说明受到分水措施的影响很大。李传哲等^[22]研究表明分水后黑河中游地区的耗水量减少；水量的减少导致河道损失量和下渗量减小, 从而导致地下水补给量也相应的减少, 最终导致产水量的下降。

2.6 地表水资源量

地表水资源量为河川径流量总和, 包括临泽县自产的水量和上游河道的来水量。1980—2012年地表水资源量的平均值是(13.17±0.60)亿m³, 最大值出现在1998年, 为21.02亿m³；最小值出现在1980年, 为3.79亿m³。地表水资源量在过去的近33年间表现为微弱的增加趋势, 增长速率为每年0.14亿m³；但没有通过0.05显著性检验。临泽县的地表水资源量主要来自上游的来水, 1980—2012年间的平均来水量为(11.09±0.58)亿m³, 最大值出现在2012年, 为18.95亿m³；最小值出现在1980年, 为2.86亿m³。临泽县自身的产水量多年平均为(2.07±0.04)亿m³。很明显, 上游来水量对地表水资源量的贡献要比临泽县自产量大的多, 其多年平均贡献率的比值约为5:1(图 7), 而且在过去的近33年里波动较小。临泽县属于耗水的中游区域, 水资源主要还是依赖上游, 特别是祁连山区的产流。这说明上游山区对中游区域的生态环境和经济有着重要作用。

3 结论

本文基于黑河中上游土地利用数据、气象数据、水文站实测径流数据、土壤等数据, 利用SWAT模型对临泽县1980—2012年的水文过程进行了模拟, 得到如下结论:

(1) 临泽县1980—2012年均降雨量、实际蒸发量、下渗量、壤中流补给量、地下水补给量分别为117.99、130.53、49.78、68.95、48.00 mm。年均产流量为53.47 mm。产流系数为0.16。产流量主要由壤中流和地下径流控制, 二者的贡献比大约是3:4。

(2) 从空间分布上看, 临泽县的实际蒸发、地下水补给量和产水系数空间差异较大, 三者都呈现“南高北低”的趋势, 即南部灌区的高, 北部非灌区的低。

(3) 临泽县各水量要素年际变化的突变点都在2000年之后的5年内, 说明受到黑河分水措施的影响很大, 特别是下渗量的减少。临泽县2007—2012年出现“蒸发悖论”现象, 根据蒸发互补理论分析, 这种现象是实际蒸发持续增加所导致的结果。

(4) 1980—2012年地表水资源量的平均值是(13.17±0.60)亿m³。上游来水量与自产量贡献率的比值约为5:1。临泽县可用地表水资源更多的依赖于上游的产水。