王辉源, 宋进喜,3, 吴 琼
(1.西北大学 城市与环境学院, 西安 710127; 2.陕西省地表系统与资源环境承载力重点实验室, 西安 710127; 3.西北大学 秦岭研究院, 西安 710127)
生态系统服务是人类赖以生存和发展至关重要的资源与环境基础[1-2],水源涵养功能在生态系统服务中占核心地位[3],具有清洁供水、水量调节、养分循环及系统生产力等多项功能[4],良好且稳定的水源涵养功能是区域生态安全、水安全及经济社会可持续发展的基础性保障。然而,由于人类不合理的经济行为,植被破坏、植被种类改变、水生态空间不足、山水林田湖草沙格局不匹配等导致水源涵养功能下降。此外,全球气候变暖及极端气候强度显著增加叠加生态水文效应,改变了水循环速率,导致水源涵养功能不确定进一步增加[5]。改善区域生态系统水源涵养功能,是当前森林生态水文学研究的热点,研究成果具有理论和实践双重价值。
定量化表征水源涵养量是水源涵养功能研究的基础性工作。目前,定量计算水源涵养量的方法主要包括试验法、单统计法和水文模型模拟法。在这些方法中,水文模型模拟法量化评估水源涵养量得到广泛应用。例如Li等[6]采用InVEST产水量模型定量分析了丹江流域水源涵养时空动态,研究水源涵养量对气候、土地利用和土壤变化的响应。田菲[7]基于SWAT模型定量表征了祁连山生态工程建设对水源涵养量的影响,研究了气候变化背景下未来水源涵养量增贮潜力。聂忆黄[8]基于地表能量平衡和SCS模型提出了一种定量评估水源涵养量的方法,划定了祁连山水源涵养能力空间分布范围。然而,目前针对水源涵养功能研究中,对水源涵养量和水源涵养功能内涵认识还存在一定不足,许多研究将水源涵养量与水源涵养功能划为等号,也有研究者采用水源涵养指数作为水源涵养功能评价指标,水源涵养指数只关注降水量与水源涵养量之间关系,忽略了水源涵养量与生态系统本身特征之间相关关系[9]。本文基于生态系统水碳循环作为水源涵养功能评价理论基础,采用水源涵养量与NPP(植被净初级生产力)的比值作为水源涵养功能定量化评价指标,综合考虑气象因素与生态系统结构对水源涵养功能的影响,为不同流域水源涵养功能定量化比较做出初步尝试。
SWAT(soil land water assessment tool)模型一般用于分析计算流域内复杂水文过程长期变化,物理机制较强,因此,本文采用SWAT模型计算水源涵养量。针对SWAT模型对不同季节降雨特征月均流量模拟差异较大问题,本文区别于传统全序列率定,将灞河流域月均流量分为春季、夏季、秋季和冬季4个季节序列分别进行率定验证,然后合并分季节模拟的月均流量率定验证结果,提高模型模拟精度。
在分季节月均流量率定验证结果基础上,首先,定量计算1959—2017年灞河流域水源涵养量;其次,分析近59 a来水源涵养量演变规律及突变特征;最后,构建水源涵养功能计算新方法,建立近18 a来灞河流域水源涵养功能演变特征。研究结果对深化水源涵养功能认识及量化水源涵养功能提供有益的尝试,为流域水资源管理提供数据基础。
灞河流域是西安市重要的生态补偿区,连接秦岭山区和渭河,流域水资源供需矛盾突出[10]。流域面积2 459.31 km2,平均坡度17.09°,地理范围介于108.97°—109.78°E,33.89°—34.43°N,高程354~2 433 m。气候为暖温带大陆性季风气候,年平均气温13.22℃,年平均降水量717.60 mm。流域上游植被茂密,是秦岭生态保护区,土壤以棕壤性土为主;流域中游土地肥沃,农业发达,土壤以褐土性土为主;流域下游,以城市建设用地为主,人为活动剧烈。
文中栅格数据空间分辨率统一为30 m,投影设置为Krasovsky_1940_Albers,坐标GCS_Krasovsky_1940。数据来源如下:(1) 数字高程模型(DEM)来源于ASTER数字地形产品;(2) 土地利用数据采用中国科学院数据中心,1990年土地利用遥感监测数据;(3) 土壤类型数据来源于地理科学生态网;(4) 马渡王站点月均流量数据来源于陕西省水文局;(5) 蓝田气象数据来源于国家气象数据中心,天气发生器数据采用CFSR气象数据模拟[11]。(6) 植被净初级生产力数据(NPP)可由MOD17A3产品获得,由栅格尺度获取小流域尺度NPP值。
SWAT(soil and water assessment tool)是美国农业部(USDA)的农业研究中心Jeff Arnold博士1994年为确定土地管理、植被变化、地下水抽取和水库管理对水质、水量的影响开发的[12-13]。因该模型具有较强的物理基础,目前,在水源涵养功能和水源涵养量研究方面,模型在晋江、祁连山、渭河等地具有成功的应用。SWAT-CUP是用于校准SWAT模型的开源程序,该程序将5种算法链接到SWAT模型,其中SUFI-2算法在大型模型校正中发现非常有效[14]。
2.2.1 水源涵养量的计算方法 在定量表征水源涵养量计算方法中,除水量平衡计算方法外,其他计算方法都存在一定不足[9]。水量平衡法仅考虑流域生态系统流入和流出,认为降水量与蒸散发量以及其他消耗水量的差值即为水源涵养量[15]。SWAT以HRU(水文响应单元)为最小模拟单元,同一个HRU具有相同的土地利用、土壤类型及地形特征。SWAT模型水源涵养量计算方法为:
W=P-E-R
(1)
式中:W为单位时间水源涵养量(mm);P,E分别为单位时间平均降水量(mm)、单位时间实际蒸散发量(mm)和单位时间径流深(mm)。
2.2.2 水源涵养量时间演变规律 本文采用Mann-Kendall非参数统计法来确定年水源涵养量突变点。采用Morlet连续复小波变化(Cmor)分析水源涵养多时间尺度变化特征,通过小波变换方差确定水源涵养变化主周期及周期特征。
由DEM生成河网,经反复调试,河流最小汇水面积500 hm2,土地利用、土壤类型及坡度阈值设定为13%,20%,20%时,最终生成的275个子流域和1 285个HRU(水文响应单元)能够较为准确地刻画研究区范围。马渡王水文站1959—2017年月均流量数据作为率定验证数据,1959—1961年月均流量数据为模型预热期,1962—1990年为率定期,1991—2017年为验证期。
选取决定性系数(R2)和效率系数(NSE)作为模型率定验证评价指标。NSE主要判断水文模型模拟结果拟合程度。
R2计算方式为:
(2)
NSE计算方式为:
(3)
本文采用SWAT-CUP对SWAT模型月均流量模拟进行参数率定,选取14个与流量相关的参数参与参数率定验证,迭代模拟次数300次,具体参数见表1。
4.1.1 基于SWAT模型全序列月均流量率定验证 全序列模型模拟结果见表2,1962—1990年率定期月均流量拟合结果(R2=0.77,NSE=0.76),1991—2017年验证期月均流量拟合结果(R2=0.73,NSE=0.73),率定期的月均流量拟合效果略优于验证期的月均流量拟合效果。Ramanarayanan等[16]建议模型率定验证效果如果R2>0.6,SEN>0.5,认为结果是可以接受的或令人满意。为了辨别不同降水特征下模型模拟精度,把全序列率定期和验证期的率定结果,分别按照春、夏、秋和冬季4个季节分别进行提取,与相对应的季节月均流量观测值进行比较。由表2可知,不同季节模型模拟准确度差别比较大,秋季模型模拟结果较好,率定期与验证期模拟结果在0.82以上,优于其他季节,尤其是秋季率定期模拟结果在0.89以上,模拟准确度较为满意。其他季节模拟结果比较可信,但准确度一般,春季和夏季模拟准确度较为相似,冬季模拟准确度较差,尤其,冬季验证期NSE=-0.69,一般认为NSE<0的时候,模拟结果不可信。
表2 基于SWAT模型全序列月均流量模拟结果Table 2 Simulation results of monthly average flow of the entire series based on SWAT model
4.1.2 基于SWAT模型分季节序列月均流量率定验证 影响月均流量形成的参数季节差异较大,但是,在模型率定时,率定参数统一幅度变化,弱化了不同季节水文过程的差异,导致月均流量模拟效果随季节变化而起伏变化,一般情况下,降水量越少时间段,流量模拟结果越差。为了尽可能提高模型模拟准确度,模型预热期仍设置为1959—1961年,将1962—1990年和1991—2017年月均流量分为春季、夏季、秋季和冬季月均流量,将1962—1990年和1991—2017年分季节月均流量分别进行率定与验证,然后合并分季节率定与验证结果,组成完整月均流量。
由表3和图1可知,分季节模型模拟结果优于全序列模型模拟结果,分季节月均流量率定期与验证期准确度为R2=0.85,NSE=0.84;R2=0.86,NSE=0.85,高于全序列月均流量率定期与验证期准确度(分季节模拟结果约0.85,全序列模拟结果约0.76)。但是,不同季节模拟结果差异仍然比较大,秋季模拟结果仍然高于其他季节,相较于全序列模拟结果,分季节模拟结果在不同季节表现不同,分季节春节率定期模拟结果略微低于全序列春季率定期,其他时间段,分季节率定结果相较于全序列率定结果都有不同程度的提升,尤其,分季节冬季验证期流量模拟结果相较于全序列冬季模拟结果提升幅度较大。分季节与全序列夏秋两季模拟结果都高于春季和冬季,夏秋两季降雨丰沛,多强降雨,流域会出现较突出的超渗产流,春冬两季降水较少,降水强度减弱,产流以蓄满产流或基流为主。一定程度上能够说明,模型更容易体现降水较多时候的流量。
图1 分季节与全序列月均流量模拟值Fig. 1 Simulated monthly average flow values by season and full series
表3 基于SWAT模型分季节月均流量模拟结果Table 3 Based on the SWAT model, the average monthly flow rate simulation results by season are carried out
由图2可知,1959—2017年灞河流域年水源涵养量呈下降趋势,变化率为k=-1.4 mm/a。近59 a来年均水源涵养量为179.46 mm,1983年和1988年年水源涵养量呈现峰值,分别达到414.43,402.47 mm,1995年与1997年年水源涵养量出现谷底,分别为-6.54,15.80 mm。水源涵养量时间演变规律与降水量息息相关,1959—2017年水源涵养量与降水量相关系数为0.87,表明降水量是水源涵养量非常重要的因素。
图2 1959-2017年灞河流域年水源涵养量和降水量Fig. 2 Annual water conservation and precipitation in the Bahe Basin from 1959 to 2017
由图3可知,基于Mann-Kendall非参数统计法,1959—1968年和1982—1992年水源涵养量呈上升趋势,但上升趋势并不十分明显,其余时间段水源涵养量呈下降趋势。1992—2004年UF与UB在0.05置信区间有多个交点,为了判断水源涵养突变点,结合水源涵养量距平数据,1988年为水源涵养量突变点。
图3 灞河流域年水源涵养量M-K突变点检验Fig. 3 The M-K mutation point test of the annual water conservation in the Bahe River Basin
基于1959—2017年灞河流域年水源涵养量距平数据,采用Morlet小波进行小波分析,获取小波方差可以明确水源涵养变化周期。信号强弱用小波系数表示,虚线围住的蓝色部分表示水源涵养量偏小,蓝色颜色越深,水源涵养量越小;实线围住的红色部分表示水源涵养量偏多,红色颜色越深,水源涵养量越大;等直线为0对应该时间尺度突变点。小波方差值表示多时间尺度下时间序列变化周期,波峰表示对应时间尺度周期,波峰值越大周期变化越明显。由图4—6可知,在35 a主周期下存在非常明显的23 a左右周期性水源涵养增加减少趋势。
图4 小波方差Fig. 4 Wavelet variance diagram
图5 35 a主周期Fig.5 35-year main cycle diagram
图6 小波系数实部等值线Fig. 6 Contour plot of real part of wavelet coefficients
植被净初级生产力(net primary production,简写NPP)指植物单位时间和单位面积上生物量(或碳)的净增益[17]。NPP量化了大气中CO2转化为植物生物质的过程,反映了气候变化和人类活动对生态系统陆地植被综合作用结果,是生态系统功能状况的重要指标[18-19]。降水的变化必定引起生态系统功能的变化,水循环与碳循环是生态系统两个重要的生物物理过程。基于此,本文提出了水源涵养功能评价新方法,即:水源涵养量与NPP的比值作为水源涵养功能评价基础。该方法不仅考虑了传统水源涵养功能评价只考虑降水量对生态系统水源涵养功能的影响,同时也考虑了生态系统自身结构、组分和景观特征与水源涵养之间的关系。
由图7—8所示,利用最小汇水面积提取小流域NPP值,获取灞河流域NPP空间分布特征。在区域水热调节综合作用下,灞河流域植被NPP空间分布具有较强的规律性,以2017年为例,灞河流域NPP是485.79 gC/m2,NPP值流域空间分布下游<灞河<灞河中游,灞河下游水域宽阔、地形平坦、经济发达、人为干扰严重;灞河中游地貌以丘陵和台地为主,土地肥沃,水热组合较好,土地利用以耕地和林地为主,植被覆盖较高;灞河上游是秦岭山地保护区,土层较薄、质地较差,林地和草地覆盖度较高。如图9所示,2017年灞河小流域水源涵养量随着坡度和高程增加而增加,2017年灞河流域水源涵养量为165.02 mm,灞河与渭河交汇处水面宽广,蒸发量巨大,水体具有极差的水源涵养功能,水源涵养量为-575 mm,库峪上游植被茂密,降水量丰富,水源涵养功能较强,水源涵养量为352 mm。
图7 2017年灞河流域栅格尺度NPPFig. 7 Grid scale NPP of Bahe River Basin in 2017
图8 2017年灞河流域小流域NPPFig. 8 Bahe Basin NPP in 2017
图9 2017年灞河流域水源涵养量Fig. 9 Water conservation in the Bahe River Basin in 2017
如图10所示,小流域平均水源涵养功能0.33 mm/g,最小是-1.36 mm/g,最大是0.89 mm/g。随着坡度上升,水源涵养功能增强;库峪上游、汤峪、辋川峪上游和倒沟峪水源涵养功能较强,浐河干流、灞河干流和荆沟峪水源涵养功能次之,灞河汇入渭河区域,水源涵养功能最弱。
图10 2017年灞河流域水源涵养功能Fig. 10 Water conservation function of Bahe River Basin in 2017
基于2000—2017年NPP和水源涵养量数据,定量分析2000—2017年水源涵养功能,如图11—12所示,结果表明近18 a来,灞河流域水源涵养功能年际变化总体上呈减少趋势,年际变化率为-0.013 5 mm/g,波动幅度0.678 mm/g,年均水源涵养功能0.291 mm/g,2013年水源涵养功能最低,2003年水源涵养功能最高。近18 a来灞河水源涵养量呈减少趋势,但是,NPP值呈增加趋势,同时,随着温室气体排放增加,气候变暖,植被NPP值会进一步增加,水源涵养功能存在进一步减弱的趋势。
图11 灞河流域2000-2017年水源涵养功能Fig. 11 The water conservation function of the Bahe River Basin from 2000 to 2017
图12 灞河流域2000-2017年水源涵养量及NPPFig. 12 Water conservation and NPP in the Bahe River Basin from 2000 to 2017
(1) 水源涵养是一个具有复杂性和动态性概念[20],认识水源涵养内涵,量化水源涵养过程和能力,有助于提高水源涵养计算结果合理性。水源涵养量是生态系统对水存贮量及形成过程的定量化研究,水源涵养量主要与气象因子相关。水源涵养功能体现的是生态系统自身结构、组分和景观对保持水分的贡献,既要体现自然因素对生态系统作用,更要体现人类不合理经济行为,植被破坏、植被种类改变,水生态空间不足等导致水源涵养功能发生改变。
(2) 为了提高模型月均流量模拟准确度,本文采用分季节率定方法,总体上,分季节率定结果优于全序列,可以提高模型月均流量模拟准确度。但是,以四季作为月均流量特征划分标准,存在一定随机性和不确定性。因此,未来研究考虑引入InVEST模型中季节参数z作为降水及流量特征划分依据[21],以此提高流量模拟精度。与此同时,分季节率定方法引入过多参数,SWAT模型存在“异参同效”问题[22],后期应加强模型参数不确定性研究分析。
(3) 植被净初级生产力的计算大体可以分为野外测量法和模型模拟法[23],野外测量法一般认为是真值,但时间序列较短,也不适用于大尺度NPP计算;模型模拟法适合大尺度大范围长时间序列NPP估计,但是人类活动对生态系统植被的影响一直是研究的难点;遥感-过程耦合模型,考虑了人类活动对生态系统植被的影响,但是现有的遥感数据(NDVI)一般时间序列较短且空间分辨率较大。水循环与碳循环是生态系统两个重要的生物物理过程,长序列水源涵养量数据易于获取,本文选取的植被净初级生产力数据(NPP)可由MOD17A3产品获得,空间分辨率约为467.32 m。采用MOD17A3计算植被净初级生产力在陕西省、黄土高原、秦巴山区等有大量应用。2000—2017年灞河流域年NPP均值为466.195 g/m2,NPP呈增加趋势,年际变化率为7.79 g/m2,灞河流域NPP呈增加趋势与黄土高原、陕西省和秦巴山地NPP演变趋势一致,增加速率高于黄土高原和陕西省,接近秦巴山地[24-26]。但现有的NPP数据时间序列较短,为了建立长序列水源涵养功能评价指标,未来需要加强长序列NPP遥感-过程耦合模型方面研究。
(1) 影响流量形成的参数季节差异较大,但是,在模型率定时,率定参数统一幅度变化,弱化了不同季节水文过程的差异,导致流量模拟效果随季节变化而起伏变化。本文在全序列率定基础上,采用分季节率定及验证,然后合并分季节率定及验证模拟结果,组成完整流量。分季节率定期与验证期模拟效果优于全序列,可以大幅提高流量模拟精度。
(2) 近59 a灞河流域年均水源涵养量179.45 mm,年际水源涵养量呈下降趋势,变化率为k=-1.4 mm/a;2013年开始年水源涵养量呈明显下降趋势;近59 a水源涵养量有多个突变时间点,突变为1988年;采用Morlet小波对灞河流域年水源涵养量距平时序数据进行小波分析,在35 a主周期下水源涵养变化周期为23 a左右。
(3) 本文尝试采用水源涵养量与NPP的比值作为水源涵养功能评价基础,综合考虑生态系统水碳循环影响,为不同区域不同尺度水源涵养功能比较做初步的尝试。结果表明近18 a来,灞河流域水源涵养功能年际变化总体上呈减少趋势,年际变化率为-0.013 5 mm/g,波动幅度0.678 mm/g,年均水源涵养功能0.291 mm/g,2013年水源涵养功能最低,2003年水源涵养功能最高。水源涵养功能在空间上的表现为随着坡度上升,水源涵养功能增强。