基于SWAT-MODFLOW的黄河中游区径流过程模拟及对黄土高原变绿的响应

张洪波, 支童, 卫星辰, 党池恒, 夏岩, 高文冰

(1.长安大学 水利与环境学院,陕西 西安 710054; 2.长安大学 旱区地下水文与生态效应教育部重点实验室,陕西 西安 710054;3.中国石油长庆油田公司勘探开发研究院,陕西 西安 710021; 4.低渗透油气田勘探开发国家工程实验室,陕西 西安 710021)

近年来,降水、气温等气候要素的趋势性变化以及极端事件(干旱、暴雨)发生频率与强度的变化,对流域水文循环过程造成了一定的影响,同时大规模水利工程建设、水土保持和植树造林工程、直接取用水等人类活动也在加剧流域下垫面的改变,影响流域截留、入渗、蒸散发、径流等产汇流过程[1],从而显著减少了流域地表径流量及其动能和势能[2]。如在我国黄河中游的黄土高原区,退耕还林还草政策的实施,使得1999—2010年间黄土高原累计造林面积达到1 890.6×104hm2 [3]。有研究表明,2000—2014年黄土高原植被覆盖指数持续显著增大,增速达6.93%/(10 a)(P<0.01)[4]。可见,“黄土高原变绿”已成为不争的事实。然而数据显示,在近50年来，黄土高原区生态效应显著提高的同时,黄河流域的径流量却呈现锐减趋势[5]。对比1956—2000年与2001—2017年两个不同时期,可发现黄河中游后一时期的平均降水量比前一时期增加了仅3.1%,但地表水资源量却减少了21%[6]。由此可见,降水因素并不是造成黄河水资源大幅减少的直接原因,黄土高原林草植被改善和大规模梯田淤地坝的运用应该是造成径流量显著减少且不容忽视的影响因素。

有关土地利用/覆被变化水文效应的研究自20世纪90年代兴起,重点对焦变化环境下的水循环问题。由于该问题的复杂性与交叉性,目前乃至未来几十年仍将是全球研究的热点和前沿问题[7]。对于黄河径流演变及其与变化环境的关系,国内外学者已开展了大量的研究,其方法可主要归纳为3类,即趋势分析法、水文模型法和概念性方法[8]。趋势分析法常基于径流和影响因素的回归关系,例如累积量斜率变化率比较法和双累积曲线法[9-10];水文模型法主要依托集总式水文模型和分布式水文模型[11-12];而概念性方法多基于Budyko假设和Tomer-Schilling框架展开,如弹性系数法和水文敏感性分析法等[13-14]。其中,分布式水文模型凭借其较强的物理机制,已广泛应用于流域径流过程模拟及对环境要素的响应分析中。SWAT模型是典型的分布式水文模型之一[15],在我国多个流域都取得了较好的应用效果。但由于黄土高原地区的地下水补给较为强烈,基流量占到了径流总量的50%～80%,而传统的SWAT模型又在地下水模拟方面能力略显不足,故在该地区的应用相对有限。为解决这一问题,本文拟引入在地下水模拟方面表现较为出色的MODFLOW模型,通过构建SWAT-MODFLOW地表水-地下水耦合模型,弥补地表水模型主要对地表及非饱和土层中水流运动过程进行模拟,而对含水层中的水流运动过程简单处理(一般基于水量平衡和简单调蓄作用来计算)的缺点,同时规避地下水模型在模拟土壤水下渗过程(包括水流在包气带、潜水层及承压水层中的运动)中涉及与地表水的交换过程时,对交换水量的计算往往采用经验参数进行简单计算的不足。通过对黄河中游区径流过程的模拟,重点探索“黄土高原变绿”情况下气候和下垫面对流域水文系统的关键性驱动过程,为变化环境下的流域水资源管理和水安全保证提供技术支持与数据基础。

延河是黄河中游的一级支流,也是陕北黄土丘陵沟壑区的重要水源,其重要性不言而喻。受气候变化、人类取用水活动以及“黄土高原变绿”的影响,延河的径流量呈现逐年减少趋势[16]。鉴于此,本文以黄河中游区延河流域延安水文站以上控制流域以及西川河子流域为研究区,通过模拟1964—2019年杏河站、安塞站、枣园站和延安站4个水文站的近天然径流(即径流不受人类活动直接或间接的影响)和高原变绿条件下的自然径流(即径流仅受间接人类活动引起的下垫面变化的影响,而不受人类直接取用水活动的影响),研究“黄土高原变绿”对流域径流过程的影响,旨在为变化环境下的流域水资源评价与科学管理提供参考。

1 研究区概况与数据来源

延河发源于靖边县天赐湾乡周山,由西北向东南流经志丹、安塞、延安区域,最终于延长县南河沟凉水岸附近注入黄河。流域地势西北高、东南低,水系结构呈树枝状。径流年内分配不均匀,主要集中于夏季(6—9月),冬季径流最少。延河流域属中温带大陆性半干旱季风气候,年平均气温8.8 ℃,年平均降水量505 mm[17]。

本文研究区为延河流域延安水文站以上流域及其一级支流西川河流域(以下简称延河流域)。延安水文站的控制流域面积为3 188 km2,多年平均径流量为1.21×108m3。枣园水文站位于西川河下游,系流域控制站,控制流域面积为713.2 km2,多年平均径流量为0.22×108m3。研究区水系分布及各水文站点位置如图1所示。

采用的数据资料主要包括:30 m分辨率的数字高程模型(DEM,来源于地理空间数据云);1 km分辨率的全国4期土地利用类型数据(19世纪70年代末、1995年、2005年、2015年,来源于中国科学院资源环境科学数据中心);土壤类型空间分布数据和属性数据(来源于世界土壤数据库,HWSD);延安站、靖边站、横山站、绥德站、延长站1962—2019年的气象数据(日降水量、日最高/最低气温、日平均风速、日相对湿度,来源于中国气象数据网中国地面气候资料日值数据集);杏河站(1959—2016年)、安塞站(1981—2016年)、枣园站(1971—2014年)、延安站(1959—2014年)的实测日径流数据(来源于黄河水利委员会水文数据年鉴)。

图1 研究区水系分布及站点位置

2 研究方法

SWAT-MODFLOW耦合模型是目前较为流行的地表水-地下水耦合模拟模型,具有输入数据量小、易获取、耦合操作简单等优点,已受到国内外水资源领域学者的广泛关注[18-20]。该模型是由美国农业部开发出的SWAT模型和美国地质调查局开发的MODFLOW模型耦合而成。其特点在于充分利用了MODFLOW模型在地下水数值模拟与SWAT在地表水数值模拟方面的优势,通过GIS以数据共享的方式将其耦合,并实现数据交互,不仅大大提高了水资源模型的模拟精度,也大大降低了耦合操作的难度。

2.1 SWAT模型

SWAT模型属于典型的分布式水文模型[21],主要应用数值分析来建立相邻水文响应单元(HRUs)之间的时空关系,即在每个HRU上应用传统的概念性模型来推求净雨,再进行汇流演算,最后模拟出口断面的流量过程[22]。

基于SWAT的流域水文过程模拟主要涉及水循环的陆面部分和水面部分。陆面部分(即产流和坡面汇流)主要控制各子流域内主要河道水的输入量。水面部分(即河道汇流)决定水由河网向流域出口的输移过程。其模型水文循环模拟简易单元结构如图2所示。

图2 SWAT模型水文循环简易单元结构

SWAT模型的水文模拟模块主要基于水量平衡方程:

(1)

式中:SWt为土壤最终含水量,mm;SW0为第i天的土壤初始含水量,mm;Rday为第i天的降水量,mm;Qsurf为第i天的地表径流量,mm;Ea为第i天的蒸散发量,mm;Wseep为第i天离开土壤剖面底部的渗透水流和旁通水流的水量,mm;Qgw为第i天的回归流水量,mm。

其中,水循环的陆面部分主要包括地表径流计算、壤中流计算、地下水计算和蒸发计算,其中地表径流计算一般采用SCS径流曲线法[23],入渗采用存储演算方法计算,地下水计算应用潜水和承压水的水量平衡方程,潜在蒸发能力一般采用Penman-Monteith法[24]。水循环的水面部分多采用曼宁公式和马斯京根法进行汇流计算。

2.2 SWAT-MODFLOW模型

SWAT-MODFLOW耦合模型的基本过程是通过水文响应单元(HRUs)和有限差分网格(Cells)之间的映射关系,将SWAT模拟的深层入渗量作为面状补给分配到MODFLOW模型的相应网格单元上,以作为地下水模拟的边界条件。同时,将MODFLOW模型所计算出的地下水排泄量,通过映射关系添加到SWAT模型的相应子流域流量中,从而实现SWAT模型与MODFLOW模型之间的数据双向传递和耦合。

SWAT-MODFLOW耦合模型主要通过QGIS软件的QSWATMOD插件来实现。2017年,Seonggyu Park和Ryan T Bailey采用FORTRAN语言编写并开发了SWAT模型与MODFLOW模型的耦合程序,即“SWAT-MODFLOW.exe”,QSWATMOD插件则正是基于这个耦合程序而创建的。

2.3 模型校准

SWAT模型通常采用SWAT-CUP实施模型校准。SWAT-CUP是由瑞士联邦水质科学技术研究所专门针对SWAT模型开发的,实现了自动参数率定,可给出最优参数值和最优参数的取值范围。其中,SUFI-2方法由于考虑了模型结构、监测数据、模型参数等不确定来源,因此本文拟选择SUFI-2方法进行流域参数率定。

MODFLOW模型校准多以基础水文地质参数作为初始值,在满足水文地质条件定性认识的基础上,不断调整水文地质参数值,直至得到既能符合水文地质条件的定性认识,又能最佳拟合径流观测值资料的水文地质参数。

本文选取决定系数R2和纳什效率系数NS作为模拟流量和实测流量的拟合程度评价指标,其计算公式为:

(2)

(3)

式中:Qm,i为实测流量,m3/s;Qp,i为模拟流量,m3/s;Qm,avg为多年实测平均流量,m3/s;Qp,avg为多年模拟平均流量,m3/s;n为实测时间序列长度。一般认为,当R2>0.60且NS>0.50时,模型拟合程度满意。

3 模型构建

3.1 模型数据库

SWAT模型运行所需的数据主要包括流域空间数据、属性数据、气象数据和水文数据。其中,空间数据包括流域DEM、流域土地利用类型数据和流域土壤类型数据。DEM主要用于生成水系、出水口点和子流域划分等;土地利用数据和土壤数据主要为生成水文响应单元(HRUs)和模型产流汇流计算服务;土壤属性数据和气象数据需要按照模型输入手册的格式进行数据库构建。

3.1.1 土地利用数据

本文选择4期土地利用类型数据(20世纪70年代末、1995年、2005年、2015年)作为模型现状下垫面条件下径流模拟的输入,利用ArcGIS裁剪出研究区土地利用类型数据,并进行土地利用类型重分类。

各类型土地利用面积统计见表1。由统计数据可知,在时域上即从20世纪70年代末期到2015年,延河流域的耕地面积呈现较为明显的减少趋势。其中,20世纪70年代末期到2000年延河流域耕地面积的变化相对较为微弱,2000年以后减少趋势逐渐凸显;林地面积从20世纪70年代末到2015年呈现连续增长趋势,同样在2000年以后林地面积变化显著;同期草地面积呈现类似的减少趋势。

表1 延河流域土地利用类型面积统计 km2

延河流域各期土地利用类型分布如图3所示。由图3可以看出,在空间上,耕地、林地和草地的面积变化较为均匀,未出现聚集性特征。时间上,延河流域在2000年后退耕还林还草效果明显,流域下垫面变化显著。

图3 延河流域各时期土地利用类型空间分布

3.1.2 土壤数据

模型运行所需的土壤数据包括土壤类型的空间分布数据和土壤的属性数据。土壤属性表主要字段包括土壤名、顶层土壤质地、土壤参考深度、土壤有效持水量、土壤相位、到土壤底部存在障碍的深度分类、土壤含水量特征、碎石体积百分比、砂石含量、淤泥含量、黏土含量、土壤容重、有机碳含量、酸碱度、土壤的阳离子交换能力、碳酸盐或石灰含量、硫酸盐含量、电导率等参数,其中以T开头的属性字段表示上层土壤属性(0～30 cm),以S开头的属性字段表示下层土壤属性(30～100 cm)。由于FAO-90土壤分类系统的土壤分级采用的是美国农业部USDA分级制,不必要再对土壤粒径进行转换,可利用GIS掩膜工具裁剪出研究区土壤类型图,实施土壤类型重分类,最后依据SWAT土壤属性数据库格式要求建库。其中,土壤有效持水量、土壤湿密度、土壤饱和导水率使用SPAW软件计算;土壤侵蚀力因子和最小下渗率采用经验公式计算。

延河流域各类型土壤面积统计见表2,土壤类型空间分布如图4所示。

表2 延河流域土壤类型面积统计表

图4 延河流域土壤类型空间分布

由表2可知,整个流域内黄绵土所占面积比重最大,其次是珊瑚砂土和新积土,所占面积比重最小的是积钙红黏土和钙质粗骨土。

由图4可知,在空间分布上,河道之外以黄绵土为主,而在河道内以珊瑚砂土和新积土为主,且新积土主要分布在杏子河和西川河支流上。积钙红黏土和钙质粗骨土占据面积极为有限。

3.2 子流域划分与模型耦合

经过反复调试,设定延河流域集水面积阈值为6 500 hm2,添加具有监测数据的杏河、安塞站为控制点,延安和枣园水文站(西川河)为流域出口点,最终将延河流域划分为33个子流域,如图5所示。设定延河流域土地利用面积阈值为5%,土壤类型面积阈值为5%,坡度阈值为5%,最终将延河流域划分为615个HRU,如图6所示。

图5 延河流域子流域划分

图6 延河流域HRUs划分

目的含水层设定为埋藏在80 m以下、100 m以上的地下水含水层。MODFLOW模型则根据区内地质、水文地质特征,对研究区进行均匀网格剖分。延河流域设置网格大小为300 m×300 m,活动单元格数为44 905个,如图7所示。

图7 延河流域MODFLOW网格剖分

SWAT-MODFLOW模型通过构建SWAT模型的HRUs和MODFLOW模型的cells的映射关系来达到两个模型的耦合。在创建完4个链接文件的基础上,自动生成swatmf_link.txt文件,即可运行耦合后的SWAT-MODFLOW模型。

3.3 模型率定和验证

输入气象数据后,设置率定期为1964—1968年,验证期为1969—1970年,模拟时间尺度为月尺度,预热期长度为2 a,进行SWAT模型的率定和验证。杏河站、安塞站、枣园站、延安站率定的最佳参数见表3。

表3 延河流域参数率定

MODFLOW模型的调参主要是通过参数反演和人工调参相结合,对各水文地质分区的渗透系数(K)及给水度(μ)进行调整,最终结果见表4。

表4 MODFLOW模型水文地质参数率定

对于SWAT-MODFLOW耦合模型,流域月均流量模拟值为SWAT模型地表模拟结果和MODFLOW地下水补给河流模拟结果之和。由于安塞站实测数据从1981年开始,枣园站实测数据从1971年开始,因此在延河流域率定期和验证期(1964—1970年)仅采用杏河站和延安站模拟月均流量与实测月均流量进行模型评价,所得站点模拟结果的评价指标R2和NS见表5。由表5可知,各站点模拟评价指标均符合拟合要求,证明所构建的SWAT-MODFLOW耦合模型可用。其中,位于支流上的杏河站模拟结果较好,这与支流上人类活动对径流影响较小相符合。另外,模型在流量峰值模拟上表现欠佳,这可能是由黄土高原地区复杂的产汇流机制所导致的。两个站点率定期和验证期的实测值和模拟值拟合情况分别如图8和图9所示。

表5 延河流域模型模拟结果

图8 杏河站率定期和验证期的模拟值和实测值对比

图9 延安站率定期和验证期的模拟值和实测值对比

4 径流模拟

4.1 近天然径流模拟

采用1964—2019年长序列气象数据和20世纪70年代末期的土地利用数据,模拟流域近天然条件下的月流量过程(即径流不受人类活动直接或间接的影响)。杏河站、安塞站、枣园站、延安站天然条件下年径流量模拟值如图10所示。由图10可见,近天然条件下的年径流年际变化基本呈现丰枯交替的平稳趋势,但受降水量减少的影响,近天然径流有微弱的下降趋势。

图10 天然条件下各站点年径流模拟值

4.2 高原变绿条件下的径流模拟

采用1964—2019年长序列气象数据和1970s及1995、2005、2015年的土地利用数据,对流域高原变绿条件下(即现状下垫面条件下)的月流量过程进行模拟(即径流仅受间接人类活动引起的下垫面变化的影响而不受直接人类取用水活动的影响)。其中,1991年之前的模拟,采用1970s的土地利用数据;1991—2000年的,采用1995年土地利用数据;2001—2010年的,采用2005年土地利用数据;2011—2019年的,采用2015年土地利用数据。高原变绿条件下,杏河站、安塞站、枣园站和延安站的年径流量模拟值如图11所示。由图11可见,高原变绿影响下,延河流域各站的年径流量均呈现出明显的下降趋势,说明间接人类活动导致的流域下垫面变化确系流域年径流量减少的主因之一。

图11 高原变绿条件下各站点年径流模拟值

4.3 高原变绿影响的定量分析

绘制延河流域杏河站、安塞站、枣园站、延安站近天然条件、高原变绿条件下的模拟年径流量序列以及实测年径流量序列的对比箱图，如图12所示。由图12可发现,延安站径流量的年际变化较大,且延河流域在近天然条件、高原变绿条件下的模拟年径流量和实测年径流量呈现逐级消减态势。这说明间接人类活动导致的流域下垫面变化和直接人类取用水活动均对流域径流量的缩减有所贡献。

基于上述不同站点径流量变化间的差值,可统计分析不同人类活动影响因素的贡献量。这里需要说明的是,由于模型模拟结果中考虑了气候因素的变化,因此所得的差值仅反映直接和间接人类活动的影响量。其中,间接人类活动对年径流量的影响可设定为近天然条件下的年径流量与高原变绿条件下的年径流量的差,直接人类活动对年径流量的影响可设定为高原变绿条件下的年径流量与实测年径流量的差。遂可得间接人类活动和直接人类活动对流域人类活动所驱动的年径流量减少的贡献量,详见表6。

图12 各站点模拟年径流量及实测年径流量分布箱图

表6 延河流域间接和直接人类活动对年径流量减少的贡献

由表6可知,不同站点控制流域范围内的间接和直接人类活动的影响不同,其对径流缩减的贡献也不同。杏河站、安塞站和延安站的间接人类活动的影响略大于直接人类活动的影响,枣园站所表征的西川河流域直接人类活动影响强烈,贡献了人类活动所驱动的年径流量缩减的70.97%。以上结果表明,高原变绿对流域径流量的影响不容小觑,已成为驱动黄河中游区径流衰减的重要因素,应给予更多的关注,并开展有效的支撑性研究。

5 结语

本文以黄河中游区延河流域延安水文站以上流域及西川河子流域为研究区,构建了SWAT-MODFLOW耦合模型,模拟了1964—2019年杏河站、安塞站、枣园站和延安站近天然条件下与高原变绿条件下的径流变化,分析了“黄土高原变绿”对流域径流过程的定量影响,得到的结论如下:

1)空间数据分析表明,研究区土地利用类型以耕地、林地和草地为主。从20世纪70年代末期到2015年,研究区耕地和草地的面积呈现减少趋势,且2000年后减少趋势显著;林地面积呈现连续增长趋势,同样2000年后增加显著。这说明,研究区在退耕还林还草工程的影响下,流域下垫面变化显著。

2)模型校验结果表明,SWAT-MODFLOW耦合模型在研究区具有适用性,纳什系数均在0.71以上,可有效开展黄河中游区特别是黄土高原区的水文模拟研究。

3)模型模拟结果表明,近天然条件(不考虑人类取用水与下垫面变化的影响)下,研究区年径流基本呈现丰枯交替的平稳趋势;而在高原变绿条件(考虑下垫面变化的影响,不考虑人类取用水的影响)下,年径流则呈现出较为明显的缩减变化,从而验证了研究区间接人类活动所导致的下垫面变化系研究区年径流减少的重要主因。

4)定量分析表明,杏河站、安塞站和延安站的径流受间接人类活动的影响略大于受直接人类活动的影响,而西川河流域中直接人类活动的影响相对强烈,贡献了人类活动所驱动的年径流量缩减的70.97%。

通过以上结论可知,变化环境下的流域水循环是一个极其复杂的过程,特别是在黄河中游的黄土高原区,气候变化、下垫面改变和人类直接活动交织影响着区域水文循环过程,驱动着地表水资源量的变化。且由于影响因素的非平稳性变化,导致区域径流的不确定性日益增大。如何认识、掌握黄河中游黄土高原区气候-下垫面-人类活动-水文系统间的驱动过程,已成为保障区域水安全和促进黄河流域高质量所迫切需要解决的问题之一,亟须开展科学深入的基础性研究,以期为变化环境下的区域水资源管理和社会经济高质量提供参考依据。