鄭明國
(1. 水利部黃土高原水土流失過程與控制重點實驗室,鄭州 450003;2. 廣東省生態(tài)環(huán)境技術(shù)研究所廣東省農(nóng)業(yè)環(huán)境綜合治理重點實驗室/廣東省面源污染防治技術(shù)工程中心/華南土壤污染控制與修復(fù)國家地方聯(lián)合工程研究中心,廣州 510650;3. 中國科學(xué)院地理科學(xué)與資源研究所陸地水循環(huán)及地表過程實驗室,北京 100101)
河流水沙變化會帶來諸多的環(huán)境問題。從全球角度言,河流水沙及其挾帶的養(yǎng)分物質(zhì)的變化會導(dǎo)致生物地球化學(xué)循環(huán)過程的相應(yīng)變化;從區(qū)域角度言,河流水沙變化不僅關(guān)系到河流本身的健康,而且在水利規(guī)劃、防洪減災(zāi)、水資源利用保護(hù)、水土保持及生態(tài)環(huán)境建設(shè)等方面也是須考慮的重要因素。作為全球變化研究的一個部分,氣候變化和人類活動影響下的流域水沙過程研究已成為一個全世界廣泛關(guān)注的重要科學(xué)問題[1]。
作為典型的多沙河流,黃河的泥沙問題為中國歷代政府所重視。自20 世紀(jì)80 年代中期以來,黃河水沙情勢發(fā)生了巨大變化[2],引起了社會各界的廣泛關(guān)注。黃河水沙變化研究已成為中國水科學(xué)領(lǐng)域的重大科學(xué)問題之一,有眾多的科研項目對此問題展開了研究(詳見文獻(xiàn)[3])。這些項目研究中,基本都設(shè)置有水沙變化貢獻(xiàn)分割相關(guān)研究內(nèi)容,即定量評估氣候因素和各種人類活動對觀測的水沙變化的貢獻(xiàn),藉此確定黃河水沙變化的主導(dǎo)因素,從而為治黃決策服務(wù)。然而,這些項目的貢獻(xiàn)分割結(jié)果并不一致,如冉大川等[4]和劉萬銓[5]比較了2000 年以前主要的水土保持減水減沙效益計算成果,發(fā)現(xiàn)各成果之間差異較大,難以作為宏觀決策的依據(jù)[5]。
針對黃河水沙變化貢獻(xiàn)分割研究,已提出有大量方法,但仍以水文法和水保法為主[3]。水文法首先建立基準(zhǔn)期(即非治理期)的降雨和產(chǎn)沙關(guān)系,然后將治理期的降雨因子代入,得到僅降雨發(fā)生變化時的產(chǎn)沙量,其與基準(zhǔn)期產(chǎn)沙的差值即為降雨對產(chǎn)沙變化的貢獻(xiàn),降雨以外的其他貢獻(xiàn)全部歸為人類活動結(jié)果。該方法的計算精度主要取決于所建立的年降雨產(chǎn)沙回歸關(guān)系的精度。黃河中游年降雨產(chǎn)沙關(guān)系通常表現(xiàn)出一定程度的離散,利用該關(guān)系式計算年產(chǎn)沙會有一定的誤差,但回歸方程計算結(jié)果為隨機(jī)變量的條件均值,因此利用該關(guān)系式計算多年均值時誤差會小很多。這樣,如忽略基準(zhǔn)期和治理期兩個時段的雨強(qiáng)變化,水文法是一種可靠的產(chǎn)沙變化貢獻(xiàn)分割方法[6],劉萬銓[5]甚至建議以水文法為基準(zhǔn)對水保法結(jié)果進(jìn)行檢驗和校正。但水文法只能區(qū)分人類活動和氣候變化的貢獻(xiàn),無法進(jìn)一步區(qū)分不同人類活動,如不同水土保持治理措施的貢獻(xiàn)。
與水文法不同,水保法可以區(qū)分不同水土保持治理措施對流域產(chǎn)沙變化的貢獻(xiàn)。該方法首先根據(jù)坡面小區(qū)觀測結(jié)果,確定各水土保持措施的減蝕率,然后統(tǒng)計各措施治理面積,在此基礎(chǔ)上得到各項水土保持措施的減沙貢獻(xiàn)。水保法存在諸如小區(qū)代表性(包括空間、時間以及措施強(qiáng)度代表性等)等問題[6],但更為嚴(yán)重的是,水保法完全忽略了減沙效益的空間尺度變異,理論基礎(chǔ)存有嚴(yán)重缺陷。世界上許多地區(qū)的研究表明[7-8],盡管坡面治理后坡面侵蝕強(qiáng)度變小,但水流進(jìn)入溝道后會繼續(xù)沖沙,可能會導(dǎo)致溝道侵蝕增強(qiáng),使得流域出口產(chǎn)沙量變化小于坡面甚至保持不變。如對巴西一個1.19km2小流域的研究表明,水土保持措施實施后坡面對流域出口的產(chǎn)沙貢獻(xiàn)降低,但同時溝道侵蝕的貢獻(xiàn)卻變大,導(dǎo)致流域出口的洪水平均含沙量治理前后保持不變[8],同樣的現(xiàn)象在黃土高原也有報道[9-10]。以坡面減蝕率直接代替流域出口的減沙率,顯然會導(dǎo)致對水土保持措施減沙貢獻(xiàn)的高估。
水沙變化貢獻(xiàn)分割的精度取決于數(shù)據(jù)的質(zhì)量以及方法的合理性兩個方面[5,11]。目前,隨著遙感等新技術(shù)的進(jìn)步以及投入經(jīng)費的增大,數(shù)據(jù)質(zhì)量問題已有了很大改善,但仍缺乏有效的貢獻(xiàn)分割方法[6],使得難以準(zhǔn)確評估黃土高原不同治理措施的減沙效益[6,12-13]。黃土高原水土保持治理已取得了巨大成效,但減沙貢獻(xiàn)分割方法的滯后使得治理效果的歸因目前存有較大爭議,這嚴(yán)重影響了治黃策略的制定。
黃土區(qū)的水保治理措施可分為坡面(主要包括植被和梯田措施)和溝道措施(主要指淤地壩)兩大類?;趯Ρ仍囼灹饔蚬灿?5a 的觀測資料,鄭明國等[10]認(rèn)為,在流域尺度上,植被、梯田等坡面措施僅通過減少產(chǎn)流來減沙,并不改變流域出口的水沙關(guān)系?;诖?,本文提出了一種新的黃土區(qū)流域減沙貢獻(xiàn)分割方法。首先分析了在水沙關(guān)系框架下,不同水保措施減沙機(jī)制的異同;在此基礎(chǔ)上,給出了不同水保措施及降雨的減沙計算公式;最后,通過與水文法計算結(jié)果以及實測數(shù)據(jù)的比較,對該方法進(jìn)行了驗證。
黃河的泥沙主要來自中游的黃土高原。黃土高原為典型溫帶半干旱氣候,降雨量集中在400~600 mm 之間,大部分地區(qū)黃土厚度超過100 m,土質(zhì)疏松,抗蝕性差,在大規(guī)模水土保持治理前,許多地區(qū)的土壤侵蝕強(qiáng)度超過10 000 t/km2·a。長期的高強(qiáng)度侵蝕使得本地區(qū)溝壑密布,地形起伏大,重力侵蝕如崩塌等頻繁發(fā)生[14]。
陜北是黃土高原淤地壩建設(shè)最為集中的區(qū)域。本文選擇陜北3 條典型河流—延河、清澗河、大理河為研究區(qū)(圖1),用于驗證淤地壩減沙效益計算結(jié)果。1970年以前,當(dāng)?shù)剞r(nóng)業(yè)活動強(qiáng)烈,水土保持治理措施較少。1970 年以后開展大規(guī)模水土保持治理,治理措施主要包括植樹、種草、修建梯田和淤地壩等,2000 年以前流域減沙以淤地壩為主,各流域歷年治理面積見表1。根據(jù)所搜集到的淤地壩攔沙數(shù)據(jù),具體驗證工作在甘谷驛、延川、綏德、李家河、曹坪5 個水文站點展開(圖1 和表2)。
為驗證坡面措施減沙效益計算結(jié)果,本文也選擇山西離石站和甘肅天水站的兩對水土保持對比溝為研究區(qū)。山西離石站位于黃土丘陵溝壑區(qū)第一副區(qū),1956—1970 年期間曾對離石縣王家溝流域內(nèi)兩個毗鄰小流域—羊道溝和插財主溝進(jìn)行過對比觀測。羊道溝流域面積0.21 km2,完全未經(jīng)治理,植被稀疏,農(nóng)地面積占總面積的59.9%。插財主溝流域面積0.19 km2,自1956 年開始治理,治理措施有梯田、地梗、草田輪作、造林以及封禁等,以林草為主,無淤地壩措施,觀測后期治理面積為78.3%。由于嚴(yán)禁人為破壞,植被結(jié)構(gòu)已和天然林區(qū)近似,林下有枯枝落葉層和地衣苔蘚。羊道溝觀測期間年均產(chǎn)沙模數(shù)為20 811 t/km2,插財主溝為8 503 t/km2,流域治理的減沙率為59.1%。
圖1 大理河、清澗河、延河流域監(jiān)測站點位置 Fig.1 Locations of the monitoring stations in the Dalihe, Qingjianhe, and Yanhe river basins
表1 延河、清澗河及大理河流域不同水土保持措施治理面積 Table 1 Cumulative areas of various soil and water conservation measures in the basins of Yanhe, Qingjianhe and Dalihe River km2
甘肅天水站位于黃土丘陵溝壑區(qū)第三副區(qū),2 個對比觀測流域—橋子西溝和橋子?xùn)|溝,空間上相鄰,自然特性也相似。橋子西溝流域面積1.09 km2,為非治理流域。橋子?xùn)|溝為治理流域,流域面積1.36 km2,于1987年開始治理,2006 年以前以造林為主,隨后開始修建有淤地壩[16]。在2006 年,橋子西溝坡耕地面積占比36.4%,林地為49.3%,橋子?xùn)|溝坡耕地僅占4.5%,林地面積比例高達(dá)88.8%[17]。1987—2010 年期間,橋子西溝年均產(chǎn)沙模數(shù)為5 756 t/km2,橋子?xùn)|溝為2 005 t/km2[17],流域治理的減沙率65.2%。
本文研究所使用的降雨、徑流和泥沙數(shù)據(jù),橋子?xùn)|溝和橋子西溝為天水站觀測資料;羊道溝和插財主溝為山西省水土保持科學(xué)研究所觀測,摘錄自《山西省水土保持科學(xué)研究所徑流測驗資料(1955—1981)》;曹坪站1970 年以前摘錄自《黃河流域子洲徑流試驗站水文實驗資料(1959—1969)》,其他來自《中華人民共和國水文年鑒》。橋子?xùn)|溝和橋子西溝數(shù)據(jù)年限為1987—2006年,羊道溝和插財主溝為1956—1970 年;延河、清澗河、大理河數(shù)據(jù)為1950s—1989 年期間觀測,各站點觀測年限均超過30 年。為方便,以下稱1970 年以前時段為P0 時段,1970—1979 年為P1 時段,1980—1989 年為P2 時段,1970—1989 年為P3 時段。本文中“產(chǎn)沙”一詞指觀測時段內(nèi)通過河流某一斷面的懸移質(zhì)總量,產(chǎn)沙模數(shù)指單位流域面積的產(chǎn)沙量。
為保證數(shù)據(jù)精度,本文搜集的淤地壩攔沙量數(shù)據(jù)均為逐壩測量結(jié)果,列于表2 中Δ 列。Δ 中也包括水庫淤積量,但所占比例一般較小,故下文中稱其為淤地壩攔沙量。Δ 計算公式為:Δ=aλV,V 為淤地壩和水庫內(nèi)淤積體體積,λ 為淤泥干容重,取1.35 t/m3。由于Δ 包括了推移質(zhì)及淤地壩修建前的填土墊方量,因此需按一定系數(shù)折減,a 為折減系數(shù)。對淤地壩,a 取0.9;對水庫,a 取0.95。V 數(shù)據(jù)的施測時間為表2 中對應(yīng)時段的終止年份。在1989 年,陜西省水土保持局對陜北淤地壩展開了一次全面普查。延河、清澗河、大理河和李家河的V 數(shù)據(jù)均為該次普查期間測量結(jié)果,數(shù)據(jù)摘錄自惠養(yǎng)瑜等[15]和張勝利等[18]。岔巴溝的V 數(shù)據(jù)摘錄自徐建華等[19],為黃河水利委員會吳堡水文總站和子洲縣水利局測量結(jié)果。
表2 中ΔSYCD為根據(jù)Δ 計算的對應(yīng)時段內(nèi)壩庫年均攔沙量。對延河、清澗河、大理河流域,惠養(yǎng)瑜等[15]根據(jù)歷年淤地壩控制面積和降雨估算了各時段的淤積比例,然后將Δ 分配到P0、P1、P2 各時段,本研究按相同比例進(jìn)行Δ 在各時段的分配;李家河站的分配比例直接采用綏德站的結(jié)果;曹坪站P0 時段的淤積量采用徐建華等[19]的計算結(jié)果。下文將以ΔSYCD為標(biāo)準(zhǔn)對淤地壩攔沙量計算結(jié)果進(jìn)行驗證,并稱之為“實測攔沙量”。
表2 研究區(qū)各時段降雨、徑流、泥沙數(shù)據(jù)及坡面措施減水和淤地壩攔沙量 Table 2 Rainfall, runoff , sediment yield, reduced runoff by slope measures, and trapped sediment by check dams for the selected catchments
表2 中坡面措施的年均減水量,即ΔRSM全部為地表徑流(surface runoff)組分,無基流組分?;蒺B(yǎng)瑜等[15]根據(jù)具有長期觀測記錄徑流小區(qū)資料,確定了不同質(zhì)量等級的梯田、林地和草地的減水效益,并按豐、平、枯水年進(jìn)行了校正,最終確定了延河、清澗河、大理河的梯田、林地和草地不同時期的減水量,對該成果數(shù)據(jù)進(jìn)行匯總后即得到ΔRSM;李家河流域按大理河的數(shù)據(jù)估算;岔巴溝數(shù)據(jù)為徐建華等[19]計算結(jié)果。
下文計算過程中假定P0 時段無水保措施,故利用淤地壩攔沙量(即表2 中ΔSYCD)和水保措施減水量(表2中ΔRSM)對P0 時段的SSY、Q 數(shù)據(jù)進(jìn)行了還原。不過,由于P0 時段的水保措施很少,還原后SSY 和Q 變化不大,還原與否對最后的計算結(jié)果影響很小。
1.2.1 基流分割
本文中地表徑流數(shù)據(jù)利用總徑流減去基流獲取,基流分割采用廣泛應(yīng)用的Lyne and Hollick 濾波方法[20]。文獻(xiàn)[21]推薦使用該方法時,濾波參數(shù)a 取值為0.98,并連續(xù)進(jìn)行三輪濾波處理(一次前向濾波,一次后向濾波,最后再進(jìn)行一次前向濾波)。本文使用該方法時,a 取值為0.925,且僅對日徑流時間序列數(shù)據(jù)進(jìn)行了一次前向濾波處理。黃河中游8 條支流的數(shù)據(jù)驗證結(jié)果表明,按此設(shè)置的計算結(jié)果和Zheng[22]提出的水沙關(guān)系法計算結(jié)果更為接近。Zheng[22]認(rèn)為,徑流輸沙關(guān)系中的截距項即為研究時段內(nèi)的平均基流量。濾波法把所有的低頻組分全部歸類為基流,因此會系統(tǒng)地高估基流。水沙關(guān)系法的計算結(jié)果比濾波法更為接近真值,但當(dāng)徑流輸沙關(guān)系不好時,如枯水年或少沙年份,該方法難以應(yīng)用。
1.2.2 減水貢獻(xiàn)估算
本文計算過程中需要定量評估降雨以及各水保措施對地表徑流的影響。地表徑流主要產(chǎn)生于汛期,故首先建立各研究區(qū)在P0 時段的年汛期降雨(PFa,mm)和年地表徑流(Ra,mm)的回歸關(guān)系,設(shè)為Ra=a0PFa+b0。分析表明,與將其他月份(如6-9 月或5-9 月)作為汛期比,采用5-8 月為汛期建立的Ra和PFa回歸方程的R2普遍較高,故本文采用5-8 月降雨量為PFa值。設(shè)P0 時段的多年平均汛期雨量為PF0(mm/a),流域治理后某一時段Pi 的多年平均汛期雨量為PFi(mm/a),則與P0 時段比,Pi 時段由于降雨變化導(dǎo)致的地表徑流變化量ΔRP(mm /a)可按下式計算
各坡面措施的減水量直接參考相關(guān)文獻(xiàn)結(jié)果,詳見表 2 中ΔRSM列及1.1 部分說明。在2000 年以前,黃土區(qū)的植被改善和梯田修建面積很有限[11],故也可假定ΔRSM=0,下文也針對此情形進(jìn)行了相關(guān)計算。設(shè)R0和Ri分別為P0和Pi 時段的年均地表徑流(mm/a),ΔR= Ri-R0。排除ΔRSM和ΔRP后,ΔR 中剩余部分可認(rèn)為是Pi 時段內(nèi)為淤地壩的年均減水量ΔRCD(mm/a):
1.2.3 利用水文法計算人類活動減沙貢獻(xiàn)
首先建立P0 時段的年汛期降雨P(guān)Fa和年產(chǎn)沙模數(shù)SSYa(t/km2)的回歸關(guān)系,設(shè)為SSYa=d0PFa+e0。按水文法,降雨變化對產(chǎn)沙變化的貢獻(xiàn)為d0(PFi-PF0)。如降雨外的因素全部歸于人類活動的結(jié)果,則人類活動的減沙貢獻(xiàn)ΔSSYH(t/km2·a)可估計為:
式中ΔSSY = SSYi-SSY0,SSY0和SSYi分別為P0 和Pi 時段的年均產(chǎn)沙模數(shù)(t/km2·a)。
表3 中給出了研究區(qū)各流域P0 時段的PFa-Ra和PFa-SSYa回歸關(guān)系。其中PFa為面數(shù)據(jù),雨量站分布見圖 1。在年尺度,各研究流域內(nèi)不同站點雨量差別不大,因此直接根據(jù)算術(shù)平均法計算面雨量值。圖2 給出了延河各時段的PFa-Ra和PFa-SSYa關(guān)系。
表3 研究區(qū)各流域P0 時段汛期降雨P(guān)Fa 與年地表徑流Ra、年產(chǎn)沙模數(shù)SSYa 的回歸關(guān)系式 Table 3 Regressions between Ra and PFa, and between annual specific sediment yield (SSYa ) and PFa for the reference period
圖2 延河流域甘谷驛站不同時段年汛期降雨量(PFa)與年地表徑流量(Ra)、年產(chǎn)沙模數(shù)(SSYa)的關(guān)系 Fig.2 Relationships between flood-season precipitation (PFa) and annual surface runoff (Ra) and annual specific sediment yield (SSYa) for different periods in the Yanhe basin
流域產(chǎn)沙主要發(fā)生在地表徑流時段,基流產(chǎn)沙很少。如忽略基流產(chǎn)沙,則任一流域任一時段的水沙關(guān)系總可以表示為:SSY=CR,其中C 和R 分別為該時段的地表徑流平均含沙量(kg/m3)和地表徑流量(mm/a)。在該水沙關(guān)系框架下,任一因素對SSY 的影響機(jī)制可概括為3 類:1)僅通過影響R 來影響SSY;2)僅通過影響C來影響SSY;3)通過同時影響C 和R 來影響SSY。黃土高原的水保措施主要包括植被、梯田和淤地壩3 類,其中植被和梯田為坡面措施,淤地壩為溝道措施。在以上水沙關(guān)系框架下,不同水保措施的減沙機(jī)制存有顯著差異。
對植被等坡面措施減沙機(jī)制的多尺度分析表明[9-10],在坡面尺度,植被等既減少R 值,也降低C 值,導(dǎo)致減沙率顯著高于減水率;但在流域尺度,植被等坡面措施減沙機(jī)制與坡面顯著不同。其原因為是黃土區(qū)溝谷重力侵蝕極其發(fā)育,而重力侵蝕基本上與植被無關(guān)[23]。在黃土區(qū),即使喬木,根系也主要分布在地表以下80 cm 深度范圍內(nèi)[24-27],而現(xiàn)代侵蝕溝(指切溝和沖溝)的深度一般為幾米到十幾米甚至幾十米,遠(yuǎn)超植被根系的作用范圍。由于植被等不能有效地抑制溝谷重力侵蝕,坡面水下溝后會繼續(xù)沖刷,這導(dǎo)致流域出口的C 值和水沙關(guān)系并沒有發(fā)生變化,洪水事件(文獻(xiàn)[9-10,17])和年尺度(圖3)均如此。因此,在流域尺度上,植被等坡面水保措施僅通過減R 來減沙,使得減沙率和減水率非常接近。目前文獻(xiàn)中支持該結(jié)論的觀測證據(jù)主要有:1)插財主溝和羊道溝1956—1970 年期間的對比觀測[9,10];2)橋子?xùn)|溝和橋子西溝1987—2006 年期間的對比觀測[17];3)甘肅天水站3 個水土保持治理流域(即橋子?xùn)|溝、呂二溝和羅玉溝)1980s—2010 年期間水沙關(guān)系隨流域治理的演變[16-17,28-29]。這3 個流域的長時間序列觀測數(shù)據(jù)均表明,如無淤地壩措施,植樹、修建梯田等坡面變化均未導(dǎo)致水沙關(guān)系的變化。如橋子?xùn)|溝1987—2006 年期間,流域坡耕地面積占總面積比例從54.2%下降到8.3%,植被和梯田面積比例從40.1%增長到90.5%,下墊面的巨大變化造成了流域降雨-產(chǎn)流和降雨-產(chǎn)沙關(guān)系的顯著改變,但并沒有改變流域的徑流-產(chǎn)沙關(guān)系,導(dǎo)致單位徑流輸沙量即C 值變化不大[16]。只有當(dāng)2006 年后流域內(nèi)建設(shè)了大量淤地壩后,水沙關(guān)系才有明顯改變,C 值才顯著降低[16]。Zhang 等[30]最近的計算也表明,在黃河中游,植被主要通過減少產(chǎn)流來減沙,而徑流含沙量的變化主要與淤地壩建設(shè)有關(guān)。
在黃土區(qū)大規(guī)模治理前,由于泥沙來源充沛,絕大部分水流都可以達(dá)到?jīng)_刷限制下的極限含沙量(水流此時不能沖刷,但輸沙能力未必達(dá)到飽和)[31-33]。其他發(fā)生在坡面的人類活動,如修路、耕作等,和植被一樣,也不能改變黃土區(qū)泥沙來源充沛的特點,因此也應(yīng)僅僅通過改變地表徑流量R 來影響流域產(chǎn)沙。這一特性在大流域應(yīng)更為突出,因為小流域受坡面過程影響顯著,隨流域面積增大,坡面過程對水沙輸移影響會逐漸減小[9]。在黃土高原,也有許多中大尺度流域,如岔巴溝流域[34]、大理河流域[35,36]、窟野河流域[9]、河龍區(qū)間[37]、渭河流域[9]等,報道有流域治理未改變徑流/輸沙關(guān)系或流域治理的減水率和減沙率接近的情況。
美國亞利桑那州Walnut Gulch 試驗流域也有和王家溝、天水試驗站類似現(xiàn)象[31]。該流域坡面小區(qū)的觀測表明,灌木覆蓋小區(qū)的平均含沙量為草地覆蓋小區(qū)的2 倍以上[38]。但在流域尺度,草地覆蓋流域比灌木覆蓋流域產(chǎn)沙模數(shù)偏小37.5%,產(chǎn)流模數(shù)偏小33.2%,兩者非常接近,這導(dǎo)致兩者的徑流平均含沙量也非常接近(1.34% vs. 1.36%)[31]。該流域泥沙來源也非常充沛[38],因此認(rèn)為,坡面變化僅通過改變R 來影響流域產(chǎn)沙,該結(jié)論可能是泥沙來源充沛地區(qū)的一個普遍規(guī)律。
圖3 羊道溝(非治理流域)和插財主溝(治理流域)年徑流深(Ra)與年產(chǎn)沙模數(shù)(SSYa)關(guān)系的比較 Fig.3 Comparison of the Ra-SSYa relationships between Yangdaogou (non-managed) and the Chacaizhugou (managed) watersheds
圖4 表明,在羊道溝小流域,降雨對洪水平均含沙量的影響很小。盡管低的含沙量基本發(fā)生在小雨量或小雨強(qiáng)情況下,但小雨量或小雨強(qiáng)時也可能出現(xiàn)很高的含沙量。對主要的產(chǎn)沙事件(定義為洪水平均含沙量高于500 kg/m3的事件,其產(chǎn)沙量占總量的86.9%),洪水平均含沙量不僅與雨量和雨強(qiáng)完全不相關(guān)(岔巴溝等流域也報導(dǎo)有類似現(xiàn)象[39-40]),而且變異較小,變差系數(shù)(CV)為16.9%。其他許多黃土流域也有類似的特點,這導(dǎo)致洪水產(chǎn)沙模數(shù)和徑流深關(guān)系可以用比例函數(shù)很好地擬合[41-42]。因此認(rèn)為,降雨對流域產(chǎn)沙的影響機(jī)制與植被相同,也僅通過影響R 來影響產(chǎn)沙。
淤地壩的減沙機(jī)制取決于其是否有良好的排水、排洪設(shè)施。無排水設(shè)施的淤地壩,如悶葫蘆壩,如未發(fā)生漫溢和垮壩,會全部攔蓄上游來水來沙。在黃土區(qū)大規(guī)模治理前,C 值不僅歷次洪水之間變化不大,空間上變異也很小[31-32],如在無定河黃土區(qū),C 值的空間變差系數(shù)僅為7%左右[32]。因此可認(rèn)為,當(dāng)流域某一范圍的來水來沙全部為無排水設(shè)施的淤地壩攔截時,流域出口的含沙量并不會改變,無排水設(shè)施的淤地壩僅通過減R 來減沙。有良好排水設(shè)施(包括溢洪道和放水建筑物)的淤地壩,壩體的壅水和攔截效應(yīng)會使得大量泥沙落淤,導(dǎo)致排放到下游的洪水含沙量大大變小。此類淤地壩如蓄水,則通過減R 和減C 兩種機(jī)制減沙,如不蓄水,則僅通過減C來減沙。這幾類淤地壩在黃土區(qū)流域一般均存在,因此可認(rèn)為淤地壩既通過減R 來減沙,也通過減C 來減沙。
圖4 羊道溝流域洪水事件平均含沙量(Ce) 隨降雨量(Pe)和最大30min 雨強(qiáng)(I30)的變化 Fig.4 Mean sediment concentrations of individual runoff events (Ce) varying with rainfall amount (Pe) and maximum 30-min intensity (I30) in the Yangdaogaou watershed
在黃土高原流域治理早期,大部分淤地壩排水設(shè)施配置均較差,淤地壩的減沙機(jī)制會以減R 為主;后期有良好排水設(shè)施的淤地壩數(shù)量逐漸變多,因此通過減C 而實現(xiàn)的減沙效應(yīng)會逐漸增強(qiáng)。圖2a 比較了延河流域P0、P1、P2 時段的降雨徑流關(guān)系,結(jié)果表明,相同降雨情況下P1 和P2 時段的R 值顯著小于P0 時段,R 值的降低顯然與淤地壩的攔蓄作用密切相關(guān);但3 個時段的C 值分別為460、469、341 kg/m3(表2),這表明只有到了P2 時段,淤地壩才通過降低C 值這種方式來減沙,清澗河流域有相同趨勢。在大理河流域,盡管C 值在P1時段就已開始降低,但降幅有限,僅在P2 時段才開始顯著降低(表2)。1980s 期間新建淤地壩數(shù)量較少[43],出現(xiàn)這種情況的原因很可能是由于黃河中游1977—1978 年大洪水期間,發(fā)生了嚴(yán)重的垮壩現(xiàn)象,為防止類似現(xiàn)象再發(fā)生,其后對許多淤地壩補建了溢洪道等排水設(shè)施[44]。
圖5 表明,在P0 時段,和事件尺度(圖4)類似,主要輸沙年份的地表徑流平均含沙量年際之間變異也較小,CV 為15%左右,這些主要輸沙年份貢獻(xiàn)了約90%的總產(chǎn)沙。假定某一流域在P0 時段的地表徑流平均含沙量為C0,根據(jù)以上分析可知,除非該流域后期修建有較多的排水良好的淤地壩,則該流域地表徑流平均含沙量可始終估計為C0。
圖5 P0 時段主要輸沙年份的地表徑流平均含沙量(Ca) Fig.5 Mean sediment concentrations of surface runoff (Ca) among the major sediment-yield years in P0 period
給定某一黃土區(qū)流域,假定多年平均降雨、地表徑流量、地表徑流平均含沙量和產(chǎn)沙模數(shù)在無治理期P0 時段的分別為P0、R0、C0、SSY0,流域治理后某一時段Pi分別為Pi、Ri、Ci、SSYi,設(shè)ΔP=Pi-P0,ΔR=Ri-R0,ΔC = Ci-C0,ΔSSY=SSYi-SSY0。忽略水土保持治理以外的其他人類活動,并假定流域治理措施僅包括坡面措施和淤地壩,則可認(rèn)為:ΔR=ΔRP+ΔRSM+ΔRCD,ΔSSY=ΔSSYP+ ΔSSYSM+ΔSSYCD,其中ΔRP、ΔRSM、ΔRCD分別為降雨、坡面、淤地壩對ΔR 的貢獻(xiàn),ΔSSYP、ΔSSYSM和ΔSSYCD分別為降雨、坡面措施和淤地壩對ΔSSY 的貢獻(xiàn)。我們的目標(biāo)為:根據(jù)P0 和Pi 時段的降雨、徑流和產(chǎn)沙觀測資料,確定ΔSSYP、ΔSSYSM和ΔSSYCD。
可用R、C、SSY 作為描述流域水沙狀態(tài)的3 個變量。在P0 時段,流域狀態(tài)為(R0,C0,SSY0),由于降雨變化和后期流域治理的實施,在Pi 時段流域狀態(tài)變?yōu)椋≧i,Ci,SSYi)。按第三部分分析,降雨、坡面治理和淤地壩均可影響R,但僅淤地壩會影響C 值。流域狀態(tài)從(R0,C0,SSY0)變?yōu)椋≧i,Ci,SSYi)可分解為以下4 個步驟:
1)降雨變化ΔP,導(dǎo)致R 變化了ΔRP,由于降雨不改變C 值,故C 值仍等于C0,經(jīng)過此步驟,SSY 變化了C0ΔRP;
2)坡面措施進(jìn)一步使R 變化ΔRSM,同上,C 值仍等于C0,經(jīng)過此步驟,SSY 變化了C0ΔRSM;
3)含沙量為C0的水流到達(dá)淤地壩,其中有ΔRCD的徑流被淤地壩攔截,經(jīng)過此步驟,SSY 變化了C0ΔRCD;
4)排出淤地壩到達(dá)流域出口的徑流量變?yōu)镽i,這些徑流在匯入淤地壩時含沙量為C0,到達(dá)流域出口時為Ci,經(jīng)過此步驟,SSY 變化值為RiΔC,完成從狀態(tài)(R0,C0,SSY0)到(Ri、Ci、SSYi)的變化。
以上第一個步驟中SSY 的變化完全由降雨引起,故可認(rèn)為:
同理,第二個步驟中SSY 的變化可認(rèn)為是坡面措施的貢獻(xiàn),故有:
第三步驟中SSY 的變化是淤地壩通過減R 導(dǎo)致的貢獻(xiàn),第四個步驟中SSY 的變化是淤地壩通過改變C 值導(dǎo)致的貢獻(xiàn),故:
式(4)—(6)等式右端諸項之和剛好為ΔSSY。ΔSSY的分割公式因此為:
式中ΔSSYH=ΔSSYSM+ΔSSYCD,為人類活動對ΔSSY 的貢獻(xiàn)。根據(jù)公式(7),ΔSSYH可估算為:
按以上分析,如某一坡面措施m 使得R 變化了ΔRm,則其對ΔSSY 的貢獻(xiàn)為C0ΔRm。如欲進(jìn)一步區(qū)分不同坡面措施的貢獻(xiàn),則:
式中ΔSSYm為第m 種坡面措施的貢獻(xiàn),∑ΔSSYm即ΔSSYSM,為各坡面措施貢獻(xiàn)之和。
以上公式中,計算結(jié)果中的負(fù)值表示減沙效應(yīng),正值表示增沙效應(yīng)。式(6)中,如ΔC>0,這表示淤地壩修建后并沒有改變C 值,而非淤地壩具有增沙效應(yīng)。實際計算時,如出現(xiàn)ΔC>0,可置ΔC=0,并將SSY 的變化全部歸因于R 的變化。
本部分將利用以上提出的方法(以下簡稱水沙關(guān)系法),計算各研究流域的減沙效益,并對結(jié)果進(jìn)行驗證,表4 給出了延河、清澗河、大理河流域的計算結(jié)果。驗證方案分兩種:1)比較水文法和水沙關(guān)系法計算的ΔSSYH是否一致,以此驗證公式(8)和(4);2)通過與觀測值比較驗證公式(4)—(6)。
分別按水沙關(guān)系法和水文法,計算了表2 中各流域各治理時段的ΔSSYH。水文法被認(rèn)為是一種較可靠的貢獻(xiàn)分割方法[5-6],圖6 表明,水文法和水沙關(guān)系法計算的結(jié)果非常接近。這證明了公式(8)的有效性,也間接證明了公式(4)的有效性。
對某一特定流域,如下墊面無變化,則歷年產(chǎn)沙的變化可認(rèn)為完全是降雨的變化所導(dǎo)致?;诖思僭O(shè),利用羊道溝的數(shù)據(jù)繼續(xù)對公式(4)進(jìn)行了驗證,結(jié)果見圖 7。ΔSSYP的觀測值和計算值非常接近。
圖8 利用羊道溝/插財主、橋子?xùn)|溝/橋子西溝兩對對比流域的觀測數(shù)據(jù)對公式(5)進(jìn)行了驗證。研究時段內(nèi),治理流域插財主溝和橋子?xùn)|溝均只有坡面治理措施,而無淤地壩。在此情形下,治理流域和非治理流域產(chǎn)沙模數(shù)的差異即為ΔSSYSM的觀測值(真實值)。圖8 表明,在多年時間尺度,公式(5)計算值對觀測值的偏差非常小,兩個地點都僅約6%;在年時間尺度,計算精度有所下降,但計算值和觀測值的R2仍高于0.8。
表4 水沙關(guān)系法計算結(jié)果 Table 4 Results of runoff-sediment relationship method
圖6 水文法和水沙關(guān)系法計算的人類活動減沙量比較 Fig.6 Comparison of the calculated sediment-reduction effects of human activities between the H (hydrology) method and the RSYR (runoff-sediment yield relationship) method
圖9 計算了表2 中各流域除P0 外其他時段的淤地壩攔沙量,按忽略(圖9a)和考慮(圖9b)坡面措施減水量兩種方式進(jìn)行計算,并分別與表2 中給出的實測值進(jìn)行了比較。當(dāng)考慮坡面措施減水量時,計算值和實測值非常吻合。然而,即使忽略坡面措施減水量,計算值和實測值也偏差不大,其原因是在研究時段,流域減沙效益以淤地壩為主, 坡面措施貢獻(xiàn)有限。根據(jù)表4 計算結(jié)果,對P1、P2、P3 時段,各流域坡面措施對減沙貢獻(xiàn)為8.4%~17.9%,平均為12.4%,而淤地壩為78%~101%,平均為86.9%。圖 9a 中兩者的偏離有隨淤地壩攔沙量變大而變大的趨勢,主要原因是淤地壩攔沙量較大的一般對應(yīng)大的流域,而大流域的坡面措施面積一般也大,其減沙量的絕對值也大,忽略坡面措施的貢獻(xiàn)因此在大流域會造成更大的絕對誤差。一旦考慮了坡面措施的貢獻(xiàn),這種趨勢就完全消失。
圖7 利用羊道溝數(shù)據(jù)驗證降雨的減沙效應(yīng)計算公式(公式(4)) Fig.7 Verifying the calculated effect of rainfall on sediment yield by equation (4) using annual observations of the Yangdaogou (non-managed) watersheds
圖8 利用對比觀測流域數(shù)據(jù)驗證坡面水保措施的減沙效應(yīng)計算公式(公式(5)) Fig.8 Verifying the calculated effect of slope measures on sediment yield by equation (5) using the paired-watersheds observations
圖9 利用淤地壩攔沙量實測數(shù)據(jù)驗證淤地壩的減沙效應(yīng)計算公式(公式(6)) Fig.9 Verifying the calculated effects of check dams on sediment yield by equation (6) using the measured values of the trapped sediment by check dams
本文淤地壩攔沙效益的計算公式,主要在多年時間尺度進(jìn)行了驗證,但也有兩次年時間尺度的結(jié)果,即曹坪站1977 和1978 年。1977 和1978 年黃河中游發(fā)生特大暴雨,估算的曹坪站淤地壩攔沙量與實測值的偏差僅為4.2%和11.1%(表4)。1978 年汛期降雨為487 mm,遠(yuǎn)超P0 時段多年均值(296 mm),計算結(jié)果表明,降雨導(dǎo)致產(chǎn)沙模數(shù)增加了23 163 t/km2,淤地壩攔沙導(dǎo)致其減少了22 784 t/km2,兩者大致相等,因此曹坪站該年的產(chǎn)沙模數(shù)(21 000 t/km2;表2)保持和P0 時段的平均值(21 901 t/km2)接近。但這并不意味著本文方法也能運用于枯水年,因為枯水年含沙量低的原因,基本不是由于淤地壩的攔沙作用,而是因為地表徑流流量很小,侵蝕和搬運泥沙的能力較弱的緣故。坡面措施減沙效益計算公式也不適用于枯水年,其原因是枯水年以小洪水或小流量為主,離石和天水站的觀測均表明,流域進(jìn)行坡面治理后,小洪水或小流量情況下含沙量明顯變小[10,17]??傊?,坡面措施和淤地壩減沙效益計算公式均不適用于枯水年。但由于枯水年對流域產(chǎn)沙貢獻(xiàn)很小,故本文方法能適用于多年時間尺度。
如上述,坡面治理可顯著降低小洪水或小流量情況下的含沙量,因此植被等坡面措施的減沙機(jī)制可更精確地概括為“植被等坡面措施不改變流域主要輸沙時段的水沙關(guān)系”。當(dāng)流域治理達(dá)到很高程度,以至于地表徑流以小洪水或小流量為主時,植被等坡面水保措施的減沙機(jī)制也可能發(fā)生變化。在黃土高原未治理前,當(dāng)流量小于某一閾值時,同流量下的含沙量變幅極大,當(dāng)流量超過該閾值時,含沙量的變幅明顯縮小[45],趨近于沖刷限制下的極限含沙量[31-33]。如坡面治理能將洪峰流量都控制在小于該流量閾值范圍內(nèi),則洪水的含沙量會有非常顯著的降低。但這一閾值非常小(參考文獻(xiàn)[41]圖2),將主要輸沙時段的流量控制在該閾值以下會非常困難。
坡面措施減沙效益的計算公式(公式(5))只在小流域尺度進(jìn)行了驗證。但大的河流,其水沙特性會比小的河流更難以改變,若流域治理未能改變小流域的C 值,則可以推斷,流域治理更不可能改變大流域的C 值[9]。因此認(rèn)為,公式(5)也應(yīng)該可以很好地運用于大流域的植被減沙效益計算。
2000 年以后,隨著退耕還林政策的實施,黃河中游植被覆蓋度顯著提高,同時,也新修了大量排水良好的淤地壩。本文對淤地壩減沙效益計算公式(公式(6))的驗證,未包括2000 年以后的數(shù)據(jù),因此需要進(jìn)一步驗證。任何一種治理措施,其減沙效益均通過改變R 或C值實現(xiàn),因此按公式(6)計算結(jié)果應(yīng)包含淤地壩的全部減沙效益。除通過壩體直接攔沙外,淤地壩尚有其他的減沙途徑,如淤積泥沙抬高了侵蝕基準(zhǔn)面,從而能穩(wěn)定坡谷、減少重力侵蝕的發(fā)生。退耕政策實施前,溝谷中泥沙淤積厚度可能尚有限,淤地壩的減沙效益主要通過壩體攔沙來實現(xiàn),故本文方法計算結(jié)果和實測攔沙量非常接近。但隨著淤積年限的增加,通過穩(wěn)定坡谷這一途徑實現(xiàn)的減沙效益應(yīng)該會越來越重要,這可能導(dǎo)致公式(6)計算結(jié)果顯著超過實測攔沙量,兩者間的差異可認(rèn)為是淤地壩通過穩(wěn)定坡谷而實現(xiàn)的減沙效益。
基于降雨和不同水保措施的減沙機(jī)制差異,本文提出了一種新的流域減沙貢獻(xiàn)分割方法。驗證結(jié)果表明,該方法可較好進(jìn)行降雨、不同坡面措施以及淤地壩等的減沙貢獻(xiàn)分割。常規(guī)的水保法基于坡面小區(qū)觀測結(jié)果,其計算結(jié)果為減蝕率,即水保措施在坡面減少了多少侵蝕。本方法考慮了減沙效益的空間尺度變異,計算結(jié)果為減沙率,即水保措施在流域出口減少了多少泥沙。本方法有較好的理論基礎(chǔ),可以有效地估算氣候變化和不同流域治理措施或不同人類活動對入黃泥沙的影響,對厘清黃河泥沙銳減的主導(dǎo)驅(qū)動因素和未來治黃策略的制定具有重要意義。