1.   引言

水文模型是人类认识水循环过程的重要手段，也是水资源管理的重要决策支撑工具。传统的水文模型或重点模拟地表水，而将地下水简化处理，如SWAT（Arnold等，1998）、TOPMODEL（Beven等，1979）、PRMS（Leavesley等，1983）等模型，或重点模拟地下水，而将地表水简单处理，如MODFLOW（Harbaugh，2005）模型。在许多地区，如中国西北干旱、半干旱的内陆河流域，地下水和地表水在生产、生活和生态环境中都发挥重要作用（Schmidt et al., 2013；Simpson et al., 2013；Yin et al., 2013），单一的地表水或地下水模型已不能满足水资源研究和管理的需要。在流域水资源管理过程中，人们逐渐认识到地下水、地表水是一个统一整体，它们会相互影响、相互转化。地表-地下水耦合模型（integrated surface water-groundwater model）可以综合考虑降水、地表水、土壤水、地下水之间的相互作用，全面分析水资源管理问题，包括人类活动和气候变化的影响（Sophocleous，2002）。近年来，越来越多的耦合模型得以开发，代表性的模型包括Hydrogeosphere（Brunner et al., 2012），MIKE-SHE（Graham et al., 2001）、MODHMS（Panday et al., 2004）、PARFLOW（Kollet et al., 2006）、CATHY（Weill et al., 2011）、PAWS+CLM（Shen et al., 2013）、GEOtop（Rigon et al., 2006）、GSFLOW（Markstrom et al., 2008）等。这些耦合模型通常由三维或二维地下水模拟模块和一维或二维的地表水模拟模块组成，能够精细模拟水资源时空分布的演变规律（Wu et al., 2014）。

虽然地表-地下水耦合模型正在迅速发展，但将之真正用于流域尺度水资源管理的案例却不多（Sebben et al., 2013）。原因之一是流域尺度地表水、地下水耦合模型需要大量数据来校准参数和评估模型，这些数据包括地下水水位、河道流量、蒸散发、土壤含水量等。而在大部分流域这些数据十分稀缺。原因之二是耦合模型全面考虑了水循环的主要环节，而模型计算大部分是基于物理过程的，因此计算量十分巨大，计算一次通常需要几个小时到几天时间，故此类模型目前还很难作为日常的管理工具。

近百年来，地球气候系统正经历着显著的变化，气候变化带来的水文效应问题受得各国学者的广泛关注。探讨变化环境下的水循环过程及水资源演变机理，系统评价各种驱动因子的相对作用与贡献对于水资源规划管理具有重要的现实意义。对于西北干旱地区，由于缺水，评估气候变化对水资源的影响尤为重要。

未来中国环境地质研究的重点领域之一就是加强以土壤、包气带、含水层为重点的地球关键带的综合调查监测，充分融合地质科学与地理科学、海洋科学、气象科学、水文科学、生物科学等相关学科，解决更综合、更复杂的生态问题（张永双等，2017）。本为以黑河中游盆地为例，基于地表水-地下水耦合模型，评估区域水资源对气候变化的相应，预测未来气候情景下区域水资源变化趋势，为西北干旱区水资源管理提供参考。

2.   研究区概况

黑河流域是中国第二大内陆河流域，流域总面积约13万km²，跨越青海、甘肃、内蒙古三省。黑河中游盆地是干旱、半干旱流域的典型代表，地表水、地下水之间的互动频繁（Hu et al., 2007），生态环境脆弱。研究黑河中游盆地水循环过程和水资源演变规律对我国面积广阔的干旱、半干旱地区的水资源开发和生态环境保护具有重要意义。

黑河中游盆地位于河西走廊中段（图 1），行政上大部分地区属于甘肃省张掖市，包括山丹、民乐、甘州区、临泽、高台等县市，小部分地区属于酒泉市。盆地范围介于东经98.9°~101.0°、北纬38.4°~ 39.8°，东与马营盆地接壤，西与酒泉西盆地接壤，南起祁连山山麓，北至龙首山、合黎山（统称北山），面积9097 km²。盆地南北宽约60 km，东西长约150 km，是东南-西北走向的狭长盆地。盆地东南海拔2200 m，西北海拔1290 m。主要的土地利用类型是位于祁连山北麓的冲积扇上的荒漠和位于北山南麓黑河两岸的绿洲农田。冲积扇地质结构松散，厚度大（一般厚300~700 m），为戈壁倾斜平原地貌；绿洲地层较薄，主要是细土平原地貌（Wen et al., 2005）。气候属于大陆性干旱气候，具有昼夜温差大、降雨量小、蒸发量大等特点。本地区平均年降雨量约190 mm，主要集中于6—9月，由东向西、由南向北呈递减趋势。平均潜在蒸散发为1325 mm，空间分布与降雨量相反，由东向西、由南向北呈递增趋势（仵彦卿等，2010）。

图  1  黑河中游盆地研究区概况

Figure  1.  General situation of the middle reaches of the Heihe River

下载: 全尺寸图片 幻灯片

盆地的地表水绝大部分来源于祁连山区，本地降水形成的地表产流量很少。汇入张掖盆地的主要河流有黑河干流、梨园河等十余条河流，流出盆地的河流只有黑河干流。梨园河大部分流量被引入灌区，只有在洪水季有流量流入黑河干流；其他河流由于流量小，出山后或是完全渗入山前的冲积扇，或是被引入灌区，没有注入黑河。由2000— 2008年之间的观测可知，全区地表水汇入量大约为每年22.2亿m³，其中17.2亿m³来自黑河干流，2.4亿m³来自梨园河。黑河干流的出山口为莺落峡，梨园河的出山口是梨园堡。黑河干流从正义峡流出盆地，2000—2008年的平均年径流量为9.7亿m³。

盆地地下水流向为东南至西北，南部山前地下水埋深可达200 m以上，绿洲的埋深逐渐变为几十米至几米。盆地地表水和地下水转换频繁，黑河出山后流经山前冲积扇大量入渗，补给地下水，至细土平原地区后转为接受地下水补给（Chen et al., 2006）。除了河道接受地下水补给外，绿洲区还有大量泉水溢出，被用于灌溉或汇入河道。前人的研究表明，地下水补给地表水的量大约为每年9.9亿m³，地表水补给地下水的量大约为每年12.3亿m³（Hu et al., 2007）。盆地的地下水补给除了地表水外还有祁连山和北山的地下水侧向补给，其补给强度由于缺少观测数据尚无定论。流域出口正义峡附近为基岩，因此盆地内无显著地下水流出。

3.   模型概况

本研究选择GSFLOW（Coupled Ground-Water and Surface-Water Flow Model）作为地表-地下水耦合模型的构建平台。GSFLOW耦合了PRMS（Precipitation- Runoff Modeling System）和MODFLOW两个模型。PRMS是经典的分布式地表产汇流模型，MODFLOW是经典的地下水模型。GSFLOW模拟了地表、地下水系统中的主要水文过程，如坡面漫流、地表径流、蒸散发、入渗、壤中流、非饱和带水分运移、饱和带地下水流动、地下水排泄、河道和地下水的相互补给等等。

GSFLOW中地表部分计算单元是水文响应单元HRU（Hydrologic Response Units），地下部分计算单元是有限差分网格。在地表土壤带和地下饱和带之间是包气带，包气带由UZF模块（Unsaturated Zone Flow package）来计算。土壤水向下渗漏到重力水水库中，当地下水水位低于阈值时，重力水水库向包气带释放水量，这些水会向下运动到饱和带；当地下水水位高于阈值时，饱和带的水会补充到重力水水库，最终成为土壤水或壤中流的一部分。UZF模块使用垂直方向的一维理查德方程求解，主要参数为包气带的水力传导系数和饱和含水率。

数字高程（DEM）数据用来进行地表流域划分和定义地下网格顶部和河床底部的高程。土地覆盖和土壤类型数据用来进行HRU的划分和HRU参数的确定。水文地质信息，如剖面图，水文地质图等，用于进行地下网格的分层和水文地质参数的分区。气象信息包括降雨、日最高气温、日最低气温。使用6个气象站点（图 1）的气象数据进行插值，把气象数据分配到每个HRU上。引水信息和抽水信息按灌区收集。共有4个水文站，莺落峡和梨园堡两个站位于模型的边界上，其流量数据作为模型的边界条件。高崖和正义峡两个站的流量数据用于模型校正。地下水水位数据共35组，观测井位置如图 1所示。地下水水位数据也被用于模型校正。

建模范围参考了Hu et al.（2007）的工作。如图 1所示，模型主要包括张掖盆地和酒泉东盆地两部分，以张掖盆地为主。模型东部与马营盆地接壤，两盆地之间的地形有隆起，是天然的地下水隔水边界，也是地表水分水岭。马营盆地与张掖盆地在南部通过浅层地下水相连，有微弱的地下水侧向流动；在北部通过山丹河相连，但山丹河已经断流，所以没有水力联系。因此模型东部的边界条件是无流边界和定流量边界。模型西部与酒泉西盆地接壤。在两盆地接壤地界是一条南北方向的地下断裂，也可以看作是地下水隔水边界。因此模型西部的边界条件是无流边界。模型南北两侧是山区，接受地下水和地表水的侧向补给。地表水的边界条件根据河流规模选择了两条主要河流——黑河干流和梨园河，在模型校正的过程中不对其流量进行调整。其他河流流量整合到地下水浅层边界入流部分。地下水边界条件根据边界附近的水力梯度、水力传导系数和地层厚度进行估算（仵彦卿等，2010），作为模型的初始值，且在模型校正时对其进行调整。模型出口是正义峡，从正义峡到模型西南角的地下部分都是基岩，因此没有地下水出流，地下水在正义峡附近河道全部转化为地表水，然后由正义峡流向下游。

使用DEM数据把地表部分共划分了104个子流域，588个HRU。地下部分的分层参考了前人的研究（张荷生等，2003；Wen et al., 2007；仵彦卿等，2010）。地下部分分为5层，其中第1、3、5层为含水层，第2、4层为隔水层。每层有150行，172列，9106个活动网格。网格大小为1 km×1 km。最顶部的第1层根据水文地质条件划分为21个参数分区，第2、4层分为12个参数分区，第3、5层分为14个参数分区。每个参数分区对应一套水力传导系数、储水系数或给水度。子流域划分情况见图 2。

图  2  子流域划分、数字河道及HRU情况

Figure  2.  Sub-basin division, digital rivers and HRU

下载: 全尺寸图片 幻灯片

模型校正的过程使用手动调参的方式进行。经过敏感性分析，筛选出主要敏感参数进行调参，包括最大蒸发深度、河床水力传导系数、河床糙率、地下水平水力传导系数、垂向水力传导系数、给水度、地下水侧向补给强度等。校正量为河道出流量、地下水水位。调参结果如图 3、图 4所示。

图  3  盆地出流量（正义峡）校正结果（亿m³）

Figure  3.  The calibration result of basin outflow (Zhengyixia station) (10⁸ m³)

下载: 全尺寸图片 幻灯片

图  4  地下水水位校正结果

a—埋深；b—等水位线；c, d—井A与井B（位置见图 1）两口典型监测井水位动态

Figure  4.  The calibration result of groundwater level

a-Average depth of groundwater table; b- Contour of groundwater levels; c, d-Seasonal variation of the groundwater level at well A and well B

下载: 全尺寸图片 幻灯片

图 5是区域水平衡中几个重要通量的模拟值。可以看出，盆地的水资源主要来源于上游来水，其次是降雨。主要耗水项为蒸散发，其次是流域出流量。盆地地表、地下水交换频繁，地表水向地下水的补给（包括河道入渗、降雨入渗及灌溉入渗）可达每年9.7亿m³，地下水向地表水排泄（包括河道出漏、泉水出漏）可达每年10.6亿m³。盆地目前地下水储量总体呈现减少的趋势，速度为每年1.7亿m³。

图  5  区域水循环通量

Figure  5.  Water budget of the study area

下载: 全尺寸图片 幻灯片

对模型的主要输入参量和输出量进行不确定性分析，分析输出量变化范围及输入量的贡献度。分析方法为Sobol方法。该方法最早由Sobol于1993年提出, 其核心是将模型输出的总方差可分解为每个参数的方差和参数相互作用的方差之和, 然后通过参数对输出方差的贡献比例进行敏感性分级。

通过不确定性分析可知，对于不同的输出变量，其不确定性的主要来源不同，但起主导作用的是少数（本案例中不超过3个）。蒸散发模拟值的主要不确定性来源为降水数据和灌溉数据；土壤含水量的主要不确定来源为最大蒸发深度；出口河流流量的主要不确定性来源为引水及河道入流量；地表、地下水交换量的主要不确定性因素为边界流量、河道入流量和河道渗透系数。最后，对于±20%的输入变量不确定性，不同输出量的变异系数不同。蒸散发的变异系数为5.6%，土壤含水率为8.9%，出口河流流量为14.5%，地表、地下水交换量为23.9%。这表示模型不同的输出量，其不确定性存在显著差异。

4.   未来气候分析

全球气候模式和区域气候模式是进行过去和未来气候变化研究的重要工具。近几年来随着计算机条件的改进，全球气候模式的分辨率得到了进一步的提高，为未来气候预测与水资源模拟的研究提供了方便。

耦合模型比对项目（CMIP5，https://esgf-node.llnl.gov/projects/cmip5/）是在世界气候研究计划（WCRP）框架下，综合对比全球各个气候实验室的大气-海洋环流模型输出结果的项目。CMIP5提供各个气候模式的诊断、验证、比对、文档和数据下载服务。在IPCC（联合国政府间气候变化专门委员会，Intergovernmental Panel on Climate Change）第五次报告中，对于碳排放预算是提出了一种新的场景假设，即代表路径浓度（RCP）。排列分别为RCP2.6、RCP4.5、RCP6.0和RCP8.5，RCP后面的数字表示到2100年辐射强迫水平, 单位为Wm^-2。笔者选取了9个CMIP5的气候模式（表 1），分析两种人类社会发展情景——中等温室气体排放情景（RCP4.5）和高排放情景（RCP8.5）下，黑河中游盆地未来百年间（至2100年）温度、降水的变化情况。

表  1  9种全球气候模式及机构

Table  1.  Descriptions of the 9 models from CMIP5 applied in this study

下载: 导出CSV 

| 显示表格

表 2为9个气候模式的历年预测降水量平均数据与实际观测数据（酒泉、张掖站平均值）的对比情况。预测降水量基本能够反应实际降水量的年际变化，但绝对值的大小有所差距。预测降水量的平均值为128 mm，观测值的平均值99 mm，误差为29%。

表  2  历年降水量与历年预测降水量

Table  2.  Precipitation and predicted precipitation over years

下载: 导出CSV 

| 显示表格

图 6和图 7显示未来百年的两种温室气体排放模式下，降水、温度的变化情况。灰色实线代表 9个独立气候模式的预测，黑色实线表示9种模型的平均值，黑色虚线表示9种模型平均值的线形回归值。总体来看，未来降水和气温都程上升趋势。在中等排放情景下，降水平均每年上升0.6 mm，温度平均每年上升0.03℃；在高排放情景下，降水平均每年上升0.8 mm，温度平均每年上升0.06℃。对比IPCC第五次评估报告中的数据，全球平均气温在中等排放情景下每年上升幅度约为0.02℃，在高排放情景下每年上升幅度约为0.04℃，降雨则根据地区的不同有所差别，潮湿与干旱地区之间的对比会增加。由此可见，本地区的气候变化烈度是超过全球平均水平的。全球气候模式模拟的日数据或月数据本身是没有意义的，对未来气候预测有意义的是模拟值的统计量。

图  6  未来百年降水变化情况

a—中等排放强度（RCP4.5）；b—高排放强度（RCP8.5）

Figure  6.  Precipitation changes in the coming one hundred years (a-RCP4.5; b-RCP8.5)

下载: 全尺寸图片 幻灯片

图  7  未来百年温度变化情况

a—中等排放强度（RCP4.5）；b—高排放强度（RCP8.5）

Figure  7.  Temperature changes in the coming one hundred years (a-RCP4.5; b-RCP8.5)

下载: 全尺寸图片 幻灯片

5.   未来水资源变化趋势

根据9个气候模式下未来气候变化情况，取平均变化量，利用耦合模型计算2019—2100年黑河中游盆地地下水、地表水资源的变化情况。地表水引水量、地下水抽水量、农田灌溉量维持2012年的统计数据。

图 8a表示两种气候情景下，地下水储量的变化情况。总体来看，两种气候情景下地下水都程下降趋势，这表明降水的增加不足以补给地下水的减少量。在中等排放强度下，地下水储量平均每年减少0.38亿m³。在高排放强度下，地下水储量平均每年减少0.34亿m³。在2091年之前，中等排放情况下，地下水减少趋势比高排放情况更快，但这种情况在2091年左右颠倒。原因在于降水的增加可能使地下水存量增加，但气温的增加有可能增加地下水的蒸发，从而减少地下水存量。在未来气候情景下，这两种情况同时发生，要看降水增加和气温增加哪个因素占据主导地位。此时需要使用耦合模型来计算地下水综合的变化情况。根据模型，两种情景下都是气温增加占主导地位，所以地下水会减少。但随着地下水蒸发越来约多，可供蒸发的地下水会减少，此时降水的影响会增加，减缓了地下水减少的速度。

图  8  两种气候情景下水资源变化情况（a—地下水储量；b—研究区出流量）

Figure  8.  Changes in water resources under two climate scenarios (a-Groundwater storage; b-Basin outflow)

下载: 全尺寸图片 幻灯片

图 8b表示两种气候情景下，研究区出流量的变化情况。总体来看，两种情况出流量随时间变化都不明显，中等排放情景下出流量呈现微弱减少趋势，平均值在15亿m³/a左右，高排放情景下出流量呈现微弱增加趋势，平均值在13.7亿m³/a左右。这种差别的原因是来自上游来水的变化。在中等排放情景下，上游山区产流量大于高排放情景。这表明对与干旱区山区产流，增温带来的蒸发的影响大于降雨的影响。

综合分析，在未来百年的时间尺度下，高排放的情景对地下水储量更有利，能减缓地下水减少趋势，但流域出流量更少；从更长远的时间来看，结果会相反。

图 9显示的是两种气候情景下，两个典型井的水位年际变化情况。井位A在上游地下水补给区，在两种排放情景下水位均呈现上升趋势，且高排放情景的水位上升趋势小于中排放情景，原因在于地下水补给区埋深较深，受到蒸发影响小，受到降水影响大，造成水位上升；同时受到河道流量的影响，高排放情景下上升趋势较小。井位B在中游地下水排泄区，在两种排放情景下水位均呈现下降趋势，且高排放情景的水位下降趋势大于中排放情景，原因在于地下水排泄区埋深浅，受到蒸发影响大，造成水位下降；同时受到河道流量的影响，高排放情景下下降趋势较大。

图  9  两种气候情景下典型井水位变化情况

a—井位A；b—井位B

Figure  9.  Annual variation of the groundwater level under two climate scenarios

a-Well A; b-Well B

下载: 全尺寸图片 幻灯片

6.   结论

本文建立了黑河中游盆地地表-地下水耦合模型，分析了未来（2019—2100年）气候情景下盆地区降水、气温的变化情况，并根据模型分析了未来气候情景下地下水、地表水资源的变化情况。

（1）黑河中游盆地地表-地下水模型的建立和校正结果表明，GSFLOW模型能很好的模拟黑河中游盆地复杂的水循环过程。盆地地表、地下水交换频繁，地表水向地下水的补给可达每年9.7亿m³，地下水向地表水排泄可达每年10.6亿m³。盆地目前地下水储量总体呈现减少的趋势。莺落峡出山口一带的地下水水位下降明显，最大下降速度可达2.57 m/a。

（2）全球气候模式结果分析表面，盆地区降水、气温在未来变化显著，且变化幅度超过全球平均水平。在中等排放情景下，降水平均每年上升0.6 mm，温度平均每年上升0.03 ℃；在高排放情景下，降水平均每年上升0.8 mm，温度平均每年上升0.06 ℃。

（3）结合黑河中游盆地地表-地下水模型和全球气候模式结果，在目前的农业活动强度下，盆地区未来水资源变化不容乐观。在中等排放强度下，地表水流量平均在15亿m³/a左右，地下水储量平均每年减少0.38亿m³。在高排放强度下，地表水流量13.7亿m³/a左右，地下水储量平均每年减少0.34亿m³。

本文探讨了结合地表-地下水耦合模型和全球气候模式的可能性。二者结合可以模拟未来气候变化下水资源的变化情况。本文的灌溉、引水、抽水强度设定为2012年的强度，以后的研究可以分析在不同的人类活动条件下未来水资源的变化情况，为水资源管理提供参考。