郑 丽,金 鑫,3,*,金彦香,3,傅 笛,翟婧雅

1 青海师范大学地理科学学院,西宁 810016 2 青海省自然地理与环境过程重点实验室,西宁 810016 3 高原科学与可持续发展研究院,西宁 810016

地下水在干旱区内陆河流域水循环中占有主导地位,是自然系统中最重要的基础性资源之一[1—3]。植被是土壤和大气间水量交换的关键驱动因素,其通过降雨截留和消耗土壤水分等方式间接影响地下水的补给,在地下水循环中具有重要作用[4—5]。相关研究表明,在冰岛,当草地及荒地被人工林代替后,植被覆盖增加,区域地下水补给量降低了10%[6]。在荷兰,当农田被林地代替后区域地下水补给量减少约485 mm/a[7]。澳大利亚昆士兰州东南部部分流域的稀疏草地被稀疏林地替代后,地下水补给量将减少25%；当稀疏草地被茂密林地替代后,地下水补给量将减少 48%[8]。我国黄土高原泾河流域的研究表明,平均森林覆盖率每增加10%,将导致流域年地下水补给量低11.1 mm[9]。这些区域尺度的研究已经证实,植被覆盖增加对地下水补给的影响非常显著。

一般来说,干旱/半干旱区的地下水补给比湿润地区更容易受到地表覆盖条件的影响[10—11]。由于降水少且季节性强、补给量通常很小等原因,干旱/半干旱区地下水补给量的估算具有挑战性,一般的土壤水资源平衡法[3]、水文资料分析法[8]等难以适用。目前,水文模型因其可对自然界中复杂水循环过程近似描述和再现,可作为区域植被变化之水文效应研究的有效手段[12—14]。如SWAT(Soil and Water Assessment Tools)模型,在我国干旱/半干旱区内陆河流域的水循环及其影响因素研究中具有广泛应用[15—17]。区域植被覆盖的增加或减少可能由人类活动导致的土地利用类型变化引起,也可能因气候变化导致的植被生长状态改变等因素引起。SWAT模型中的土地利用更新模块可在一定程度上表达土地利用类型变化情况[10],但其缺乏对于其它因素导致的植被覆盖度时空变异性的有效考虑：在SWAT模型中,叶面积指数(LAI)作为表征地表植被生长状况、植被覆盖度等的重要指标,是连接植被与水循环的重要桥梁[18]。但是,SWAT仅基于温度计算LAI,忽略了降水、地形等因素。这在湿润的热带、亚热带地区水文过程模拟中被更多被关注[19]。但在降水时空异质性强、地形复杂的干旱区内陆河流域,这一问题被忽视。

基于卫星遥感的长时间序列LAI数据产品因能较好展示LAI的空间异质性、季节性变化等,在评估区域植被生长状况、植被覆盖条件等方面具有重要应用。如,中分辨率成像光谱仪叶面积指数MODIS(MODerate-resolution Imaging Spectroradiometer) LAI[20]、全球地表卫星叶面积指数GLASS (Global Land Surface Satellite) LAI[21]、全球长时间序列叶面积指数GLOBMAP (Long-term Global Mapping) LAI[22]等。这些产品在不同区域的时空精度具有差异,且部分产品在中国西北部分区域存在空值区[21—23]。相对来说,GLASS LAI 产品在中国范围内具有相对较高的准确性[21,23],并且具有较好的时空连续性和完整性[21]。基于该数据集进行区域生态-水文过程模拟更为可靠。

SWAT模型仅采用简单的线性方程来表达地下径流变化与潜水面埋深变化速率之间的关系,无法真实表达流域地下水过程[24]。利用成熟的地下水模型MODFLOW取代SWAT模型的地下水模块,二者通过同参变量的传输和反馈进行耦合,既发挥SWAT模型在作物生长、地表水模拟中的优势,又结合了MODFLOW对于地下水模拟的长处。土壤和水评价工具与模块化有限拆分地下水流耦合模型(SWAT-MODFLOW)耦合模型已在美国科罗拉多州南普拉特河流域[25]、非洲林波波河流域[26]、巴音河流域[10,27]的地下水-地表水转换关系、地下水补给、地下水位预测等方面得到了应用且效果优于原始SWAT模型。

巴音河位于柴达木盆地东北缘,是盆地第四大内陆河。近年来,在人工植被恢复以及气候暖湿化背景下,巴音河流域植被覆盖条件明显变好[10]。本文为探究干旱区内陆河流域植被覆盖变化对地下水补给的影响,以巴音河中下游为例,改进SWAT模型,采用GLASS LAI数据代替其LAI计算模块,再结合流域土地利用/覆被类型变化,准确刻画区域植被覆盖度变化。将改进后的SWAT模型与MODFLOW模型耦合。在准确模拟植被覆盖变化下流域地下水循环过程的基础上,基于控制变量法分析巴音河中下游地区植被覆盖变化对地下水补给的影响。本研究可为干旱区内陆河流域生态环境保护与可持续发展提供科学依据。

1 材料与方法

1.1 研究区概况

巴音河位于青海省德令哈市,是柴达木盆地东北部最大的内陆河,也是德令哈市境内最大的河流。该河流中下游段为德令哈市区北部黑石山水库以下到可鲁克湖、尕海湖入湖口以上的流域[27],该流域中游段由北向南穿过德令哈市市区,下游段由东向经郭里木、戈壁乡最终汇入克鲁克湖-托索湖湖区(图1)。巴音河出山后河水大量下渗,河道径流量较少,部分河段河水全部潜入地下。气候寒冷、干旱,水资源短缺,日照时间长,属于高原荒漠半荒漠气候,主要的土地覆盖类型为荒漠和草地[28]。

图1 研究区概况图Fig.1 The study area

1.2 GLASS LAI数据集

本文采用时空分辨率为 8 d/250 m 的最新 GLASS LAI v5 数据,时间跨度为 2001—2019年。该数据集由北京师范大学全球变化处理与分析中心(http://www.glass.umd.edu/Download.html)发布。该数据集采用了加权线性组合的方式,首先将 MODIS LAI和基于生物地球物理参数的叶面积指数(CYCLOPES LAI) 融合,产生二者的融合 LAI[20],再与经过预处理的中分辨率成像光谱仪/高级甚高分辨率辐射计(MODIS/AVHRR) 地表反射率数据,建立广义回归神经网络训练样本,训练出地表反射率与 LAI 值的关系模型,最后将 MODIS/AVHRR 地表反射率数据作为关系模型输入数据,得到全球长时间序列的GLASS LAI 产品[29]。使用该数据之前进行了特异值去除并将有效值归为 0—10。

1.3 土地利用/覆被数据的获取

本研究基于携带专题制图仪的美国陆地卫星系列第五颗卫生Landsat5 TM、携带陆地成像仪的美国陆地卫星计划第8颗卫星Landsat8 OLI、全球环境与安全监测计划的第二颗卫星(Sentinel-2A) 影像数据获取巴音河中下游2001—2019年土地利用/覆被数据。首先在研究区针对耕地、林地、草地、水体、城镇用地、裸地六种地类开展野外调查。利用 Landsat5 TM(2011 年 11 月之前)、Landsat8 OLI(2013 年 4 月之后)、以及 Sentinel-2A(2015 年 8 月之后)影像数据,创建归一化植被指数(NDVI)、增强型植被指数(EVI)和归一化水体指数(NDWI)时间序列数据集,结合野外调查结果,实现各地类光谱特征、植被指数特征、时空变化特征等特征提取,进而利用随机森林算法对研究区五种地类进行分类。运用混淆矩阵对获得的土地利用/覆被数据进行精度评价,其总体精度达到94%,kappa系数为0.96。图2显示了本研究获取的2001年及2019年研究区土地利用/覆被数据。

图2 2001年及2019年巴音河中下游土地利用/覆被数据Fig.2 Land use/ land cover maps of the Middle and Lower Reaches of the Bayin River in year 2001 and 2019

1.4 SWAT模型改进

本研究利用GLASS LAI代替原始SWAT模型基于理想叶面积发育模型计算的LAI。过程如下：将基于GLASS数据集的250 m 分辨率、8 d LAI 数据集,映射至SWAT模型每个基本计算单元(HRU) 上。由于 LAI数据集与HRU的空间尺度不同,首先对二者进行空间叠加。叠加完成后,对于 1 个 HRU 包含多个网格的情况,对所有网格的 LAI 值进行平均后,赋予相应的 HRU。由于SWAT 模型模拟水文过程的时间步长为 1 d,故在 LAI 数据集与动态 HRU 完成对应后,对每个 HRU上的8 d LAI 数据进行插值,使其时间分辨率变为 1 d,从而实现GLASS LAI与SWAT 模型的耦合(如图3)。生物量及流域产流、产沙等后续过程均基于改进后的作物生长模块计算。由于改进后的SWAT模型更能有效刻画区域植被动态变化,故本研究将其命名为DVSWAT。

图3 GLASS LAI与SWAT -MODFLOW模型的耦合Fig.3 Coupling of the GLASS LAI and SWAT -MODFLOWHRU：水文响应单元 Hydrologic response unit；GALSS：全球地表卫星 Global land surface satellite；LAI：叶面积指数 Leaf area index；SWAT：土壤和水评价工具 Soil and water assessment tools；MODFLOW：模块化有限差分地下水流耦合模型 Modular finite difference groundwater flow model

1.5 DVSWAT -MODFLOW模型及构建

通过GIS平台建立SWAT模型的基本计算单元与MODFLOW模型的基本计算单元之间的映射关系,然后将每个HRUs分配到MODFLOW模型的相应网格单元上,将DVSWAT模型计算出的地下水补给量、蒸发量等数值用来作为地下水模拟的边界条[27]。同时,利用映射关系把MODFLOW模型所计算出的地下水排泄量添加到DVSWAT模型的相应地表水计算单元中,以此来实现DVSWAT模型与MODFLOW模型之间的耦合(图3)[30]。其中,由DVSWAT模型负责计算降水对地下水的直接补给以及灌溉对地下水的补给。由DVSWAT模型负责计算河水水位并结合MODFLOW模型的河流子模块计算河水渗漏对地下水的补给。

本文选取了德令哈站的逐日降水、温度、相对湿度、风速及太阳辐射数据作为DVSWAT模型输入数据,结合30 m分辨率的DEM数据(http://gdex.cr.usgs.gov/gdex/)、全国1∶400万土壤类型数据、2001、2010、2015、2019年土地利用类型数据以及2001—2019年逐日GLASS LAI数据建立DVSWAT模型。土壤属性数据查阅自《青海土壤》[31]及《德令哈市志》[32]。

采用研究区钻孔数据、地下水位观测数据、河网数据、数字高程模型(DEM)数据等建立MODFLOW模型,建模方法主要参考已有研究[27,33]。以巴音河中下游流域边界以及浅层含水层底板作为隔水边界,流域入口及出口作为定水头边界,河网作为河流边界。将研究区概化为250 m×250 m的规则网格,全区共251行、272列。本研究根据需要仅模拟浅层含水层地下水循环过程,且将浅层含水层概化为单层非均质各向同性,将地下水运动概化为二维非稳定流[27,33]。其中,给水度、渗透系数、河流渗透系数等参数主要参考已有研究[33],再结合研究区地下水位观测数据进行微调。

采用了基于遥感的全球陆地蒸散发阿姆斯特丹模型GLEAMv3(最新)蒸散发数据(https://www.gleam.eu/)、地下水位观测数据对DVSWAT-MODFLOW模型进行校准。其中,GLEAM数据集空间分辨率为0.25°,时间分辨率为1天。该数据集在中国区域内精度较高[34],被广泛应用于水文模型校准及验证[35—36],具有较高可靠性。此外,利用决定系数(R2)、纳什系数(NSE)、误差百分数(PBIAS)、绝对误差(AE)评价模型的适用性[10,37]。

2 结果与讨论

2.1 巴音河中下游土地覆被类型变化

基于2001年和2019年土地利用/覆被数据(图2)提取了研究区土地利用转移矩阵(表1)。从2001年至2019 年,土地覆被类型发生了明显变化。其中,草地、建设用地、林地和水体面积增加,2019年建设用地、林地和水体地面积约是2001年的5.41倍、2倍和2.35倍,草地面积增加了98.96%；裸地和耕地面积减少,分别减少了30.75%和52.88%。本研究主要关注植被覆盖增加对流域地下水补给的影响。故将2001及2019年土地利用/覆被数据进行空间叠加,提取研究区植被恢复斑块,结果如图4所示。其中,草地恢复主要在流域北部及东南部山区,林地恢复主要在流域中部农业区。

表1 2001—2019 年土地利用转移矩阵/km2Table 1 Land use transfer matrix from year 2001 to 2019

图4 研究区植被恢复情况Fig.4 Revegetation of the study area

2.2 巴音河中下游植被覆盖度时空变化

图5显示了基于GLASS LAI数据集的2001年及2019年研究区7月平均LAI空间分布情况。从空间上看,研究区北部高海拔山区、中部农业灌区、中南部地下水位高值区植被覆盖状况相对较好,这些区域2019年7月植被覆盖度明显高于2001年同期。2001年研究区7月(植物生长旺季)平均LAI值为4.37,2019年7月平均LAI值为5.63,比2001年增加了28.83%。这与巴音河流域人工植被恢复及气候暖湿化有关[27,33]。

图5 2001 年及2019 年巴音河中下游年平均LAI空间分布Fig.5 Distribution of the annual average LAI in the Middle and Lower Reaches of the Bayin River in year 2001 and 2019

图6 基于植被动态变化的土壤和水评价工具与模块化有限拆分地下水流耦合模型(DVSWAT-MODFLOW)模型月蒸散发模拟效果Fig.6 The monthly ET simulated by the DVSWAT-MODFLOWNSE：纳什系数 Nash-Sutcliffe efficiency；PBIAS：误差百分数 Percent bias；R2：拟合优度 Absolute error；ET：蒸散发 Evapotranspiration

2.3 DVSWAT-MODFLOW蒸散发模拟效果评价

巴音河流至中下游后大量渗漏,补给地下水。故河道径流量较小且无观测数据。基于河道径流数据的传统参数校准方法难以适用。本研究基于GLEAM v3蒸散发数据集校准DVSWAT-MODFLOW模型蒸散发相关参数。GLEAM v3蒸散发数据集在资料缺乏区水文模型校准中具有重要应用[33]。图6显示了DVSWAT-MODFLOW模型的月蒸散发模拟效果。各子流域R2值在0.83以上,NSE值在0.68以上,PBIAS绝对值在22%以内。可见模型对于月蒸散发的模拟效果较好[38]。图6亦显示了DVSWAT-MODFLOW模型模拟的各子流域多年平均蒸散发分布情况。流域下游子流域地下水位较高,年蒸散发量相对较大,超过200 mm。个别子流域年蒸发量超过300 mm,这是相应子流域仅包含耕地且被大水漫灌所致。

2.4 DVSWAT-MODFLOW地下水位模拟效果评价

图7 巴音河中下游月地下水位实测及模拟值对比Fig.7 The simulated and observed groundwater level in the Middle and Lower Reaches of the Bayin River1—10号观测井观测数据涵盖2019年各月；11号观测井观测数据涵盖2009—2011年各月；12—14号观测井观测数据涵盖2013—2015年各月；15号观测井观测数据涵盖2004—2005年各月；16—20号观测井观测数据涵盖2001—2014年各月

图7显示了DVSWAT-MODFLOW模型月地下水位模拟效果。20个地下水位观测井的空间分布较为均匀(图1)。流域南部观测井较为密集,月地下水位模拟效果相对优于其它区域,R2达到0.95以上,绝对误差(AE)小于0.35 m。1号及17号观测井处绝对误差较大,接近1 m。总体上来说,各观测井月地下水位模拟效果较好。

图8 不同地表覆盖类型对应的年平均地下水补给量 Fig.8 The annual average groundwater recharge corresponding to the different land cover types

2.5 植被覆盖增加对地下水补给的影响

图8显示了研究区不同地表植被覆盖类型对应的年平均地下水补给量。为剔除气象因素、坡度、土壤类型等对地下水补给的影响,本研究对具有相同气象条件、土壤类型及坡度的HRUs所对应的土地覆被类型进行了统计。在年平均降水量为244.23 mm的情况下,耕地对应的年平均蒸散发量、地下水补给量最大,分别为393.12 mm及373.46 mm这主要是由大水漫灌所致。裸地对应的年平均蒸发量为178.14 mm,年平均地下水补给量为55.68 mm。林地对应的年平均蒸发量为193.04 mm,年平均地下水补给量为39.53 mm。草地对应的年平均蒸发量为179.79 mm,年平均地下水补给量为50.14 mm。可见研究区地下水补给量耕地>裸地>草地>林地。表明巴音河中下游耕地、裸地的大量减少和草地的大面积增加使地下水补给量降低。

为探究巴音河中下游植被覆盖度增加对地下水补给量的影响,本研究选取年平均LAI增加较明显的26号子流域为典型区(图9),将2001年LAI与2019年LAI分别输入DVSWAT-MODFLOW模型。结果表明,2001年该典型区平均LAI在0—0.48之间,在年平均降水量为244.23 mm的情况下,对应的年平均蒸散发量为195.53 mm,对应的年平均地下水补给量为38.14 mm。2019年该典型区平均LAI在0—0.79之间,在年平均降水量为244.23 mm的情况下,对应的年平均蒸散发量为208.10 mm,对应的年平均地下水补给量为26.92 mm(图9)。可见植被覆盖度增加使蒸散发量增加,地下水补给量减少。

图9 不同植被覆盖度对应的地下水补给量Fig.9 The annual average groundwater recharge corresponding to the different vegetation coverage

为在全流域尺度综合探究植被覆盖增加对流域地下水补给的影响,首先基于研究区2001年对应的土地利用/覆被数据以及LAI数据建立DVSWAT-MODFLOW1。在保证它地表特征参数及气象数据等输入参数不变的前提下,基于图4所显示的植被恢复斑块建立研究区草地、林地面积增加(其它土地利用类型不变)后的土地利用/覆被数据。再结合2019年LAI数据建立DVSWAT-MODFLOW2。图10显示了DVSWAT-MODFLOW1及DVSWAT-MODFLOW2模拟的年际尺度河水渗漏补给地下水量。表明植被覆盖增加将使2001—2019年河水渗漏补给量减少。减少量最小值为2 mm,出现在2006年,对应的年降水量亦最小,为143 mm。最大减少值为10.34 mm,出现在2019年,对应的年降水量为256.5 mm。图10显示了模型模拟的地下水其它补给量(包括降水直接补给、灌溉直接补给等),可见植被覆盖增加后,地下水其它补给量减少4.1—16.18 mm。最小减少量出现在2006年。最大减少量出现在2002年,对应的年降水量为302 mm。可见年降水量可在一定程度上决定植被覆盖增加对地下水补给量的影响强弱。在月际尺度上,植被覆盖增加将使河水渗漏补给量减少0—7.18 mm(图11)。减少量为0的月份一般处于较为寒冷、干旱的植物休眠期或河流枯水期。河水渗漏补给最大减少量出现在2012年7月,对应的月降水量为109 mm。此外,植被覆盖增加将使地下水其它补给量减少0—7.85 mm(图11)。减少量为0的月份同样处于较为寒冷、干旱的植物休眠期或其它降水较少时期。最大减少量出现在2016年8月,对应的月降水量为113 mm。植被覆盖增加对地下水补给量的影响较强的月份集中于植物生长旺盛且降水较多时期。综合来看,2019年植被覆盖情况对应的年际及月际尺度地下水补给量较2001年分别减少了6.1—26.52 mm以及0—15.03 mm。

图10 植被覆盖增加对年地下水补给的影响Fig.10 The impacts of the increasing vegetation coverage on annual groundwater rechargeDVSWAT：耦合植被覆盖变化的SWAT模型 Soil and water assessment tool coupled with dynamic vegetation coverage

图11 植被覆盖增加对月地下水补给的影响Fig.11 The impacts of the increasing vegetation coverage on monthly groundwater recharge

3 结论

(1)DVSWAT-MODFLOW模型的月蒸散发模拟效果较好,各子流域R2值在0.83以上,NSE值在0.68以上,PBIAS绝对值在22%以内。月地下水位模拟效果亦较好,各观测井R2值在0.92以上,绝对误差在0.28—0.97 m之间。

(2)2019年巴音河中下游林地及草地面积较2001年增加明显,分别增加了5.41倍及98.96%。耕地及裸地面积减少明显。此外,2019年植物生长旺季平均LAI值为5.63,比2001年增加了28.83%。

(3)年际尺度上,巴音河中下游植被覆盖增加导致河水渗漏补给地下水量减少2—10.34 mm,导致地下水其它补给量减少4.1—16.18 mm。植被覆盖增加对年际尺度地下水补给量的影响强弱在一定程度上取决于年降水量。

(4)月际尺度上,巴音河中下游植被覆盖增加导致河水渗漏补给地下水量减少0—7.18 mm,地下水其它补给量减少0—7.85 mm。植被覆盖增加对月际尺度地下水补给量影响较强的月份集中于植物生长旺盛且降水较多时期。