岔巴溝流域汛期徑流模擬及地表產流特征分析

景珂星, 毛 歡, 宋進喜, 黃 鵬, 吳 瓊, 龐國偉

(1.西北大學 城市與環境學院, 西安 710127; 2.陜西省地表系統與環境承載力重點實驗室,西安 710127; 3.陜西省黃河研究院, 西安 710127)

地表徑流作為水文循環中不可或缺的基本組成部分，是流域水文過程演變的重要指標[1]。在全球氣候變化的大背景下，干旱半干旱地區短歷時暴雨強度和極端強降水日數急劇增加[2]，加上人類活動影響下流域土地覆被的空間差異性使得流域徑流事件和地表產流能力形成了高度的時空變異性[3]，造成了極端洪澇事件的頻繁發生[4]，人民的生命財產安全受到了極大的威脅。為了準確預測徑流流量，明確產流規律，流域水文模擬已成為國內外水文科學研究的重要方法。

中國黃土高原地區特殊的土壤條件以及降水變差大且集中的氣候特征造成了嚴重的水土流失[5]。從20世紀50年代末，國家開始在黃土高原采取了一系列的水土保持措施，如修建淤地壩、梯田和植樹造林等，使得水土流失得到有效控制，生態環境明顯改善[6]。但21世紀以來，隨著大量淤地壩淤滿，水保措施抵御高強度降水事件的能力下降[7]，當地居民的生命財產安全受到了嚴重的威脅。因此探討不同土地覆被條件下的地表產流特征有助于深入了解黃土高原地區降雨—徑流的變化趨勢和產流規律，可初步為探索流域產水產沙與水土保持服務的權衡關系奠定基礎，為推動黃河流域生態保護和高質量發展提供參考。

目前針對黃土高原徑流變化及其影響因素的研究，一類是通過徑流小區的降雨—徑流試驗，以控制變量法分析不同地形和降水條件影響下的產流情況[8-9]；另一類是基于流域尺度的水文模型進行徑流預測及變化的歸因分析[10-11]。然而，由于黃土高原地區流域水文條件時空異質性導致地表徑流及產流能力在不同季節和降雨事件內是不斷變化的[12]，小范圍的降雨—徑流試驗對物理機制考慮不足，對徑流產生和變化的空間規律描述也有一定局限。而SWAT(Soil and Water Assessment Tool)模型不但考慮物理機制且將一些物理機制難以描述的過程用“經驗模型”來刻畫，這使其廣泛應用于不同環境下徑流模擬研究。SWAT模型在黃土高原的徑流模擬也有廣泛的應用，多位學者利用SWAT模型分別在黃土高原汾河[13]、桑干河[14]、渭河[15]、北洛河流域[16]及整個黃土高原[17]的氣候變化及土地利用對徑流變化的影響進行定量分析；龐佼等[18]利用SWAT模型對黃土高原安家溝流域的月徑流進行模擬；Zhao等[19]在涇河流域通過SUFI_2算法對SWAT模型模擬徑流的參數進行不確定性分析。這些研究證明SWAT模型在黃土高原區的適用性較好，但多以SCS徑流曲線法作為產流方法，與黃土高原超滲產流[20]的產流方式存在較大的差異，且研究尺度多集中于月尺度和年尺度，日尺度徑流模擬效果不佳，基于日徑流事件的模擬相對缺乏。因此，本文以黃土高原岔巴溝流域為研究區域，采用與黃土高原超滲產流原理相似的Green-Ampt下滲法驅動SWAT模型進行日尺度的汛期(6—10月)徑流模擬，并以地表徑流系數為指標對降雨強度和前期土壤含水量對不同土地覆被地表產流的影響進行分析，以期推進黃土高原地區日尺度水文模擬的發展，并為揭示黃土高原不同土地覆被產流機制的變化提供參考。

1 研究方法與數據來源

1.1 研究區概況

岔巴溝流域位于陜西省榆林市子洲縣北部(109.5°—110.1°E,37.5°—37.8°N)，是大理河的一級支流(圖1)，岔巴溝河全長26.6 km，流域面積205 km2，曹坪水文站控制面積187 km2。河流最上游發源于劉新窯村，年平均徑流總量約1×107m3，年平均流量約為0.35 m3/s。岔巴溝流域呈對稱樹葉狀，其海拔范圍為900～1 284 m，平均海拔大約1 080 m，流域屬溝壑密度大的黃土丘陵溝壑區，土層顆粒較小且較為疏松。流域年平均降水量約為474 mm，降水量年內分配極為不均，約70%降雨集中在7—9月，多為歷時短、雨強大的暴雨；降水年際變化較大，極易引起山洪泥石流等自然災害，如岔巴溝2017年“7·26”特大洪水等[21]。目前政府已投入大量資金進行山洪災害的防控，但相關研究仍落后于治理工程。

圖1 岔巴溝流域概況

1.2 數據來源

1.2.1 水文氣象數據 本研究所用的徑流數據選取位于流域出口的曹坪水文站的汛期日徑流數據；降水數據主要選取流域內13個雨量站的汛期降水摘錄數據。降雨徑流數據來源于《黃河流域水文資料—黃河中游區上段(河口鎮—龍門)上下冊》，時間范圍為2006—2013年。由于流域內沒有氣象站，根據需要選取流域周邊綏德、橫山、榆林3個國家氣象站的逐日氣溫、蒸發、風速和相對濕度數據，由中國氣象科學數據共享服務網(http:∥data.cma.cn)提供。水文站和雨量站點流域內分布情況見圖1。

1.2.2 遙感數據 構建SWAT模型需要的遙感數據主要包括流域數字高程模型(DEM)、土地覆被、土壤類型等。本研究采用的流域數字高程模型是30 m分辨率的GDEM數據，由中國科學院地理空間數據云平臺(http:∥www.gscloud.cn)提供；土地覆被數據為2010年30 m分辨率的中國土地利用遙感監測數據，由中國科學院資源環境科學數據中心(http:∥www.resdc.cn)提供；土壤數據為1 km格網的HWSD(Harmonized World Soil Database)世界土壤數據庫[22]。

1.3 SWAT分布式水文模型

1.3.1 SWAT模型原理及結構 SWAT模型對流域的水文建模主要分為水文循環的陸地階段和匯流階段。SWAT可以直接從文件中讀取輸入數據，如降水、最高和最低氣溫等，從而在運行時生成模擬器。SWAT模型對水文循環的模擬基于水量平衡方程，其表達式為：

(1)

式中:SWt表示土壤最終含水量(mm);SW0表示第i天土壤初始含水量(mm);t表示時間;Rday表示第i天的降水量(mm);Qsurf表示第i天地表徑流(mm);Ea表示第i天的蒸散發量(mm);Wseep表示第i天從土壤剖面進入包氣帶的水量(mm);Qgw表示第i天基流回歸的水量(mm)。

本文以Green-Ampt下滲法代替SCS徑流曲線法作為SWAT模型的產流方法，Green-Ampt方程假設地表永遠存在多余水分進而預測下滲，模型假設在濕峰以上的土壤完全飽和，在濕峰處含水量急劇變化。Mein等[23]開發了一種方法來確定Green-Ampt方程的下滲蓄水時間，該方法被引入SWAT中作為模擬地表徑流的替代方法，需要用戶輸入時間尺度小于日的降水數據計算徑流。Green-Ampt Mein-Larson方程的下滲速率定義為：

(2)

(3)

式中:finf,t為時刻t的下滲速率(mm/h);Ke為有效水力傳導率(mm/h);Ψwf為濕潤鋒土壤基質勢(mm),是土壤孔隙率、砂粒百分比和黏土百分比的函數;Δθv為濕潤鋒土壤體積含水量的變化(mm/mm);Finf,t為時刻t累計下滲量(mm);Ksat為飽和導水率(mm/h)；CN為徑流曲線數。

1.3.2 SWAT模型構建 構建SWAT模型的前期準備為土壤和氣象屬性數據庫的建立。土壤數據庫包括岔巴溝流域各類土壤的物理和化學性質，土壤物理性質中土壤的顆粒組成在HWSD數據庫中直接獲取，并作為SPAW軟件的輸入數據計算土壤濕密度(SOL-BD)、土壤層有效持水量(SOL-AWC)和飽和導水率(SOL-K)等；土壤化學性質主要包括pH、電導率(SOL-EC)、有機碳(SOL-CBN)等，由HWSD土壤數據庫直接獲取。氣象數據庫包括氣象站氣象數據以及岔巴溝流域13個雨量站降雨數據，以.txt格式輸入模型。

通過基礎數據的準備和數據庫的構建，借助ArcSWAT 2012版本軟件建立岔巴溝流域SWAT模型。首先基于DEM數據生成河網，通過推薦閾值將岔巴溝流域劃分為31個子流域，并設置流域總出水口；其次通過疊加土地利用數據、土壤數據和坡度數據，將三者的最小閾值均設置為0%，保留流域所有具有相同土地利用類型、土壤類型和坡度的640個水文響應單元(Hydrologic Response Units)；最后將氣象數據導入模型，改變流域產流方法為Green-Ampt下滲法后運行SWAT模型，實現岔巴溝流域2006—2013年水文過程的模擬，模擬結果“output.hru”文件中的時間步長內HRU總降水量(PRECIP，mm)、產生的地表徑流量(SURQ_GEN，mm)和初始時段土壤剖面的含水量(SW_INIT，mm)等數據用于后續地表產流特征分析。

1.3.3 SWAT模型敏感性分析和參數率定 不同的參數對SWAT模型的模擬結果影響很大，因此在模型率定前需進行參數的敏感性分析，選取合適的參數以提升模擬效率。本研究利用SWAT-CUP中的SUFI-2(sequential uncertainty fitting Version2)作為優化算法進行參數敏感性分析和率定。參數敏感性分析結果采用t統計量(t-Stat)和顯著性指標p值(p-value)來衡量，t統計量絕對值越大、顯著性指標p值越接近0，參數敏感性越強[24]。表1列出了岔巴溝流域率定使用的23個參數及其詳細信息。

表1 參數含義及初始率定范圍

模型效率通過確定系數(R2)、Nash-Sutcliffe效率(ENS)和百分比偏差(PBIAS)第3個指標來進行評估。指標的計算公式如下[25]：

(4)

(5)

(6)

在3類指標中，ENS的范圍為-∞到1，較高的ENS值是首選；R2的取值范圍從0到1，數值越大表示模擬效果越好；PBIAS比較模擬輸出與觀測數據的平均趨勢。一般情況下，當R2，ENS>0.5且PBIAS≤±25%，則認為模型效果較好。

2 結果與分析

2.1 岔巴溝流域SWAT模型敏感性分析結果

為了得到模型相關各參數的敏感性信息，本研究基于SWAT-CUP軟件采用SUFI-2算法對模型23個有關參數進行了多次迭代，得到了各參數的敏感性分析結果和參數率定范圍，其中敏感性較高的參數有CN2，CANMX，SLSOIL，LAT_TTIME，SOL_AWC和SOL_K等，表2列出了敏感性參數的驗證范圍、t統計量和顯著性指標p值的詳細結果。敏感性參數同Gao等[26]在黃土高原涇河流域的敏感性分析結果相近。其中CN2(徑流曲線數)是對徑流影響最顯著的參數，CN2值越大表明流域下墊面不透水性越強，產生的徑流量越大，CN2的t統計值為-7.239，p值為0，是最敏感的參數；樹冠最大儲量CANMX影響降雨的樹冠截留量，值越大徑流量越小；SOSOIL為側向地下水流坡長，數值大小決定流域坡面特征對HRU徑流計算的影響程度；LAT_TTIME為橫向流動時間，受土壤水力性質影響很大，影響流域匯流速度；SOL_AWC和SOL_K是與土壤屬性相關的參數，其中SOL_AWC為土壤層有效含水量，值越大土壤可容納水分越多，徑流量越小；SOL_K為土壤飽和水力傳導率，值越大土壤滲透性越強，徑流量越小。

表2 參數敏感性分析結果

2.2 岔巴溝流域SWAT模型模擬結果

在岔巴溝流域SWAT模型徑流模擬過程中，將2006—2007年為模型預熱期，以2008—2010年汛期為模型率定期，2011—2013年汛期為模型驗證期。以R2，ENS，PBIAS第3個指標對模型的適用性進行評價，表3是流域總出水口曹坪水文站汛期日尺度徑流率定期和驗證期的模擬評價結果，結果表明曹坪水文站的汛期流量模擬值與實測值吻合較好，率定期和驗證期的ENS系數分別為0.76，0.74，R2分別為0.78，0.75，模擬結果較好。PBIAS的計算結果表明模型在率定期將流量低估了11.3%，而在驗證期將流量高估了17.5%，這可能是由于驗證期特別是2013年的實測日徑流量變化很大，急起急落現象突出造成的。圖2為曹坪水文站率定期、驗證期的模擬結果，由圖可以看出SWAT模型可以較好的模擬流域汛期徑流的變化趨勢。

圖2 曹坪站汛期徑流實測值與模擬值比較

模型參數的不確定性采用P-factor和R-factor這2個指標來衡量，其中P-factor表示95 PPU區間包含實測數據的百分數，R-factor定義了95 PPU的頻帶平均寬度除以相應測量變量的標準偏差。理論上，P-factor的范圍是0～100%，R-factor的范圍是0～∞，P-factor接近1，R-factor接近0是完全接近于實測數據的模擬[27]。表3中模型率定期和驗證期P-factor分別為0.52，0.51，表明50%以上的實測數據在95 PPU置信區間范圍之內，而R-factor分別為0.32，0.25，表明模型率定的準確性可以接受，模型參數的不確定性較小。總體而言，以Green-Ampt下滲法驅動的SWAT模型在岔巴溝流域表現出較好的模擬精度，模型率定期和驗證期ENS系數均大于0.70，R2均達大于0.75，PBIAS介于±25%，參數不確定性較小，可以進一步分析流域不同土地覆被下的地表產流特征。

表3 模型模擬結果評價

2.3 岔巴溝流域地表產流特征分析

為進一步理解流域地表徑流的形成機制和變化規律，引入可以綜合降雨特征、前期土壤含水量和地表特征的變量—地表徑流系數作為評價指標[28]，在汛期徑流模擬的基礎上，基于率定好的模型，以次降雨事件內水文響應單元平均降雨量大于25 mm且平均地表產流在2.5 mm以上為條件篩選出四次降雨—徑流事件，結合640個水文響應單元累計得到2 560個不同降雨量—產流量的分析對象，用于分析不同土地覆被下的地表產流特征。

2.3.1 降雨強度對不同土地覆被產流特征的影響 降雨強度是黃土高原丘陵溝壑區產流產沙的主要影響因子[8]。根據最大1 h降雨強度等級和土地利用類型將2 560個分析對象的地表產流量和徑流系數進行分類統計得到不同最大雨強等級[29]下耕地、草地、林地的地表產流量和徑流系數(圖3)。從整體趨勢來看，隨著最大雨強的增加，3種土地覆被的地表產流量和徑流系數都逐漸上升，尤其當最大雨強從8.1～15.9 mm/h提高至16 mm/h以上時，3種土地覆被下地表產流量和徑流系數均顯著增大；不論在何種最大雨強區間下，3種土地覆被地表產流量以及徑流系數的關系均為耕地>草地>林地，且隨著最大雨強的增加三者間的差距也逐漸加大。

從不同土地覆被來看，地表產流隨不同土地利用類型的變化差異較明顯，最大雨強在2.5 mm/h以下時耕地、草地和林地均沒有產生地表徑流，耕地和草地在最大雨強達到2.6 mm/h后開始產生地表徑流，最大雨強在8 mm/h以上林地才開始出現產流；隨著最大雨強的增加耕地平均地表產流量從0 mm逐漸提升至0.91 mm，4.14 mm和12.83 mm，徑流系數在最大雨強在16 mm/h以上時達到0.28，是受雨強影響最顯著的土地利用類型，其次為草地，林地變化幅度最小，這是由于森林可以通過減小雨滴動能、攔截雨量以及改變土壤結構和地表結皮來削弱雨強對地表產流的影響[30]。綜合來看，流域地表徑流系數隨最大雨強的增加呈顯著上升趨勢，且最大雨強達到16 mm/h以上時地表徑流系數顯著增加。在同等級雨強下，耕地將降雨轉化為徑流的比例最大，其次為草地和林地。

圖3 不同降雨強度下土地覆被地表產流能力對比

2.3.2 前期土壤含水量對不同土地覆被產流特征的影響 前期土壤含水是指暴雨事件發生前流域表層土壤的含水量，在溫帶氣候下的許多研究表明前期的土壤水分條件是產生徑流的一個重要因素[31]。為了進一步分析不同土地利用類型下前期土壤含水對產流的影響，選取產流能力最強的降雨強度區間，四次事件累計出現498個最大雨強在16 mm/h以上的分析對象，其中耕地214個、草地197個、林地87個，得到其前期土壤含水量與徑流系數的擬合關系見圖4。從圖4A可以看出498個分析對象的初始土壤含水與地表徑流系數擬合R2為0.35，說明前期土壤含水量的大小可以揭示大雨強下部分地表徑流系數的變化，地表徑流系數隨前期土壤含水量的增加呈上升趨勢。

將分析對象按照土地覆被類型細分后擬合關系進一步提高(圖4B，C，D)，其中林地前期土壤含水量與地表徑流系數的擬合關系最好，R2為0.47，其次為耕地和草地；從擬合方程的斜率來看，3種土地覆被下地表徑流系數對前期土壤含水量的響應存在差異，其中耕地和草地的斜率分別為0.049，0.050，大于林地的0.044，說明耕地和草地地表徑流系數對前期土壤含水量的變化敏感程度更高，即隨著前期土壤含水量的增加地表徑流系數的變化幅度比林地更大。這可能是由于森林涵蓄土壤水分以及截留降水的作用減少了地表徑流的產生，造成其對前期土壤含水量變化敏感程度較低。

圖4 不同土地覆被前期土壤含水量與地表徑流系數擬合關系

3 討 論

以Green-Ampt下滲法為產流方法的SWAT模型在參數率定后達到了較好的模擬效果，但部分學者的研究結果認為降雨數據的時間分辨率對Green-Ampt下滲法的徑流模擬效果影響很大[32]，由于研究區內的雨量站均為汛期站，僅獲得了5—9月2 h步長的降雨摘錄數據，對降雨的強度輸入有所削弱，加上流域周邊氣象站點分布較少使得氣象因子的輸入精度有一定偏差，對模型產流有一定的影響。DEM數據的空間分辨率對坡度、河網及子流域劃分等建模過程有顯著影響，進而影響到流域徑流模擬精度，有研究表明20～150 m分辨率可以基本滿足流域模型模擬精度的要求[33]，但岔巴溝流域屬于典型的黃土丘陵溝壑區，地形破碎，溝壑縱橫，在陜北黃土高原的研究發現利用DEM提取的平均坡度隨數據分辨率的降低呈線性下降的態勢[34]，因此需要更高精度的DEM數據來識別研究區的坡度和溝壑密度。本研究由于未能獲取更精細的DEM數據，故使用30 m分辨率的DEM數據進行水文分析，對徑流模擬的精度可能會造成一定影響。在之后的研究中將利用無人機實地獲取分辨率更高的DEM數據來降低地形對徑流模擬的影響。

4 結 論

(1) 模型參數敏感性分析結果表明，徑流曲線數CN2，樹冠最大儲量CANMX、側向地下水流坡長SOSOIL、橫向流動時間LAT_TTIME、土壤有效含水量SOL_AWC和土壤飽和水力傳導度SOL_K等是流域日尺度徑流模擬較為敏感的參數；

(2) 以Green-Ampt下滲法作為產流方法的SWAT模型率定期和驗證期的ENS系數分別為0.76，0.74，R2分別為0.78，0.75，PBIAS介于±25%，參數不確定性較小，說明了其在黃土高原日徑流模擬的能力，可以用于分析降雨—徑流事件的地表產流規律；

(3) 同等雨強條件下，耕地將降雨轉化為地表徑流的比例最大，其次為草地和林地，且當降雨強度大于16 mm/h后地表徑流系數顯著增加；前期土壤含水量可以揭示大雨強下部分地表徑流系數的變化，其中耕地和草地地表徑流系數對前期土壤含水量的變化敏感程度更高。