王奇花,高玉鳳,田 沉,李夢瑤,付興濤
(1.太原理工大學 水利科學與工程學院,山西 太原 030024;2.山西省水土保持科學研究所,山西 太原 030013)
晉西梁峁交錯、溝壑縱橫、呈現(xiàn)支離破碎的地貌景觀,加之溝坡陡峭、土壤特殊(離石黃土以粉砂為主,質(zhì)地疏松,遇水易溶解、崩塌)、植被稀疏、夏季多短歷時大強度暴雨,極易形成地表徑流并造成嚴重的水土流失。以往晉西的研究主要采用小流域徑流場、徑流小區(qū)實測等方法且多集中于侵蝕程度的時段性變化、溝坡侵蝕方式演化及水土保持措施的效益分析等方面[1-2],就土壤侵蝕風險評估方面研究相對較少,而準確預測產(chǎn)流產(chǎn)沙率,從而選擇合理的水土保持措施對該區(qū)水土流失治理尤為重要。
土壤侵蝕模型作為預報水土流失、指導水土保持措施配置的有效工具[3],是土壤侵蝕學科的前沿領(lǐng)域和土壤侵蝕過程定量研究的有效手段。國外許多土壤侵蝕模型被運用于黃土高原土壤侵蝕過程研究,國內(nèi)學者也在總結(jié)國外土壤侵蝕模型研究成果的基礎(chǔ)上,開發(fā)出了有針對性的土壤侵蝕模型??傮w來看,CSLE、Zheng 等經(jīng)驗統(tǒng)計模型[4-5]模擬精度雖高,但不能對侵蝕過程做出理論性解釋;物理成因模型DYRIM[6]可以進行氣候變化與土地情景分析,但輸入數(shù)據(jù)復雜;Yang等[7-8]考慮了溝道侵蝕,卻分別存在不考慮河道淤積和預測精度不高且模擬結(jié)果不穩(wěn)定的問題;WEPP和EUROSEM模型能夠較好地長時間預測徑流泥沙量[9-10],但Centeri等[11]通過將WEPP、EUROSEM和MEDRUSH模型模擬結(jié)果與人工模擬降雨試驗實測數(shù)據(jù)對比發(fā)現(xiàn),3個模型的有效性均隨土壤基本性質(zhì)的變化而變化,且統(tǒng)計分析表明EUROSEM較WEPP與MEDRUSH模型在估計土壤流失方面較為突出。
EUROSEM和WEPP模型均是在坡面徑流小區(qū)觀測資料基礎(chǔ)上發(fā)展起來的,主要用于模擬和預測坡面和小流域的水土流失。已有研究顯示W(wǎng)EPP模型在晉西地區(qū)具有較好的適用性[12],這類連續(xù)模擬模型需要一年中大量有關(guān)氣候和土地使用條件變化的輸入數(shù)據(jù),而黃土高原地區(qū)水土流失主要是由少數(shù)幾次大雨或暴雨所引起的[13],雖然該模型也可以針對次降雨事件,但只能模擬次降雨的總土壤流失量?;谑录膭討B(tài)分布式土壤侵蝕模型——EUROSEM模型,則是將單次降雨初始條件指定為數(shù)據(jù)輸入,以分鐘為基礎(chǔ)對水沙峰值的排放量大小和出現(xiàn)時間進行預測[10],這與該區(qū)的降雨特征和需求相契合。國內(nèi)外學者運用EUROSEM模型進行了諸多模擬研究,如王宏等[14]應用該模型對三峽庫區(qū)陡坡地侵蝕狀況進行了模擬,指出其在預測產(chǎn)流量上效果較好,而在產(chǎn)沙量的模擬中效果相對較差;Folly等[15]與Mati等[16]分別在荷蘭和肯尼亞的流域上分析此模型的適用性,結(jié)果均表明其能夠較好地模擬不同環(huán)境和降雨特征下的產(chǎn)流率峰值和總產(chǎn)流量;該模型在尼日利亞裸土表面進行的降雨模擬試驗也得到了良好的模擬結(jié)果[17]。該模型在產(chǎn)流時間和產(chǎn)沙率的模擬方面存在一定不足,這主要是由于其對徑流含沙量的估計不足[18],但在大多數(shù)情況下,降雨事件期間的侵蝕產(chǎn)沙總量是可以充分預測的。
本文將該模型引用到晉西有較好代表性的王家溝流域[19],為該地區(qū)土壤侵蝕風險評估提供理論依據(jù)。鑒于此,本文將野外人工模擬降雨試驗與模型模擬方法相結(jié)合,分析晉西黃綿土坡面徑流侵蝕過程及其與影響因素的相關(guān)關(guān)系,通過對比分析實測數(shù)據(jù)與模擬結(jié)果,探討該模型在本地區(qū)的模擬效果,并評價其在不同雨強、坡長條件下對次降雨產(chǎn)流、產(chǎn)沙的預報能力,為本地區(qū)土壤侵蝕風險評估及水土流失治理提供理論與應用依據(jù)。
2.1 研究區(qū)概況試驗于2019年7—9月在山西省呂梁市離石區(qū)王家溝流域野外徑流小區(qū)開展。地理坐標為東經(jīng)111°11′,北緯37°31′,氣候?qū)贉貛Т箨懶约撅L氣候,年平均氣溫8.9℃。流域最大年降雨量711.50 mm,最小年降雨量240.20 mm,多年平均降雨量490.30 mm,年內(nèi)降雨分布不均,7—9月份雨量占全年的62.3%。試驗區(qū)土壤母質(zhì)為新生界第四系中更新統(tǒng)離石黃土上或晚更新統(tǒng)馬蘭黃土與中更新統(tǒng)離石黃土的混合土,是典型的砂質(zhì)壤土,其土壤顆粒組成為砂粒84.05%、粉砂粒14.20%、黏粒1.75%,有機質(zhì)含量約13.42 g/kg,pH值8.15,總孔隙度49.05%。
2.2 試驗設(shè)計根據(jù)王家溝暴雨資料,流域內(nèi)高強度、短歷時暴雨多發(fā)生于7—9月,暴雨強度主要集中于60~90 mm/h,是造成該區(qū)強烈土壤侵蝕的主要原因;另外,山西省水文局降雨監(jiān)測資料顯示,該研究區(qū)汛期最大降雨強度可達90.30 mm/h,因此本試驗設(shè)計雨強為60、90、120 mm/h。試驗徑流小區(qū)坡度為20°,寬度為2 m,坡長分別為2、3、4 m,表面均無植被覆蓋。降雨器噴頭距離地面高度10 m 左右,根據(jù)付興濤[20]對降雨均勻性的測定及雨強標定,得出本試驗所用人工模擬降雨器的降雨均勻系數(shù)在85%以上,雨滴分布及終速等指標均符合試驗要求,可以開展試驗。在每個徑流小區(qū)進行3場不同強度的降雨試驗,共降雨9場次。
每場降雨試驗開始前,測定坡面土壤體積含水率約為25%。降雨開始后,用秒表記錄下產(chǎn)流時刻。開始產(chǎn)流后,每隔2 min 用采樣桶采集一次泥沙樣;自產(chǎn)流開始持續(xù)降雨30 min,共采集15個徑流泥沙樣。降雨結(jié)束后測量采集的渾水體積以得到徑流量,并將采樣桶靜置,直至泥沙全部沉淀,倒去上部清水,將泥沙烘干(105℃條件下烘24 h)稱重得到產(chǎn)沙量。
圖1 野外人工模擬降雨試驗
2.3 EUROSEM模型EUROSEM模型是一個基于過程的單一事件模型,通過對土壤侵蝕過程的物理描述,以分鐘為時間單位模擬次降雨條件下地塊或小流域侵蝕過程[10]。它主要涉及植被對降雨的截留、下滲、雨滴濺蝕、徑流侵蝕、徑流搬運和泥沙沉積。該模型由模塊化結(jié)構(gòu)組成,可以和地理信息系統(tǒng)進行無縫鏈接,它將自身鏈接到KINEROS模型的水沙運動結(jié)構(gòu)中來模擬侵蝕,水沙通過一系列相互連接的均勻斜坡面和溝道要素在地表運動,其中要素特征可參數(shù)化。在模型運算中,通過參數(shù)設(shè)置來描述地形及降雨特征。
試驗在晉西黃綿土坡面進行,故只考慮坡面要素,從3個坡長的徑流小區(qū)中各選擇一個來設(shè)定EUROSEM模型的參數(shù)。土壤比重(RHOS)通過比重瓶法測定,并結(jié)合環(huán)刀取土計算土壤孔隙度(POR),初始土壤含水率(THI)通過烘干稱重獲得,而雨滴沖擊土壤顆粒可分散性(EROD)、土壤聚合度(COH)、坡面糙率(RFR)、入滲滯后因子(RECS)則參考Morgan等[10]給出的參考值,飽和導水率(FEIN)、毛細管張力(G)、曼寧系數(shù)(MANN)通過不斷改變輸入?yún)?shù)數(shù)值并將模擬結(jié)果與實測值進行對比而確定(見表1)。
2.4 數(shù)據(jù)分析計算與模擬精度評價方法試驗采用Excel 處理相關(guān)數(shù)據(jù)并繪制產(chǎn)流產(chǎn)沙率隨產(chǎn)流歷時的變化曲線;采用SPSS 對雨強、坡長與產(chǎn)流產(chǎn)沙總量進行相關(guān)性分析;采用Nash模型效率系數(shù)ME[21]和相對誤差RE評價模型模擬精度。模型效率系數(shù)ME越高,說明實測值與模擬值擬合效果越好;而相對誤差RE越小,說明實測值與模擬值擬合度越高,模型適用性越好。
表1 EUROSEM模型主要參數(shù)值
3.1 產(chǎn)流分析坡度為20°,不同雨強和坡長條件下,產(chǎn)流率隨降雨時間的延長總體上呈增大趨勢,在開始產(chǎn)流后的5 min內(nèi)增速很快,達到產(chǎn)流率峰值的77%~94%,后逐漸變緩并基本趨于穩(wěn)定,模擬結(jié)果與實測結(jié)果具有相似的變化趨勢(圖2)。分析原因,降雨初期土壤入滲強度大于降雨強度,降雨全部下滲,不產(chǎn)生地表徑流。隨著降雨的持續(xù),土壤含水率不斷增大,表土孔隙也在雨滴的擊濺作用下發(fā)生改變,土壤入滲能力迅速減小至低于降雨強度,坡面開始產(chǎn)流,且產(chǎn)流率迅速增大。在產(chǎn)流開始約5 min后,由于土壤含水率及表土孔隙變化減小,土壤入滲率基本趨于穩(wěn)定,產(chǎn)流率也相應趨于穩(wěn)定;此外,隨著降雨強度的增大,坡面單位時間所承受的雨量增大,雨滴動能也增大,土壤入滲率達到相對穩(wěn)定的時間縮短,所以一定坡長條件下降雨強度越大,產(chǎn)流率越大且達到峰值所需時間越短,孫佳美等[22]也指出降雨強度越大,達到產(chǎn)流穩(wěn)定狀態(tài)的時間越短。
進一步分析次降雨產(chǎn)流量隨坡長和雨強的變化可知,產(chǎn)流量隨坡長和降雨強度的增大而增大。在降雨強度為60 mm/h時,坡長從2 m延長至4 m,產(chǎn)流量增量為0.106 m3,而在降雨強度為90和120 mm/h時,其增量分別為前者的1.281和1.401倍,表明坡長對坡面產(chǎn)流量的影響隨著降雨強度的增大而增大。而降雨強度一定時,產(chǎn)流量隨坡長的延長不呈比例增加,坡長從2 m延長至3 m,3 m延長至4 m。產(chǎn)流量增量在降雨強度為60 mm/h 時為0.046、0.060 m3;90 mm/h 時為0.073、0.062 m3;120 mm/h 時為0.086、0.063 m3,說明降雨強度越大,坡長的延長對產(chǎn)流量增量的影響越顯著。分析其原因,隨著坡長的延長,承雨面積增大,單位時間內(nèi)坡面承雨量增加,坡面下部徑流流速加快,徑流下滲機會減少,隨著降雨的進行觀察到坡面下部細溝的出現(xiàn),導致坡面徑流在短時間內(nèi)匯集于細溝內(nèi),流出出口斷面。而隨著坡長延長,降雨強度增大,使得降雨擊濺表土的機會增加,濺蝕力增強,表土孔隙改變,表土結(jié)皮經(jīng)歷周期性發(fā)展[23],而結(jié)皮的形成能在很大程度上降低土壤入滲率[24],增大坡面產(chǎn)流量,且雨滴動能越大,結(jié)皮硬度越大[25]。
圖2 產(chǎn)流過程模擬值與實測值對比
相關(guān)性分析表明(表2),產(chǎn)流量與降雨強度呈極顯著正相關(guān)關(guān)系,相關(guān)系數(shù)為0.948,而坡長與其相關(guān)系數(shù)僅為0.279,說明在降雨強度、坡長對產(chǎn)流量共同影響下,降雨強度與產(chǎn)流量的相關(guān)性較坡長大。在分別剔除坡長、降雨強度變量的影響時,降雨強度和坡長與產(chǎn)流量均呈極顯著正相關(guān)關(guān)系,偏相關(guān)系數(shù)分別為0.987 和0.878,說明兩個變量共同作用時彼此制約了對產(chǎn)流量的影響,降雨強度很大程度上掩蓋了坡長對產(chǎn)流量的影響程度。該結(jié)論與付興濤等[26]通過室內(nèi)人工模擬降雨試驗研究黃土陡坡產(chǎn)流量隨坡長的變化過程得出的結(jié)論吻合。
3.2 產(chǎn)沙分析分析不同雨強和坡長條件下實測坡面產(chǎn)沙率隨產(chǎn)流歷時的變化可知(圖3),產(chǎn)流初期產(chǎn)沙率迅速增加,并在10 min內(nèi)達到第一次峰值,之后隨著產(chǎn)流歷時延長在某一數(shù)值附近波動變化,而模型模擬結(jié)果則顯示產(chǎn)沙率在波動增長后呈穩(wěn)定趨勢??赡苡捎诋a(chǎn)流初期雨滴直接打擊土壤表面,分散剝離表層松散物質(zhì),使其隨徑流流出坡面出口處,而隨著降雨時間的延長,產(chǎn)流量增加,增強了其對坡面的沖刷能力,使得產(chǎn)沙率急劇增大。但隨著降雨的持續(xù)進行,產(chǎn)流量和徑流挾沙能力逐漸達到穩(wěn)定,且薄層水流厚度的增加和坡面結(jié)皮的形成也有效緩解了雨滴的濺蝕及徑流對表土顆粒的沖刷,因此產(chǎn)沙率逐漸趨于穩(wěn)定。然而,細溝的產(chǎn)生會引起產(chǎn)沙率的急劇增大,并表現(xiàn)出一定的波動性[27]。此外,在整個徑流過程中,水流的能量分配是不斷變化的,溝床泥沙和被搬運泥沙在水流作用下不斷發(fā)生交換,從而使得侵蝕和沉積過程交替進行[28],因此產(chǎn)沙率并不表現(xiàn)出平滑趨勢而是在某一數(shù)值附近波動,這一數(shù)值受降雨強度、坡長等因素影響。
坡長和雨強是影響坡面徑流侵蝕產(chǎn)沙的重要因素[29-30],實測坡面產(chǎn)沙量顯示,產(chǎn)沙量隨坡長的延長和降雨強度的增加急劇增加,王占禮等[31]研究黃土裸坡土壤侵蝕過程指出,坡度為20°時不同降雨強度下產(chǎn)沙總量與坡長具有顯著的冪函數(shù)關(guān)系。在60、90及120 mm/h這3種降雨強度下,隨著坡長的延長,產(chǎn)沙量呈增長趨勢,但其增量存在減小情況,坡長從2 m延長至3 m,產(chǎn)沙量分別增加了0.226、1.512、2.788 kg;而從3 m 延長至4 m 時,則增加了0.250、1.510、2.040 kg,說明雨強大于60 mm/h時,坡長由2 m延長至3 m產(chǎn)沙量明顯增大,而由3m延長至4m時增量有所減小。分析其原因,隨著坡長的延長,侵蝕面積增大,可供濺蝕的物質(zhì)增多,且徑流所具有的重力勢能也增大,轉(zhuǎn)換成的動能相應增大,即侵蝕動力增大。此外,細溝的形成和發(fā)展直接影響著坡面產(chǎn)沙過程[32],而隨坡長延長,細溝侵蝕加劇[33],所以坡面產(chǎn)沙量必然隨坡長的延長而增大。但隨著坡長延長,徑流含沙量增加,水體能量主要消耗于泥沙的搬運,徑流參與侵蝕的能力減小,所以產(chǎn)沙量增量會減小。隨坡長延長產(chǎn)流量也表現(xiàn)出這一規(guī)律,結(jié)合晉西黃土坡面室內(nèi)模擬試驗結(jié)論[20],推斷4 m坡長為該研究區(qū)產(chǎn)流產(chǎn)沙量增量減小的臨界坡長。初步建議,以4 m為間隔布設(shè)水土保持措施,以減緩坡面水土流失。當降雨強度為60 mm/h,坡長由2 m延長至4 m時,坡面產(chǎn)沙量增量為0.476 kg;當降雨強度為90和120 mm/h時,坡面產(chǎn)沙量分別增長到了60 mm/h時的6.348和10.122倍。除了雨滴能量這一因素外,降雨強度越大,坡面產(chǎn)流量越大,挾沙能力越強,而且降雨強度的增大使得徑流紊動性增強[34],侵蝕能力相應增大。各種因素相疊加,導致了坡面產(chǎn)沙量隨降雨強度的增大顯著增加。
表2 產(chǎn)流產(chǎn)沙總量與雨強、坡長的相關(guān)性分析
圖3 產(chǎn)沙過程模擬值與實測值對比
相關(guān)性分析表明(表2),產(chǎn)沙量與降雨強度呈顯著正相關(guān)關(guān)系,相關(guān)系數(shù)為0.747,與坡長之間的相關(guān)系數(shù)為0.558,而在分別剔除坡長、降雨強度變量的影響時,降雨強度和坡長與產(chǎn)沙量之間的偏相關(guān)系數(shù)分別為0.900 和0.838。說明降雨強度、坡長在對產(chǎn)沙量共同作用時,降雨強度與產(chǎn)沙量的相關(guān)性較坡長大,而此時兩者對彼此都有一定的制約,當排除一個變量的影響,研究另一個變量對產(chǎn)沙量的作用時,二者均對其有很大的影響。
3.3 EUROSEM模型的適用性評價圖2、圖3顯示EUROSEM模型可對產(chǎn)流產(chǎn)沙過程進行較為細致的模擬,且對產(chǎn)流過程的模擬效果較產(chǎn)沙好,而良好的模擬產(chǎn)流過程是模擬侵蝕產(chǎn)沙的基礎(chǔ)[10]。在產(chǎn)流過程模擬中,產(chǎn)流率峰值出現(xiàn)時間與實測值稍有偏差,且坡長為2 m、3 m時,實測值整體在模型模擬值附近波動變化。而在4 m時實測值幾乎均小于模型模擬值,但產(chǎn)流率達到峰值時其相對誤差RE在不同雨強條件下分別為3.50%、2.70%和9.85%。在產(chǎn)沙過程模擬中,產(chǎn)沙率雖然都呈增長趨勢,但產(chǎn)沙率峰值首次出現(xiàn)時間卻較實測值提前了約5 min,其主要原因可能是在模型中細溝侵蝕由預先設(shè)定好的參數(shù)進行模擬,而在實際產(chǎn)流過程中,細溝存在發(fā)育過程,該階段使得模型模擬結(jié)果與試驗結(jié)果產(chǎn)生了時間差。由于與模型建立的試驗土質(zhì)不同,實測產(chǎn)沙率峰值出現(xiàn)的時間間隔也要長于模擬值。總體而言,雨強一定的條件下,產(chǎn)流產(chǎn)沙率實測值與模擬值間的差值隨著坡長的延長而增大。
圖4 模擬值與實測值關(guān)系
在產(chǎn)流產(chǎn)沙過程模擬的基礎(chǔ)上,整體來看,EUROSEM模型模擬效果較好,在晉西黃土高原典型坡面上的應用比較成功。圖4進一步表明,坡面產(chǎn)流產(chǎn)沙的模擬值與實測值呈顯著線性關(guān)系(R2均約為0.984)。相比而言,該模型對產(chǎn)流量的模擬效果較產(chǎn)沙量好,其效率系數(shù)ME高達0.978,模擬值與實測值相對誤差RE范圍為-12.95%~9.04%,從產(chǎn)沙量的模擬效果來看,雖然ME為0.974,但其RE范圍為-14.72%~24.13%,以實際產(chǎn)沙量的20%為許可誤差,RE合格率為89%。另外,對比不同坡長下實測值與模擬值(圖5)可知,模型模擬值與相應實測值間差值總體上隨坡長的延長和雨強的增大而增大,這可能是因為模型無法應對降雨期間不斷變化的水文條件和土壤特性[15],而且坡長越長,坡面侵蝕程度在不同坡段差異越大,降雨不僅引起土壤濺蝕,還能在土壤表面形成結(jié)皮[35]。4 m坡長條件下,產(chǎn)流產(chǎn)沙總量模擬值幾乎均大于實測值,造成高估的原因可能是4 m為該研究區(qū)產(chǎn)流產(chǎn)沙量增量減小的臨界坡長,而模擬初期產(chǎn)流產(chǎn)沙率過高,后期模擬值幾乎仍略高于實測值(在雨強為90 mm/h時,產(chǎn)沙率后期波動較大,導致實際產(chǎn)沙總量大于模擬值)。
圖5 模擬值與實測值對比
由于坡面糙率、導水率、細溝形態(tài)等特征在不同場次降雨間可能發(fā)生變化,而在模型模擬過程中直接輸入修正后的參數(shù),導致模型參數(shù)輸入的精確性方面可能有所欠缺。但總體而言,效率系數(shù)ME與相對誤差RE均表明模型在模擬研究區(qū)產(chǎn)流產(chǎn)沙總量上效果良好,為精確模擬徑流侵蝕產(chǎn)沙提供了良好的基礎(chǔ)。
基于野外人工模擬降雨試驗實測數(shù)據(jù),應用EUROSEM模型對晉西黃綿土坡面徑流侵蝕產(chǎn)沙過程進行模擬,并對比分析實測結(jié)果與模擬結(jié)果,得出如下結(jié)論:(1)雨強、坡長與產(chǎn)流產(chǎn)沙總量均為正相關(guān)關(guān)系,當二者共同作用于產(chǎn)流產(chǎn)沙量時,雨強較坡長對其影響大(雨強、坡長-產(chǎn)流量的相關(guān)系數(shù)分別為0.948、0.279,雨強、坡長-產(chǎn)沙量的相關(guān)系數(shù)分別為0.747、0.558),而偏相關(guān)分析顯示二者單獨對產(chǎn)流產(chǎn)沙總量均有很大的影響(雨強、坡長-產(chǎn)流量的偏相關(guān)系數(shù)分別為0.987、0.878,雨強、坡長-產(chǎn)沙量的偏相關(guān)系數(shù)分別為0.900、0.838)。(2)雨強大于60 mm/h,坡長由3 m延長至4 m時,產(chǎn)流產(chǎn)沙實測增量較2 m延長至3 m減小,且降雨強度越大,減幅越大。因此,初步建議在該區(qū)以4 m為間隔布設(shè)水土保持措施,以減緩水土流失。在4 m坡長條件下,產(chǎn)流產(chǎn)沙總量模擬值幾乎均大于實測值。(3)產(chǎn)流率隨降雨時間的延長先增大后趨于穩(wěn)定,EUROSEM模型模擬結(jié)果與實測結(jié)果具有相似的變化趨勢,且兩者峰值出現(xiàn)時間稍有偏差;實測產(chǎn)沙率在產(chǎn)流初期迅速增加,后在某一數(shù)值附近波動變化,雨強越大、坡長越長,波動越明顯,而模型模擬結(jié)果則在波動增長后呈平穩(wěn)趨勢,且產(chǎn)沙率峰值首次出現(xiàn)時間較實測值提前了約5 min。(4)EUROSEM模型對坡面產(chǎn)流產(chǎn)沙總量的模擬效果良好,相對誤差RE在許可范圍之內(nèi),反映模型總體預測效果的模型效率系數(shù)ME=0.978,0.974。
總體而言,EUROSEM模型能較好的預測晉西次降雨坡面徑流侵蝕產(chǎn)沙情況,說明該模型在晉西黃綿土坡面上有較好的適用性。但本次試驗在裸坡面上進行,僅考慮坡長、雨強的影響,并未考慮模型中的溝道組分,也未涉及植被等的影響,今后的研究中應進一步優(yōu)化相關(guān)參數(shù),增強模型在該區(qū)的適用性。