基于地表水与地下水分割校正的漓江流域水供给服务时空格局研究

刘 佳,肖 玉,*,张昌顺,黄孟冬

1 中国科学院地理科学与资源研究所, 北京 100101

2 中国科学院大学, 北京 100049

生态系统服务指人类直接或间接地从生态系统过程和属性中获得的福祉[1],服务的供给能力直接关系着社会经济系统是否能持续运行[2]。水供给服务是非常重要的生态系统供给服务[3],一方面直接为人类生产生活提供水产品[4],另一方面为多种生态系统过程提供媒介,间接影响食物供给、气候调节、景观美学等多种生态系统服务[5]。然而,受气候、土地利用、水资源需求和养分循环变化等多重压力因素的持续影响[6—7],水供给服务表现出明显的空间异质性和时间可变性。这些压力因素将在未来几十年内持续增加[8],水供给服务的稳定性面临重大考验。因此,精确评估不同区域不同条件下的水供给服务,根据实际情况校正水供给服务的空间分布,有助于提高区域水资源的管理水平,为未来水源地空间规划和生态服务功能提升提供科学依据。

产水量是衡量水供给服务能力的关键要素,其时空分布状况直接影响着水供给服务的有效性和可获得性[9]。定量评估产水量及其时空变化特征是水供给服务研究的核心与关键[10],GIS分析方法、水文物理模型、机器学习的发展为水供给服务定量化与空间化研究提供了坚实的技术支撑[11—13]。目前,InVEST模型[14]、SWAT模型[15]、半分布式HSPF水文模型[16]等经常被用于评估特定区域的产水量[17—19]。其中,InVEST模型凭借运行参数少且空间表达清晰的优势已成为计算产水量的主流方法[20],其结果广泛用于探究西北干旱区[21]、东北黑土区[22]、青藏高原三江源区[23]、黄土丘陵沟壑区[24—26]、秦岭地区[27]、东亚季风区[28]、云南复杂山区[29]等区域水供给服务的时空变化规律。然而,针对特殊地质背景的水供给服务研究还比较少,尤其是关于喀斯特地区的研究更为少见。

作为全球典型的脆弱生态系统之一,喀斯特生态系统敏感度高且抗干扰能力弱,不合理的人类活动导致该区土壤流失不断加剧,石漠化情况日益严重。明晰喀斯特地区水供给服务的空间分布可为当地生态治理与恢复提供科学依据。近年来,逐渐有学者关注到喀斯特地区水供给服务并将InVEST模型应用于评估该地区水供给服务能力[30—31]。然而,这些研究很少考虑喀斯特特殊的地貌因素对水资源赋存和分布规律的影响[32]。Zhang等[33]关注到这个现象并基于InVEST模型探讨了喀斯特景观和非喀斯特景观的水供给服务能力差异,结果显示产水量在喀斯特区域和非喀斯特区域中没有明显的差异。事实上,喀斯特与非喀斯特因下垫面条件不同,其产流汇流特征、水源分配等方面存在不同程度的差别[34],与水相关的生态系统服务能力也因此有所差异。如夏林等[35]比较贵州省内乌江流域不同喀斯特景观水分涵养能力发现,非喀斯特、亚喀斯特和纯喀斯特景观水源涵养能力逐渐降低。可见,喀斯特地区生态系统服务水分涵养能力不仅低于非喀斯特地区,而且不同喀斯特类型的生态系统服务水分涵养能力也有巨大差异。然而,到目前为止,尚未有研究能较好地反映水供给服务能力在不同喀斯特类型间的差异。相关研究主要受两方面因素限制:一是模型忽略了喀斯特生态系统的特殊性。由于InVEST模型多采用全国性的普遍参数来估算产水量,模拟结果不区分地表径流、壤中流和地下径流,可能造成喀斯特地区产水量的模拟值高于实际值,进而无法准确评估喀斯特地区生态系统服务水供给服务能力。Zhang等[33]也认识到这个问题,但是没有深入探讨如何改进模拟结果使其更接近于喀斯特地区水供给服务的真实值。二是缺乏高精度岩溶个体及其组合形态的地貌数据。虽然全国地貌区划数据库对岩溶分布区基本地貌进行了划分,但划分结果难以体现岩溶发育类型对区域水资源分配的影响程度[36],不能很好地衡量水供给服务在不同喀斯特地貌间的差异。因此,有必要提取岩溶地貌单元数据并结合喀斯特生态系统结构与功能的特点[37],校正InVEST模型在喀斯特地区的模拟结果[38],评估不同区域不同地貌单元水供给服务能力。

喀斯特地貌又称岩溶地貌,具有特殊的二元三维结构,发育有地表与地下两种喀斯特形态[39],其水循环和水资源过程与同一气候区的非喀斯特有极大的差异[40]。同非喀斯特相比,喀斯特地下水对降雨响应迅速[41],丰沛的降雨难以在喀斯特地表汇集形成河流,大气降水和地表水会沿着岩石裂缝、漏斗和天坑等深入地下转化为地下水,并沿着纵横交错的洞穴、暗河移动至其他区域[32]。因此,喀斯特地区大部分流域拥有独特的地表、地下双层径流系统[42],地表成为喀斯特径流的形成场和分配场,地下则是喀斯特径流的输移场和调蓄场[43],区域产水量也随之发生明显的纵向分异。故而,改进喀斯特地区InVEST模型模拟结果面临的第一个问题是分割区域产水量。大气降水入渗系数是岩溶地下水系统最基础的水文地质参数,它反映了水资源在地表与地下的分配方式[44]。利用降雨入渗系数可以简单有效地区分不同岩溶地貌单元地表与地下产水量[45—46]。地表产水量是存在于地表之上包括地表径流在内水量,相对容易在原位或小流域内被利用;而地下产水量是储存于岩溶裂隙、溶洞、暗河中水量,可能会跨流域流动,很难在原位或小流域内被利用。值得注意的是,喀斯特地区水体最终都会以地表径流或地下河的形式汇入最低排泄基准面[47]。换言之,除去人工取水,地下产水量最终会因注入最低排泄面而被人类使用。因此,第二个问题是根据水资源的可用性对地下产水量进行空间校正。目前水文学常用概念模型、物理模型、数值模型和分布式岩溶水文模型来刻画喀斯特复杂的地下径流过程[48]。其中,概念模型将地下水视为一个整体,以地下水补给、径流和排泄过程的物理概念为基础,将产汇流过程抽象化后通过数学方法进行定量模拟,因其方法简单且所需参数较少而得到广泛应用[49]。鉴于此,本研究创造性地提出地下产水量分割系数并将地下水补给、径流和排泄过程概念化,结合水量平衡法构建地表水与地下水分割校正概念模型,以期校正喀斯特地区InVEST模型模拟结果。

漓江流域是喀斯特地貌的典型代表之一,峰丛、峰林、洼地等岩溶地貌单元的发育极大地影响了流域产流与汇流过程。分析喀斯特地区产水量时应遵循因地制宜原则,从地貌尺度出发探究水供给服务在不同地貌分区的差异。此外,漓江流域作为一个较为完善的水文地质单元,从子流域角度评估水供给服务有助于提高当地水资源管理水平。因此,本文基于DEM数据构建30 m精度岩溶地貌单元数据,利用ArcSWAT模型划分漓江子流域,分别采用InVEST模型和喀斯特产水量分割校正概念模型估算2000、2005、2010、2015和2020年漓江流域产水量,分析漓江流域水供给服务在栅格、地貌分区和子流域尺度的时空分布格局,探究不同地貌分区、不同子流域水供给服务能力在校正前与校正后的差异性。校正结果可为喀斯特地区水资源管理、生态补偿政策制订提供良好的数据基础。

1 研究区概况

漓江流域(24°6′39″N—25°54′56″N,110°4′59″E—111°17′35″E)坐落于桂东北山间盆地谷地岩溶地貌区。区内地势大致呈四周高中部低(图1),边缘地区多环绕有中低山和丘陵,中部以岩溶地貌为主,漓江谷地自北向南纵贯全区。中部岩溶区根据含水岩组的埋藏条件划分为裸露型、覆盖型和埋藏型(图1)。裸露型岩溶区以纯碳酸盐岩为主,地表石峰丛聚,洼地广布,漏斗、落水洞等垂直岩溶形态发育强烈。覆盖型岩溶区表层覆盖有第四系土层,大部分岩溶表面上离散分布着各种塔状、柱状或密集成林状的峰林。埋藏型含水岩组多分布于可溶性基岩之下。岩溶发育地区,地表河流稀少,但各种各样的地下洞穴纵横交错,为岩溶地下水的形成和赋存创造了良好的条件。流域属于典型的中亚热带季风性湿润气候,夏长冬短,雨热同期,年均气温为18—21℃,年均降水量约1688 mm,全年雨量充沛但季节分布不均。境域土地利用类型以林地为主,其次是耕地和草地(图1)。林地主要分布在北部和东部地区,耕地集中于河流两侧平原地带。城镇及工矿用地主要集中在兴安县、灵川县、桂林市区、阳朔县、平乐县等城市驻地。其核心干流漓江发源于越城岭山脉的最高峰猫儿山,流经灵川、桂林、阳朔至平乐,与荔浦河、恭城河相汇成为桂江上游,而后继续南流至梧州汇入西江。漓江地表径流来源于流域内的地表水和地下水,在雨洪时地表水向地下水渗透,低水和枯水期地下水补给河槽[50],地下水流向整体由西北流向东南,最终排泄至漓江河谷,成为漓江径流的重要组成部分。

图1 研究区DEM、岩溶类型和2020年土地覆被类型分布图Fig.1 Distribution of DEM, karst type and land cover type in 2020 in the study area

2 数据来源与处理

2.1 数据来源

本研究涉及的数据包括DEM、土壤、土地覆被、气象日值、水文、河流、行政边界、岩溶分布等数据。具体来源和描述见表1。所有数据的投影坐标统一为Krasovsky_1940_Albers坐标,栅格数据空间分辨率:30 m。

表1 研究数据来源与描述Table 1 Data sources and descriptions

2.2 数据处理

2.2.1ArcSWAT划分子流域

基于30 m分辨率的SRTM DEM数据,利用ArcSWAT的Automatic Watershed Delineation模块在漓江流域划分出142个小流域,而后结合流域内的自然地理特征和水系分布对小流域进行合并,最终得到18个子流域,分别用数字1—18来标记(图2)。18个子流域根据其位置被分为五个区域:漓江上游区、漓江中游区、漓江下游区、荔浦河区和恭城河区。其中,1—3号子流域位于漓江上游区,4—8号子流域位于漓江中游区,9—12号子流域位于漓江下游区,13—14号子流域位于荔浦河区,15—18号子流域位于恭城河区。

2.2.2基于DEM划分地貌形态

本研究以30 m DEM为基础信息源,结合数字地形分析方法实现了岩溶地貌峰丛区、洼地区和平原区边界的准确提取[51]。首先,基于正负地形分析方法将复杂多样的岩溶地貌降维简化成正负地形[36],有效突出岩溶地貌形态差异,采用趋势面分析方法构建平原趋势面,确定流域峰丛与平原的边界。其次,采用GIS水文分析方法提取流域山脊线与山谷线,根据山脊线和山谷线的交点确定地形特征点并通过地形坡度变率筛选出地形鞍部点[52]。最后,基于峰丛洼地的分型维特性采用反距离权重插值方法构建鞍部趋势面,进而确定洼地的边界[53]。岩溶地貌类型提取结果如图所示(图2),流域内岩溶地貌呈南北向展布,从横向上来看,自两侧向中部体现非岩溶→峰丛→洼地→平原的展布序列;从纵向来看,峰丛洼地、峰林平原并非只发育在某个岩溶区之中,而是不均衡地发育在裸露型、覆盖型和埋藏型岩溶区之上。从地貌组合类型的分布位置来看,峰丛洼地分布于盆地谷地和非岩溶斜坡地带,小部分呈岛状散布在峰林平原中。峰林平原平行发育于峰丛洼地附近,连片展布于漓江两侧,少部分平原呈条带状穿插于峰丛洼地中。

图2 子流域和地貌形态空间分布图Fig.2 Spatial distribution of sub-watersheds and landform types漓江流域划分出142个小流域,而后结合流域内的自然地理特征和水系分布对小流域进行合并,最终得到18个子流域,分别用数字1—18来标记

2.2.3地下水补给、径流和排泄边界

研究区可分为河流阶地、岩溶区和非岩溶区,岩溶区根据含水岩组的出露条件又可分为裸露型、覆盖型和埋藏型。考虑到喀斯特地貌的形成与演化多聚焦于水,不同的地貌形态可以反映出地下水的补给、径流和排泄特征[54],本文基于上述分类结合地表岩溶地貌形态,将研究区细化为11种地貌单元,并根据地貌类型与地下水补径排特征确定了地下水补给区、径流区和排泄区边界。地下水主要补给来源为降雨入渗,补给区涵盖整个流域。降落到非岩溶地区的水大部分以山前径流、地表河流等外源水形式进入岩溶区,仅有一小部分入渗形成基岩裂隙水,侧向移动补给相邻岩溶区。岩溶区地下水获得补充后沿着洞穴、地下河网开始流动,其径流方向受到地形倾向、构造分布方向和地表水系展布的制约[55],大体上自四周碎屑岩山区、峰丛洼地向中部平原流动,最终汇入漓江河谷进行集中排泄,成为漓江地表径流的重要组成部分[56]。由于峰林洼地地表之上没有明显的排水系统,大多数降水通过落水洞、竖井、漏斗等灌入地下,然后在平原上或裸露的峰林山脚下以泉或地下河的形式排泄,因此将平原划分为主要径流区[41],峰丛洼地为重要补给区。径流区汇集流域地下水补给后,最终都会汇入最低排泄基准面即漓江进行集中排泄,漓江及其支流恭城河、良丰河、潮田河等是其主要排泄通道,据此划定平原区的漓江及其支流作为集中排泄区。具体分类见表2。

表2 地貌单元及地下水区域Table 2 Geomorphological units and groundwater areas

3 研究方法

本研究将产水量视为一个整体,利用喀斯特地貌类型特征和降雨入渗系数确定产水量分割系数,依据分割系数将InVEST模型产水量分为地表产水量和地下产水量两部分,然后结合地下水补径排特征和水量平衡法建立产水量空间校正概念模型,并模拟漓江流域2000—2020年水供给服务时空分布特征。产水量空间校正模型示意图见图3。

图3 喀斯特地区产水量空间校正模型示意图Fig.3 Schematic diagram of the spatial correction model of water yield in karst areas

3.1 InVEST模型

InVEST模型Water yield模块以水量平衡原理为基础,结合地形、气候、植被、土壤等因素,在统一的栅格尺度上用降水量(输入)减去实际蒸散发量(输出)来估算该栅格的产水量[20]。产水量结果不区分地表径流、壤中径流和基流,涵盖所有径流、土壤含水量、枯落物持水量和冠层截留量。模型算法如下:

PETxj=Kcxj×ET0xj

AWCx=min(MaxSoildepth,Rootdepth)×PAWC

PAWC=54.509-0.132×sand%-0.003×(sand%)2-0.055×silt%-0.006×(silt%)2-0.738×clay%+0.007×(clay%)2-2.688×OM%+0.501×(OM%)2

式中,Yxj是第j种土地利用类型像元x的年平均产水量(mm);Px是像元x的年均降水量(mm);AETxj是第j种土地利用类型像元x的年实际蒸散发量(mm);PETxj为第j种土地利用类型像元x的潜在蒸散量(mm);ωx是一个经验参数;ET0xj是第j种土地利用类型像元x的参考蒸散量(mm);Kcxj是第j种土地利用类型像元x的植物蒸散系数;Z是经验常数;AWCx是土壤有效含水量(mm);MaxSoildepth是土壤的最大根系埋藏深度;Rootdepth是植物根系深度;PAWC是植物可利用水含量;sand%指土壤质地中砂粒比重;silt%土壤质地中粉砂比重;clay%土壤质地中粘粒比重;OM%土壤质地中有机质比重。

3.2 产水量分割系数

产水量分割系数指单位时间单位栅格面积地下产水量与总产水量的比值,其值取决于区域降雨入渗系数。基于前人对喀斯特地区降雨入渗系数的研究成果,确定不同地貌类型产水量分割系数。非岩溶区以碎屑岩为主,产水量分割系数为0.15—0.2。裸露型岩溶区地下水以灌式补给为主,分割系数为0.5—0.6。覆盖型岩溶区地下水以面状入渗补给为主,分割系数为0.2—0.4。埋藏型岩溶区地下水以其他相邻含水层补给为主,分割系数为0.2—0.3。由于同一岩溶区内地表岩溶形态各异,涵盖峰丛、洼地、平原谷地和河流阶地4个地貌单元,降雨下渗系数表现为河流阶地<平原<峰丛<洼地[57]。据此进一步对分割系数进行调整以提高产水量分割精度。对于没有研究结果的区域,通过类比移植到条件类似的地貌单元,最终得到11个地貌单元的产水量分割系数,见表3。

表3 产水量分割系数Table 3 Water yield partitioning factor

3.3 产水量空间校正模型

根据分割系数将单元栅格产水量划分为地表产水量和地下产水量,然后依据地下水补径排特征构建补给区、径流区和排泄区的水量平衡方程式,完成产水量的空间校正。具体公式如下:

GYxi=αi×Yxj

SYxi=(1-αi)×Yxj

AYxi=SYxi

BYxi=SYxi

式中,αi是地貌单元i的产水量分割系数;GYxi是地貌单元i上像元x的地下产水量(mm);SYxi是地貌单元i上像元x的地表产水量(mm);AYxi是补给区中地貌单元i上像元x校正校正后的产水量(mm),其中i=11,12,13,21,22,23,31,32,33,5;BYxi是径流区中平原地貌单元i上像元x校正后的产水量(mm),其中i=13,23,33;CYx是排泄区中像元x校正校正后的产水量(mm);Areai是地貌单元i的面积(m2),其中i=11,12,13,21,22,23,31,32,33,5;Area4是地貌单元4即平原区河流的面积(m2)。

3.4 校正模型验证

基于水文站的实测数据、InVEST模型模拟结果和校正模型模拟结果,计算不同年份校正前后的年产水深度与实测年径流深度的均方根误差(Root mean squared error, RMSE),RMSE越小代表实测值与模拟值之间的偏差越小。计算方法如下:

式中,RMSE是均方根误差;N是水文站个数,St是实测值,Pt是模拟值。

4 结果与分析

4.1 InVEST模型产水量时空分布

4.1.1基于栅格的InVEST模型产水量时空分布及模型验证

InVEST模型模拟结果显示,2000年、2005年、2010年、2015年和2020年平均产水深度分别为1069.11 mm、1024.98 mm、859.95 mm、1512.55 mm、1174.01 mm,多年平均产水深度为1128.12 mm(图4)。产水深度年际波动较大且不同时期的变化趋势有明显差异:2000—2010年流域产水深度呈轻微下降趋势,年降速为20.92 mm/a;2010—2015年上升趋势较明显,年增速为130.52 mm/a;2015—2020年又出现明显下降趋势,下降速率为67.71 mm/a。从水供给总量来看,2000—2020年漓江流域多年平均供水量为146.81×108m3,年际供水量在111.91×108m3—196.94×108m3之间,呈现先减后增再减的变化趋势,2015年供给总量最高,2010年最低。从水量平衡的角度来看,降水和实际蒸散发是决定产水深度和供水量的两个关键因素。从图4可以看出,同期平均降水量也呈现先减后增再减的趋势,基本维持在1400—2100 mm之间,2015年平均降水量最多,2010年降水量最少,极差为639.55 mm。平均实际蒸散发量整体呈下降趋势,2005年、2010年、2015年和2020年均实际蒸散发量在2000年的水平上分别下降了1.41%、2.27%、4.30%、5.87%。

图4 漓江流域降水量、实际蒸散发量、产水深度和供水总量年际变化Fig.4 Interannual variation of precipitation, actual annual evapotranspiration, water yield, and total water supply in the Li River Basin

从空间分布来看,漓江流域产水量大致呈“北高南低”的空间分布格局,产水深度从西北向东南逐渐递减,梯度变化明显(图5)。流域北部产水深度普遍高于南部是因为北部的猫儿山阻挡了南来气流,迫使气团抬升形成暴雨,使得北部的水分输入远高于南部。对于北部地区,又以西北尤其是桂林市区和兴安县等区域产水能力相对较强,主要因为这些区域地处平原地带,建设用地和农业耕地较为集中,蒸散量较低,在降水相同的情况下水分输出较少。

图5 2000—2020年基于InVEST模型产水深度空间分布图Fig.5 Spatial distribution of water yield based on InVEST model from 2000 to 2020

InVEST模型输出结果包括单位栅格年产水深度和流域年总产水量的预测值,不同的Z系数对应不同的产水量预测值[17—19]。Z系数是一个表征地区降水和水文地质特征的经验常数,数值调节范围为1—30。在其他参数确定的情况下,通过调节Z系数来校验模型模拟结果,与流域水文站实测径流量误差最小的产水量预测值对应的Z值即为模型最优系数。本研究使用流域出水口平乐站的多年平均径流量作为InVEST模型产水量结果的参考。结合2006—2019年水文数据与以往文献资源[62—63],计算出1990—2019年平乐站的多年平均径流量约为1.46×1010m3,当Z取3.5时,模拟产水量相对误差为0.31%,模拟效果较优。此外,荣检[64]基于InVEST模型评估了2000、2005、2010和2015年广西西江流域的产水功能,其中漓江流域所在地区产水深度在900—1600 mm之间变化。Wang等[65]基于InVEST模型计算了西南喀斯特地区产水量,其中桂西北地区多年产水深度高于900 mm。徐洁等[20]基于InVEST模型对1995—2010年东江湖流域的平均产水深度进行了分析,在1100—1600 mm之间变化。东江湖流域与漓江流域同属亚热带湿润季风气候,二者降水量相近但漓江流域实际蒸散发高于东江湖流域。本研究结果与荣检和Wang等的漓江区域结果较为一致,与徐洁等结果具有可比性,证明本文InVEST模型模拟的产水量较为可信。

4.1.2基于地貌单元的InVEST模型产水量

地貌是喀斯特生态系统得以存在和发展的物质基础,不同地貌组合方式决定了该区域生态系统服务的供给和维持[66—67]。研究结果显示,非岩溶区单位栅格多年平均产水深度略高于岩溶区,但仅高57.05 mm(表4)。可见,岩溶区和非岩溶区的产水能力在InVEST模型中没有明显差异,这与Zhang等[33]的研究结论一致。深入探究岩溶区产水量差异发现,随着不同岩溶含水岩层出露情况的变化,裸露型、覆盖型、埋藏型岩溶区的平均产水深度依次降低,说明裸露型岩溶区产水能力要高于覆盖型和埋藏型。随着地貌形态的变化,平原、非岩溶、河流阶地、洼地、峰丛的平均产水深度依次降低,表明平原产水能力最高,峰丛产水能力最低。具体到细化后的地貌单元,裸露型—平原区平均产水深度最高,覆盖型—平原区次之,埋藏型—洼地区最低。从供水总量来看,非岩溶区供水总量最大,多年平均供水总量高达80.26×108m3,其次是裸露型—平原区(28.37×108m3)、覆盖型—平原区(11.01×108m3)和裸露型—洼地区(11.06×108m3),覆盖型—峰丛区、埋藏型—峰区供水总量最低。供水总量可以从数量上反映区域的供水能力,但该结果受岩溶区面积影响较大。

表4 不同地貌单元的InVEST模型产水量结果Table 4 Results of water yield simulated with the InVEST model of different geomorphic units

4.1.3基于子流域的InVEST模型产水量时空分布

结果显示,子流域多年平均产水深度为1076.42 mm,其中漓江上游区(1343.66 mm)>漓江中游区(1207.70 mm)>恭城河区(982.52 mm)>荔浦河区(962.10 mm)>漓江下游区(862.97 mm)。可见,流域产水能力沿上游至下游逐渐减弱。具体到各子流域产水深度,4个子流域属于高值区(>1300 mm),6个属于中值区(1000—1300 mm),8个属于低值区(<1000 mm)。高值区包括1—4号子流域,其中4号产水能力又高于1—3号(图6)。这是因为4号子流域海拔较低,耕地和林地并重,而1—3号子流域海拔较高且以林地为主,林地植被蒸腾作用高于耕地,造成1—3号子流域水分输出较多,产水深度相对较小。中值区涵盖5—8号和15—16号子流域。剩余子流域为低值区,其中又以11—12号产水深度最低。从供水总量来看,流域分区供水总量排序为:漓江上游区>漓江中游区>恭城河区>荔浦河区>漓江下游区。这是因为汇水面积不同导致供水总量产生较大的空间变程,漓江中游区、恭城河区、漓江上游区、荔浦河区和漓江下游区的面积依次递减,供水总量也依次递减。其中,1号、16号子流域的供水总量以绝对优势占据前两位,5号子流域居第三位,多年平均供水量是1号子流域的近3/4。其后是3号、14号、2号子流域,供水总量最低的是10号、11号子流域,每年供水总量均低于1×108m3,不足1号子流域供水量的1/10。

图6 2000—2020年子流域InVEST模型产水深度空间分布Fig.6 Spatial distribution of water yield simulated with the InVEST model of different the sub-watersheds from 2000 to 2020

4.2 基于分割系数校正后的产水量时空分布及其变化

4.2.1基于栅格的校正后产水量时空分布及模型验证

本研究基于分割系数将InVEST产水量分为地表产水量与地下产水量。结果显示,2000—2020年漓江流域地表平均供水量为1.09×1010m3,约占总供水量的74.13%;另外25.87%为地下产水量,常年储存于地下。这些地表产水量与地下产水量分别沿不同的路径汇入漓江,最终到达流域出口。漓江是整个流域的最低排泄基准面,也是境内地表水和地下水的最终排泄通道[68],流域地下产水量只在流域内部进行重新分配,不发生跨境流动。因此,基于地表与地下产水量分割系数的产水量校正只对产水空间分布进行了校正,流域所有年份平均产水深度和供水总量并未发生改变。从空间分布来看,校正后的漓江流域产水量总体上表现为“西北部高于东南部,四周高于中部”的空间分布格局(图7)。高值区主要分布在非岩溶区,以流域北部、东南部和西南部的高海拔地区最为集中。原因在于该区域对雨水汇集作用较强,大部分降雨可以迅速转化为河流进而得以保存于地表。低值区分布在流域中部岩溶地区,尤其是漓江桂林—阳朔段两侧。该段以全岩溶地貌为主,流水岩溶作用强烈,水资源呈现地表水贫乏、地下水丰富的分配格局[69]。可见,校正后的产水量分布格局与水资源分配格局具有较好地一致性。

与InVEST模型模拟结果相比,2000年、2005年、2010年、2015年和2020年漓江流域每年均有96.63%区域产水深度在校正后下降,仅3.37%区域产水深度高于校正前。产水深度下降幅度可以反映产水量分割程度,即降幅越大,地下产水量越多。每年产水深度下降幅度及其时空分布具有明显差异(图8)。2000—2010年,产水深度降幅低于200 mm区域所占面积最大,其比重从48.40%稳步增长至56.24%,覆盖面积逐渐扩散至整个非岩溶地区。即2000—2010年研究区有近一半区域的地下水下渗量不超过200 mm,地下产水量呈减少趋势。2010—2020年,降幅200—500 mm的区域面积比重上升,于2015年达到峰值(61.23%)。表明地下水下渗量在200—500 mm的区域面积逐渐增大,研究区地下产水量逐渐增加。从空间分布来看,产水深度降幅程度表现出明显的圈层结构,以桂林雁山区和阳朔县交界地带为轴线,产水深度降幅量呈带状向南北两侧逐渐减少。

图8 2000—2020年研究区校正前后产水深度变化量空间分布Fig.8 Spatial distribution of water yield changes before and after correction in the study area from 2000 to 2020

本文在研究区不同位置选取了7个集水面积不等的水文站,基于实测数据和模拟结果分别计算2005、2010和2015年校正前后的年产水深度与实测年径流深度的均方根误差(RMSE)。结果显示,2005、2010和2015年校正后7个站点的年产水深度RMSE均小于校正前的RMSE(表5),表明校正后模拟值与实测值之间的偏差更小。

表5 地表与地下水分割校正后产水量结果验证Table 5 Validation of water yield results after surface and groundwater splitting correction

4.2.2基于地貌单元的校正后产水量

岩溶区是流域重要的地下水补给区,校正后的地貌单元产水深度更多反映其地表产水能力。从岩溶地貌和非岩溶地貌对比来看,非岩溶区单位栅格多年平均校正后产水深度比岩溶区高352.27 mm,表明非岩溶地区地表蓄水能力强于岩溶区。与InVEST模型结果相比,这个结论与实际情况更为相符[61]。从不同岩溶地貌单元产水量来看(表6),随着岩溶含水岩层出露情况变化,裸露型、埋藏型、覆盖型岩溶区的产水深度依次升高,反映出裸露型岩溶区地表产水能力相对最高,相同面积情况下该岩溶区地下供水总量低于其他岩溶区。随地貌形态变化,平原、峰丛、洼地校正后的产水深度依次降低,代表平原、峰丛和洼地的地表产水能力逐渐降低,但单位面积地下供水总量逐渐增多。具体到细化后的地貌单元,河流沿岸平均产水深度最高,非岩溶区、埋藏型—平原区、覆盖型—平原区地表产水深度次之,裸露型—峰丛区、覆盖型—洼地区地表产水深度最低。

表6 校正后基于地貌单元的产水量结果Table 6 Corrected water yield results of different geomorphic units

结合校正前后的地貌单元产水深度差异发现(图8),产水深度降幅最大的轴线地带是裸露型—洼地区的主要分布区,从轴线向南北方向延伸分布有裸露型—峰丛区和裸露型—平原区。作为典型的岩溶地貌区,该区域通常因大量降水下渗而导致地表干旱缺水,域内流域都存在不同程度的水土流失[40],也印证了校正后的产水深度更能反映不同地貌的实际供水能力。

4.2.3基于子流域的校正后产水量时空分布

结果表明,不同流域分区校正后供水总量排序为:漓江中游区>漓江上游区>恭城河区>荔浦河区>漓江下游区。多年平均校正后产水深度排序为:漓江中游区(1308.61 mm)>漓江上游区(1254.74 mm)>恭城河区(974.65 mm)>荔浦河区(941.37 mm)>漓江下游区(929.06 mm),说明五大流域分区中,漓江中游区产水能力最高,供水总量也最大。相应的,漓江下游区产水能力最弱,供水总量最小。具体到各子流域(图9),校正后有3个子流域产水深度属于高值区(>1300 mm),9个属于中值区(1000—1300 mm),6个属于低值区(<1000 mm),说明研究区有2/3子流域拥有较强的产水能力,每年平均产水量可以稳定在1000 mm以上。

图9 校正后子流域研究期内产水深度空间分布Fig.9 Spatial distribution of water production depth during the study period in the corrected sub-watershed

对比校正前后子流域产水深度发现(图10),漓江上游区、荔浦河区、恭城河区的子流域校正后的产水深度低于InVEST模型模拟结果,其中又以3号、12号、17号和1号子流域产水深度降幅较大。这是因为流域分区平均海拔较高导致地下水排泄区分布相对少,区域地下产水量以跨区补给为主。相对而言,漓江中游区和漓江下游区校正后的产水深度高于InVEST模型模拟值,其中校正后产水深度涨幅较大的子流域是8号、4号、10号和9号子流域,主要分布在桂林—阳朔—平乐河段。由前文分析可知,桂林—阳朔河段产水深度在校正前后降幅最大,导致校正后地表产水能力相对较弱。但该区地处流域谷地腹中,兼备地下水补给区与排泄区,区域水循环交替频繁且强烈[70],校正后不但获得了本区地下产水量的补给,还获得了来自其他地区的补给,因此,综合后的产水深度高于校正前。值得注意的是,这也是导致校正后漓江中游区产水能力高于上游区产水能力的重要原因。

图10 校正前后子流域研究期内产水深度变化Fig.10 Didfference of water yield in the sub-basin before and after correction

5 讨论

西南岩溶流域是我国长江水系和珠江水系的上游,也是我国典型的生态脆弱区,明晰该地区真实的水供给服务能力对生态恢复和重建具有重要意义。岩溶发育地区,地貌与水文相互依存又彼此制约。一方面,水是岩溶地貌演化的直接动力,塑造了峰丛洼地、峰林平原、岩溶峡谷、岩溶断陷盆地等多种岩溶地貌,特殊的地质地貌分异结构使得喀斯特地区生态系统服务具有明显的时空异质性[71]。另一方面,不同的地貌组合形成了特殊的产汇流过程,演化出独特的地表、地下双重结构[43],生态系统服务也随之发生明显纵向分异。然而,目前在喀斯特地区使用InVEST模型进行水供给服务评估时存在一定的局限性:一是因为InVEST模型简化汇流过程,没有区分地表径流、壤中流和基流,未能较为准确地反映岩溶区水供给服务真实情况;二是因为缺乏高精度岩溶个体及其组合形态的地貌数据,研究结果难以揭示不同地貌类型水供给服务能力。

针对以上限制因素,本研究探索性地提出地下产水量分割系数,将地下水补径排特征概念化,以水量平衡法为基础地构建了地表水与地下水分割校正概念模型。同时,以漓江流域为例,基于DEM数据构建了30 m精度的喀斯特地貌单元数据,从栅格、地貌分区和子流域尺度对比分析了InVEST模型和分割校正模型的水供给服务空间格局。结果表明,校正后的产水量与流域实测值更为接近,其空间分布情况与实际更为相符。就岩溶区和非岩溶区水供给能力而言,校正前岩溶区与非岩溶地区的水供给服务能力没有明显差异,校正后的非岩溶地区水供给服务能力明显强于岩溶区。事实上,岩溶发育地区,地表河流稀少但地下水资源十分丰富,而InVEST模型只考虑降水(水分输入)和实际蒸散发(水分输出),忽略了地下水下渗过程,一定程度上掩盖了岩溶区与非岩溶区的产水量差异。本文的分割校正模型在此基础上充分考虑了不同地貌类型的地下水补径排特征,校正后的水供给服务能力更贴近实际。就不同岩溶区水供给能力而言,校正前裸露型岩溶区水供给能力最高,校正后其水供给能力降为最低。变化原因在于裸露型岩溶区地表没有明显排水系统,大多数降水通过落水洞、竖井、漏斗等灌入地下形成丰富的地下产水量,而地下产水量又会迅速沿着地下河网移动至平原上或裸露的峰林山脚下并以泉或地下河的形式注入漓江,最终导致裸露型岩溶供水能力在校正后发生大幅度下降。就子流域水供给能力而言,校正后漓江中游区产水能力提高至第一,而上游区产水能力降为第二。这是因为上游区域属于非岩溶森林生态分区,几乎没有获得地下水补给。而中游地区包含漓江流域最低排泄面,在校正后获得了大量地下水补给,产水能力因此大幅提升。综上,本研究提出的分割校正模型较为准确地评估喀斯特地区水供给服务空间格局,有效提高了InVEST模型在喀斯特地区的适用性,研究结果不但可为后续水供给服务流动模拟研究与生态补偿机制建立提供更科学合理的数据基础,还有助于当地政府制定科学合理的水资源利用政策。

然而,本研究仍存在一些局限,需要在未来做更深入的研究。首先,本文采用的入渗系数主要参考相似地区的前人研究成果,可能会影响分割系数的准确性,未来采集野外数据并添加至模型可以进一步提高分析结果的可靠性。其次,由于数据缺乏和地下径流过程的复杂性,本研究没有将地下河出口、泉水出露点这部分栅格划入排泄区,仅以平原区漓江河流200 m缓冲区作为排泄区。未来还可以深入挖掘地表、地下产水量的空间变化规律,不断改进模型使其更加符合实际。最后,漓江流域降水的季节性和年际变动较大,可能导致Z系数在年际间的有较大的差异。未来需要增加雨季与旱季的产水量分析,进一步提高模拟结果准确性。

6 结论

本文综合考虑影响喀斯特地区岩溶地貌形态、地表与地下双层径流结构,利用喀斯特地貌类型特征和降雨入渗系数将InVEST模型产水量分割为地表产水量和地下产水量两部分,结合地下水补径排特征和水量平衡法建立产水量校正概念模型,校正了InVEST模型产水量模拟结果,从栅格、岩溶分区和子流域3个尺度分析了漓江流域2000—2020年产水量校正前与校正后时空分布格局。研究结果显示:(1)2000—2020年研究区基于InVEST模型模拟产水量呈现先减后增再减的趋势,产水量空间分布格局为北高南低;就地貌单元而言,岩溶区和非岩溶区产水深度相差不大,不同岩溶区产水能力由高到低排序为裸露型>覆盖型>埋藏型,不同地貌形态产水能力由高到低排序为平原>洼地>峰丛;就子流域而言,流域产水能力降序为漓江上游>漓江中游>恭城河区>荔浦河区>漓江下游,流域分区供水总量由高到低排序为漓江上游>漓江中游>恭城河区>荔浦河区>漓江下游。(2)校正后的水供给服务空间分布格局为西北部高于东南部,四周高于中部;就地貌单元而言,校正后岩溶区产水能力显著高于非岩溶区,不同岩溶区产水能力由高到低排序为埋藏型>覆盖型>裸露型,不同地貌形态产水能力由高到低排序为平原>峰丛>洼地;就子流域而言,流域产水能力降序为漓江中游>漓江上游>恭城河区>荔浦河区>漓江下游。(3)相比InVEST模型模拟结果,基于分割系数校正后的漓江流域产水总量没有变化,但除河流以外的其他区域产水深度均有不同程度下降,不同子流域区域中漓江上游区、荔浦河区、恭城河区的子流域均有减少。利用水文站数据验证发现,校正后产水量模拟值与实测径流值之间的偏差更小,特别是更真实地反映了岩溶地貌区产水量时空格局。因此,本研究更为准确的模拟了漓江流域水供给服务时空格局,研究结果可以为漓江流域水资源管理提供科学依据。未来需要通过补充观测、区分旱季与雨季以及增加地下水文数据来进一步提高模拟结果准确性。