張阿龍,于慶峰,2,于 嬋,房麗晶,高瑞忠
(1.內蒙古農業大學水利與土木建筑工程學院,呼和浩特010018;2.內蒙古農業大學職業技術學院,內蒙古 包頭014109;3.內蒙古自治區水文總局,呼和浩特010020)
近年來,氣候波動和人類活動導致極端降水事件頻繁、徑流銳減[1-3],流域產匯流過程發生變化[4-6],進而草原流域草地退化、土地沙化和沙塵暴等生態環境問題日益突出。國內外學者對于流域產匯流過程從物理模擬、數學及數值模擬進行研究分析,管曉祥[7]對VIC 模型、新安江模型、WBM 模型和GR4J 模型在黃河典型子流域徑流過程模擬效果進行對比;Hamed[8]和Martinez[9]發現,在相同的降水歷時內,徑流量和產沙量對降水強度存在正向響應;姚沖[10]通過模擬降水對南方典型黏土坡面進行研究,查明降水量、降水強度和坡度對坡面產流產沙量和土壤侵蝕過程的影響,并建立坡面土壤侵蝕模型;SWAT模型在域水文研究中應用廣泛,Aronld[11]、Van Liew[12]運用SWAT 模型分別對美國不同流域徑流和小沃希托河流域進行水文模擬,均取得了滿意的模擬效果,Demirel[13]同時建立人工神經網絡模型和SWAT模型,得出SWAT 模型顯示的均方誤差值更小,模擬效果更好。
綜上所述,人工物理模擬降水-徑流試驗對地表產匯流模擬具有科學性,SWAT 模型對水文模擬效果具有優越性,因此,本文將人工物理模擬降水-徑流試驗和SWAT 模型共同引入徑流模擬及變化歸因分析中,研究草原流域產匯流特征,定量分析驅動徑流變化的氣候波動與人類活動的貢獻,旨在探討環境變化對草原天然徑流演變的影響,進而為草原流域植被建設、生態環境保護和生態水文研究提供參考。
巴拉格爾流域(117°32'~118°30'E,43°26'~44°39'N)位于內蒙古錫林郭勒盟(圖1),是以針茅和羊草為主的典型草原內陸河流域,其流域面積2 703 km2,冬季寒冷、夏季干熱,平均氣溫2.39 ℃、年均降水量330.3 mm;流域主要經濟產業為畜牧業和工業[2,5]。
土壤植被數據通過野外調查獲取。結合流域多年降水資料、植被特征、地形地貌和草地管理方式等設計試驗方案,共布設取樣點15個(圖1),分別測得0~30 cm 土壤干容重、飽和滲透系數、飽和含水率、毛管上升含水率和田間持水率;選取50 cm×50 cm 的植被樣方,將鮮草剪去裝入保鮮袋,并在該區挖取深30 cm 土層,篩洗出根系,裝入保鮮袋,在實驗室進行烘干稱重測取地表和地下植被干物質量。

圖1 研究區位置及土壤植被采樣點分布Fig.1 Location of balaguer river basin and distribution of sampling points
降水-徑流模型由供水、降水和徑流裝置3 部分構成(圖2),供水裝置由水箱,水泵,供水槽組成,由雨強調節器控制水量大小,人工模擬雨強均勻的降水,降水器支架為角鋼焊接成的矩形框架,降水量大小利用雨量筒測定,經過多次重復試驗后在雨強調節器上標定。降水-徑流模擬試驗方案設計為在典型區隨機采集100(長)cm×50(寬)cm×30(深)cm 的低擾動土樣塊,進行不同降水強度下草原產匯流及水土侵蝕過程的試驗模擬,分別設置不同降水強度、持續時間和土地坡度等多組合模擬場景,即降水強度為60、90 和120 mm/h,降水持續時間為1、2、3、4、5、6、7、8、9、10、11、12、13、14、15 min,坡度為10°、16°、28°等。

圖2 降水-徑流模擬裝置示意圖(單位:cm)Fig.2 Schematic diagram of simulated precipitation-runoff device
研究區氣候數據包括降水、氣溫、風速、輻射和相對濕度等,來源于研究區內氣象站實時監測和中國氣象數據網(http://data.cma.cn),DEM 數據來源于地理空間數據云(http://www.gscloud.cn),土壤數據來源于寒區旱區科學數據中心(http://wastdc.westgis.ac.cn)(圖3),土地利用數據來源于美國地質勘查局(United States Geological Survey,USGS)數據共享網站(http://glovis.usgs.gov/)提供的Landsat TM/OLI 影像(時間分辨率16 d,空間分辨率為30 m)解譯獲取(圖4)。

圖3 土壤類型分布Fig.3 Distribution of the soil type

圖4 土地利用和植被覆蓋分布Fig.4 Distribution of the land use and vegetation cover
統計分析采用Excel@2010軟件完成,土壤-植被-根系的空間分布圖通過ArcGIS@10.3軟件制作。
Mann-Kendall(M-K)檢驗是一種廣泛用于時間序列趨勢和突變分析,該法計算簡便,可以識別突變開始時間及指出突變區域[14]。
克里金插值法是常用的地質統計格網化方法,其考慮了變量空間相關特性,廣泛應用于地質、水文、氣象等領域的空間分析[15]。
利用降水-徑流模擬試驗建立不同坡度降水強度與產流系數的數理回歸模型,由降水強度經驗頻率與降水強度、不同坡度空間分布面積占比的變化范圍分別劃分區間,計算各區間經驗頻率,公式如下:

式中:Ri為坡度在i區間的徑流系數;Qi為坡度在i區間的面積占比;Pj為降水強度在j區間的經驗頻率;rij為坡度在i區間,經驗頻率在j區間的徑流系數;i為坡度,0°~10°、10°~16°、16°~28°;j為降水強度,0~10、10~60、60~90、90~120 mm/h;R為模擬徑流系數。
SWAT 模型進行水文演變驅動因素分析,以Nash-Suttcliffe效率系數(Ens)和決定系數(R2)進行模型參數率定和模型驗證評價及檢驗模擬擬合過程是否合理[17]。
以流域徑流突變前的多年平均徑流量為天然時期流域徑流的基準值,則天然徑流量與徑流突變后期的多年實測平均徑流量之間的差值就是受氣候變化與人類活動的影響部分,因此,氣候變化與人類活動的影響分割公式如下[18]:
式中:QH為突變后的實測徑流;QN為突變前徑實測流均值;QA為突變前實測徑流均值;QB為突變前模擬徑流均值;QC為突變后模擬徑流均值;QCI為修正突變后模擬徑流均值;I為修正系數;ΔQC為氣候變化引起的均值變化度為;ΔQH為人類活動引起的均值變化度;WC為氣候變化引起的貢獻率;WH為人類活動引起的貢獻率。
通過M-K 檢驗法得到氣象水文要素的突變年份(表1),突變主要集中在20 世紀90年代到20 世紀末,1982年最低氣溫突變最早,2001年降水量突變最晚,1988年徑流深突變,突變前后變化率為-35.97%,而徑流系數突變年份為1994年,突變前后變化率為-38.57%,徑流深與徑流系數突變前后變化率較為接近。

表1 氣象水文要素突變前后變化Tab.1 Distinction of the hydrometeorological factors before and after the shift year
土壤物理參數(表2)中飽和滲透系數介于0.19~14.23 m/d,但均值僅為2.87 m/d,表明土壤物理性質空間異質性顯著,干容重最大值、最小值與均值接近,毛管上升含水率極差最大,飽和含水率次之。

表2 土壤物理特征統計Tab.2 Statistical parameters of the soil physical characteristics
流域下墊面特征與土壤質地、植被分布和植被根系密切相關。研究區土壤類型、草干重、草根干重空間分布特征見圖5,土壤以砂質壤土為主,粉砂質壤土相間的多礫沙質壤土和少礫砂質壤土呈島狀分布的分布特征,草干重最大值80~90 g 主要分布在上游山谷區,分布面積較小,呈現出中部最少,由南向北,由東向西遞減的趨勢,草根重以200~300 g 分布區域為主,最大值介于400~500 g,分布于上游東北部,100~200 g呈星點狀分布,草根重于草種分布存在較大差異。

圖5 土壤類型、草干重、草根干重空間分布圖Fig.5 Spatial Distribution of soil type,grass dry weight and grass root dry weight
研究區多年降水頻率分布曲線表明流域降水強度主要集中在0~10 mm/h,約占80%,地面坡度空間分布圖(圖6)中流域坡度主要集中在0~10°,約占流域面積的86%。

圖6 流域坡度空間分布圖Fig.6 Distribution of the basin slope
根據降水經驗頻率特征和不同坡度空間分布面積占比,分別將降水-徑流模擬試驗分為4組和3組(表3)進行。

表3 流域坡度面積與降水強度經驗頻率占比Tab.3 Proportion of different slope area and empirical frequency of precipitation Intensity
降水徑流模擬試驗中每場降水前土壤平均含水率保持在10%~15%之間,產匯流特征參數見表4,隨著降水強度和坡度的增加,徑流總量和徑流系數均增大,地面坡度為10°時,降水強度為120 mm/h 的產流總量是降水強度為90 mm/h 的2.19 倍,是降水強度為60 mm/h 的4.69 倍;而地面坡度為16°時,降水強度為120 mm/h 的產流總量是90 mm/h 的2.07 倍,分別是60 mm/h 的4.71 倍;而坡度為28°時,降水強度120 mm/h 的產流總量是90 mm/h的1.92倍,分別是60 mm/h的4.66倍,綜上,總產流量不僅與降水強度有關,而且受坡度影響較大。當坡度為10°時,3種降水強度產流機制主要以蓄滿產流為主,隨著坡度的增加,產流量不斷加大,產流機制由蓄滿產流為主逐漸向超滲產流過度,徑流系數加大。

表4 不同坡度和降水強度下產流特征參數Tab.4 Statistical parameters of runoff characteristics
根據模擬降水試驗結果,建立不同坡度降水強度與產流系數的多項式回歸方程(圖7),各回歸方程R2值均大于99%,擬合效果較好。

圖7 不同坡度、降水強度與徑流系數關系Fig.7 Relationship between slope,rainfall intensity and runoff coefficient
不同區間降水強度值以中值代替,不同區間坡度值以最大值代替,計算不同坡度與降水強度對應的徑流系數,結果如表5。

表5 不同坡度和降水強度下的徑流系數Tab.5 Runoff coefficient of the different slope and rainfall intensity
根據公式(1)、(2)計算得到天然徑流系數R值為0.424(表6),可看出該區坡度為16°~28°區間徑流系數相對較大,而坡度為0°~10°區間徑流系數較小,表明流域坡度對徑流系數影響較大。

表6 流域天然徑流系數計算Tab.6 Calculation of the basin runoff coefficient
基于人工降水模擬試驗建立降水強度-坡度-徑流系數的數學模型結果,可得突變前人類活動對巴拉格爾流域徑流系數變化的貢獻率為24%,氣候變化貢獻率為76%,突變后人類活動對巴拉格爾流域徑流系數變化的貢獻率為53%,而氣候變化貢獻率為47%(表7)。

表7 流域徑流系數變化歸因分析Tab.7 Attribution analysis of runoff coefficient
基于水文氣象要素實測數據,利用SWAT 模型研究流域徑流變化驅動因素,率定期模擬值和實測值的納什系數Ens為0.88,決定系數R2為0.87,驗證期模擬值和實測值的Ens為0.78,R2為0.81,由SWAT 的評價標準可知,SWAT 模型在該流域模擬效果較好,驗證期的模擬值總體略大于實測值,但徑流極大值實測值均大于模擬值(圖8)。

圖8 流域年徑流模擬值與實際值對比Fig.8 Comparison between the observed and simulated runoff
由公式(3)~(9)分析流域氣候和人類活動驅動影響徑流的貢獻,氣候變化貢獻率為79%,人類活動貢獻率為21%,徑流變化主要受氣候變化影響,而受到的人類活動影響較?。ū?)。

表8 流域徑流變化歸因分析Tab.8 Attribution analysis of basin runoff
巴拉格爾流域水文氣象要素突變主要發生在20 世紀90年代到20 世紀末,突變前后水文氣象要素差異較大,主要表現為氣溫升高,降水量、徑流量減少,生態環境惡化加劇,導致生態環境與氣象水文要素互饋陷入惡性循環,人類活動(對地下水吸奪,過度放牧,工業興起,城市擴張等)加重了惡化趨勢,造成巴拉格爾河豐水期斷流、草場退化等生態水文現象凸顯。Zhang等[2]通過對巴拉格爾流域氣象水文要素突變時序研究,發現平均最低氣溫、平均氣溫、平均最高氣溫、徑流量、降水對環境敏感性由強到弱,由氣溫的突然升高導致了巴拉格爾河流域春季徑流量的急劇下降,因春季徑流占全年的40%以上,春季徑流的突變減少和溫度的突變增加,導致降水的突然減少,突變時序與白勇[19]、高瑞忠等[20]對巴拉格爾流域水文氣象突變年份一致;土壤類型和植被的空間分布對產匯流過程影響較大,不同的土壤類型滲透性差異較大,不同的植被對降水的截留和產匯流速率的影響較大,由于研究區土壤類型的質地和土壤水力參數相近,植被分布主要以草地為主,其他地物類型星點分布,雖然對徑流系數有一定的影響,但相對流域整體而言可忽略不計。魯克新等[21]對黃土區降水產流產沙過程研究認為裸坡坡面下,認為雨強是影響徑流的主要因素,坡度是影響產沙的主要因素,忽略了模擬降水過程中蒸散發作用的影響和流域下墊面的空間異質性等因素的干擾,徑流系數可以間接反映流域內自然地理要素對降水-徑流關系的影響,由于人工模擬降水-徑流關系可排除人類活動造成的干擾,且各要素均可人為控制,根據研究區地形地貌多年水文氣象資料計算出各要素的分布頻率,再通過人工降水模擬得理想條件較接近天然狀態的徑流系數,雖然該方法對天然徑流系數的刻畫存在一定的誤差,但從模擬結果看,通過人工降水模擬試驗建立降水強度-坡度-徑流系數的數學模型得出天然徑流系數值為0.424,模擬值與突變前的徑流系數接近,較突變前多年平均實測徑流系數大,因巴拉格爾河流域產流機制主要為超滲產流,多年平均徑流系數的改變主要受下墊面條件控制,下墊面的變化主要受人類活動的影響;試驗結果表明突變前人類活動對巴拉格爾流域徑流系數變化的貢獻率為20%左右,突變后為50%左右。利用分布式水文模型SWAT對徑流變化歸因分析發現巴拉格爾流域上游氣候變化貢獻率為79%,人類活動貢獻率為21%,徑流變化主要受氣候變化影響,而受到的人類活動影響較小,與人工降水模擬試驗建立降水強度-坡度-徑流系數的數學模型得出的貢獻率存在較大差異,由于SWAT 模型假設突變前的多年平均徑流量為天然時期流域徑流的基準值,縮小了人類活動貢獻率,進而導致突變后人類活動對徑流的影響的貢獻率偏小,但SWAT 模型模擬結果與蘇輝東等[22]對黃河、黑河、長江、雅魯藏布江流域四個寒區流域得出氣候變化對徑流的影響貢獻率達到78%以上是徑流演變的主導作用的結論相近,而與高瑞忠[20]對巴拉格爾流域和王威娜[23]、于嬋[24]等對錫林河流域的研究相反,可能是使用的水文數據序列和徑流歸因分析的方法不同造成。
變化環境下徑流量的減少逐年加劇,由于突變后經濟飛速發展,人類活動范圍的不斷擴大(過度放牧對草場的破壞、地下水開采對徑流量的吸奪、礦產資源開發對地下水的輸干和對下墊面條件的破壞等)和全球持續增溫等因素的影響,導致草原生態系統退化,地下水位遞減,草場面積萎縮等環境問題突出。氣溫突變增高導致徑流減少,徑流突變減少和氣溫的突變增高共同導致了降水的突變減少,促使巴拉格爾流域水文氣象極端事件增多,干旱加劇,水文氣象要素與生態環境互饋機制改變,但氣候變化也受人類活動的干擾。張樹磊等[25]認為在氣候較為干燥的地區,徑流對氣候和下墊面變化都更為敏感,且區域差異性明顯;劉昌明等[26]和張曉明等[27]認為森林植被增加不僅能顯著減少豐水、平水及枯水期年的徑流,也改變了地下水與地表水的分配。
學界對高原內陸河草原流域降水強度-坡度-徑流系數之間關系的深入研究相對較少,需進一步探討。對不同降水模式下流域產匯流試驗研究中由于缺乏不同植被蓋度的對比研究,降水-徑流模擬設備還需進一步改進完善,需要加入測量土壤實時含水率、徑流無間斷測定等設備,考慮壤中流與地下徑流的水文過程,以及各徑流中的化學成分變化等,徑流系數模型建立未考慮蒸散發,壤中流和地下徑流等,由于降水強度取值較少,對降水強度與徑流系數的擬合曲線可能存在一定的影響,對徑流系數數學模型還需進一步研究,對徑流歸因分析方法間的差異仍需要進一步探究。
(1)流域水文氣象要素突變主要集中在20 世紀90年代到20世紀末,最低氣溫突變最早,降水量突變最晚。
(2)流域產流量不僅與降水強度有關,而且受坡度影響較大,當坡度為10°時,3 種降水強度產流機制主要以蓄滿產流為主,隨著坡度的增加,產流量不斷加大,產流機制由蓄滿產流為主逐漸向超滲產流過度,徑流系數加大。
(3)根據人工降水模擬試驗發現天然徑流系數R值為0.424,突變前人類活動對巴拉格爾流域徑流系數變化的貢獻率為20%左右,突變后為50%左右。
(4)基于SWAT 模型研究流域徑流變化驅動因素表明巴拉格爾流域突變后徑流變化的人類活動貢獻率為21%。