楊 光, 魯克新, 李 鵬,, 柳海亮, 劉桂華, 曹峰華, 馬天文, 王得軍, 文妙霞, 孫景梅
(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室, 西安 710048; 2.西安理工大學(xué) 旱區(qū)生態(tài)水文與災(zāi)害防治國(guó)家林業(yè)局重點(diǎn)實(shí)驗(yàn)室, 西安 710048)
SWAT模型是美國(guó)農(nóng)業(yè)部農(nóng)業(yè)研究中心于20世紀(jì)90年代研發(fā)的以SWRRB為框架,同時(shí)融入了CREAMS、GLEAMS、EPIC和ROTO等模型優(yōu)點(diǎn)的流域尺度長(zhǎng)時(shí)段分布式水文模型[8]。國(guó)內(nèi)外學(xué)者多使用SWAT模型進(jìn)行流域水沙研究。Zhang等[9]有關(guān)土地利用變化對(duì)于渾河流域產(chǎn)水產(chǎn)沙的影響的模擬結(jié)果表明,林地全年產(chǎn)沙量減少,水分入滲增加;草地和林地產(chǎn)流產(chǎn)沙對(duì)土地利用變化的響應(yīng)相似,但前者的水土保持效果弱于后者。邢立文等[10]研究結(jié)果表明,氣候變化對(duì)產(chǎn)沙的影響比產(chǎn)流大,產(chǎn)流、產(chǎn)沙與降雨呈正相關(guān)關(guān)系而與氣溫呈負(fù)相關(guān)關(guān)系,降雨量的變化對(duì)產(chǎn)流和產(chǎn)沙的影響比氣溫大。劉世梁等[11]有關(guān)瀾滄江中游小流域水土流失變化與NDVI時(shí)空變化的研究結(jié)果表明,瀾滄江中游小流域的徑流產(chǎn)沙量在空間尺度上隨著流域植被NDVI的增大而減小,在時(shí)間尺度上與流域植被NDVI間的關(guān)系不顯著。朱楠等[12]研究了羅玉溝土地利用結(jié)構(gòu)對(duì)流域水沙的影響,結(jié)果表明林地的減水、減沙效果最好,坡耕地的產(chǎn)流、產(chǎn)沙量最大;當(dāng)林地分布在流域上游時(shí),到達(dá)流域中游出口的水沙量最小。然而,以往研究學(xué)者主要利用SWAT模型研究氣候變化和土地利用變化對(duì)流域水沙的影響等,而對(duì)于黃河流域上游徑流侵蝕功率的空間分布規(guī)律的研究還不夠深入。
本研究以寧夏回族自治區(qū)水土流失嚴(yán)重的清水河流域?yàn)檠芯繉?duì)象,基于SWAT模型開(kāi)展清水河流域徑流模擬,并使用程圣東[4]提出的能更好表示流域年徑流過(guò)程的年徑流侵蝕功率公式計(jì)算各子流域出口斷面的年徑流侵蝕功率,進(jìn)而進(jìn)行年徑流侵蝕功率的空間分布特征和尺度效應(yīng)研究,以期從徑流侵蝕功率角度解釋流域水力侵蝕的空間分布情況,為后續(xù)在清水河流域不同區(qū)域開(kāi)展針對(duì)性的水土保持措施布設(shè)以及生態(tài)治理工作提供理論支撐。
清水河是黃河寧夏段最大的支流,地理位置為東經(jīng)105°09′—116°43′,北緯35°49′—37°30′,發(fā)源于六盤(pán)山東麓固原市,向北流經(jīng)原州區(qū)、海原縣、同心縣、中寧縣等縣(區(qū)),在中寧縣的泉眼山西側(cè)注入黃河,全長(zhǎng)320 km,流域面積14 481 km2,河道平均縱比降1.49‰。流域地貌以黃土丘陵溝壑區(qū)為主,地勢(shì)南高北低,河源海拔2 489 m,入黃口海拔1 190 m,相對(duì)高差1 299 m。
清水河流域多年平均降水量349 mm,年降水量空間分布不均,由南到北遞減,流域上、下游的年降水量分別為600,200 mm。流域內(nèi)年平均氣溫由南向北遞減,北部中衛(wèi)市年平均氣溫8℃以上,中部同心縣年平均氣溫7~8℃,南部固原市年平均氣溫5~6℃。清水河流域日照多、濕度小、風(fēng)大,水面蒸發(fā)強(qiáng)烈,多年平均水面蒸發(fā)量1 272 mm,年水面蒸發(fā)量介于900~1 300 mm[13]。清水河流域多年平均入黃泥沙量約占寧夏黃河流域入黃泥沙總量的49%[14]。
SWAT模型的輸入數(shù)據(jù)包括流域數(shù)字高程模型(DEM)數(shù)據(jù)、土地利用數(shù)據(jù)、土壤類(lèi)型數(shù)據(jù)及逐日降水量等氣象數(shù)據(jù)。根據(jù)土壤類(lèi)型數(shù)據(jù)使用SWAT模型的SPAW工具計(jì)算得到土壤下滲率、電導(dǎo)度、可持水量等土壤參數(shù)數(shù)據(jù),進(jìn)而制備SWAT模型運(yùn)行所需的土壤參數(shù)數(shù)據(jù)庫(kù)。根據(jù)逐日降水量氣象數(shù)據(jù),使用SWAT Weather工具計(jì)算各氣象站點(diǎn)的各月月降水量的均值標(biāo)準(zhǔn)差、最高和最低氣溫均值和標(biāo)準(zhǔn)差、月均濕度、月均太陽(yáng)輻射以及月均風(fēng)速數(shù)值,制備天氣發(fā)生器,為后續(xù)流域徑流模擬提供參考。
本研究使用的清水河流域空間數(shù)據(jù)包括中科院地理空間數(shù)據(jù)云30 m分辨率的DEM數(shù)據(jù)、中國(guó)科學(xué)院遙感所解譯的1980年流域土地利用數(shù)據(jù)以及世界土壤數(shù)據(jù)庫(kù)(Harmonized World Soil Database,HWSD)中的1∶100萬(wàn)比例尺流域土壤類(lèi)型數(shù)據(jù)。
本研究使用的水文氣象數(shù)據(jù)包括:1976—1986年《黃河流域水文年鑒》摘錄歷年的清水河流域22個(gè)雨量站的逐日降水量以及清水河干流韓府灣水文站的逐日平均流量;來(lái)自國(guó)家氣象科學(xué)數(shù)據(jù)中心的中衛(wèi)、中寧、同心、海原、固原和西吉6個(gè)氣象站點(diǎn)的逐日降水量、氣溫、風(fēng)速、太陽(yáng)輻射以及相對(duì)濕度等氣象數(shù)據(jù)。清水河流域水文站、氣象站和雨量站空間分布圖見(jiàn)圖1。
圖1 清水河流域DEM數(shù)據(jù)、站點(diǎn)分布
本研究使用ARCSWAT 2012版本對(duì)清水河流域進(jìn)行徑流模擬。在進(jìn)行流域劃分時(shí),本研究把集水區(qū)面積閾值設(shè)置為10 000 hm2并以清水河入黃口作為流域出口,將整個(gè)清水河流域劃分為93個(gè)子流域,詳見(jiàn)圖1B。在建立水文響應(yīng)單元(HRU)時(shí),需要輸入土地利用數(shù)據(jù)、土壤類(lèi)型數(shù)據(jù)和坡度數(shù)據(jù)。
土地利用類(lèi)型對(duì)流域產(chǎn)匯流具有重要影響,將預(yù)先準(zhǔn)備好的1980年清水河流域的土地利用數(shù)據(jù)重分類(lèi)為耕地、林地、草地、水域、城鎮(zhèn)居民用地以及未利用地6類(lèi)土地利用類(lèi)型。統(tǒng)計(jì)結(jié)果表明,清水河流域最主要的土地利用類(lèi)型為草地,面積占流域面積的56%;其次為耕地,面積占流域面積的37%,而其他4種土地利用類(lèi)型面積合計(jì)占流域面積的比例不足10%。清水河流域的耕地主要沿河道兩岸分布,而流域北部分布有大范圍的草地以及未利用地。清水河流域1980年土地利用類(lèi)型空間分布圖見(jiàn)圖2。
圖2 清水河流域1980年土地利用類(lèi)型空間分布
輸入土壤數(shù)據(jù),從HSWD的土壤數(shù)據(jù)庫(kù)中提取清水河流域的土壤數(shù)據(jù),重新分類(lèi)為23種不同的土壤類(lèi)型;將流域內(nèi)的坡度劃分為5級(jí)(0~15%,15%~25%,25%~35%,35%~45%,45%~99.99%)。
在HRU定義時(shí),將土地利用數(shù)據(jù)、土壤類(lèi)型數(shù)據(jù)以及流域坡度的閾值都設(shè)置為5%,共劃分出2 058個(gè)HRU。通過(guò)閾值的設(shè)定,可以減少面積過(guò)小的土地利用、土壤類(lèi)型以及坡度對(duì)于模型模擬結(jié)果精度的影響[15]。將整理好的氣象數(shù)據(jù)如逐日的降水、氣溫、風(fēng)速、太陽(yáng)輻射以及相對(duì)濕度數(shù)據(jù)輸入SWAT模型,隨后將其寫(xiě)入數(shù)據(jù)庫(kù)后即可運(yùn)行SWAT模型。
為了確保SWAT模型模擬結(jié)果的準(zhǔn)確性,本研究以清水河干流韓府灣水文站1976年、1977年作為預(yù)熱期,以1978—1982年歷年各月實(shí)測(cè)月平均流量數(shù)據(jù)進(jìn)行模型參數(shù)率定,以1983—1986年歷年各月實(shí)測(cè)月平均流量數(shù)據(jù)進(jìn)行模型驗(yàn)證。
影響SWAT模型模擬結(jié)果的參數(shù)眾多,為了減少模型運(yùn)行的不確定性和提高模型的運(yùn)行效率,需要對(duì)模型的相關(guān)參數(shù)進(jìn)行敏感性分析。
本研究使用LH-OAT算法對(duì)SWAT模型的參數(shù)進(jìn)行敏感性分析,以p和t檢驗(yàn)值作為參數(shù)敏感程度的評(píng)判標(biāo)準(zhǔn),p值越接近于0且t值越大,則說(shuō)明參數(shù)越敏感[16]。本研究選擇對(duì)模型模擬結(jié)果敏感性較高的參數(shù)參與后續(xù)的模型參數(shù)率定和驗(yàn)證。
SWAT模型運(yùn)行成功后,需要驗(yàn)證模型模擬結(jié)果的可靠性。本研究選取清水河干流的韓府灣水文站的實(shí)測(cè)徑流資料對(duì)模型模擬結(jié)果進(jìn)行驗(yàn)證。選取參數(shù)敏感性分析中的敏感性較高的參數(shù)參與模型率定和驗(yàn)證。使用SUFI-2算法對(duì)率定期參數(shù)進(jìn)行率定,得到最優(yōu)參數(shù)值,進(jìn)而將最優(yōu)參數(shù)值帶入驗(yàn)證期進(jìn)行參數(shù)驗(yàn)證。
本研究使用決定系數(shù)R2以及納什效率系數(shù)NS對(duì)清水河SWAT模型的適用性進(jìn)行評(píng)估。R2表示兩個(gè)變量間的相關(guān)程度,值越接近于1說(shuō)明兩個(gè)變量間的相關(guān)關(guān)系越好;NS一般用于評(píng)價(jià)觀(guān)測(cè)值與模擬值的擬合程度,值越接近于1說(shuō)明擬合情況越好,當(dāng)NS=1時(shí),觀(guān)測(cè)值和模擬值完全相等。一般研究[17-21]認(rèn)為,當(dāng)R2和NS均大于0.5時(shí),模型模擬結(jié)果可信。如果模型率定期和驗(yàn)證期的模擬結(jié)果均可信,則認(rèn)為最優(yōu)參數(shù)值適用于研究流域,模型模擬結(jié)果能較好地表示研究流域的實(shí)際徑流情況。
程圣東[4]將徑流侵蝕功率由場(chǎng)次尺度推廣到年尺度的,考慮年徑流深為Hy,年內(nèi)最大月平均流量為Qm,得到年尺度徑流侵蝕功率計(jì)算公式如下:
(1)
(2)
(3)
使用清水河干流韓府灣水文站1978—1986年的實(shí)測(cè)月平均流量數(shù)據(jù)對(duì)清水河流域SWAT模型進(jìn)行參數(shù)率定和驗(yàn)證。以1979—1982年的月平均流量為率定期數(shù)據(jù),采用LH-OAT方法進(jìn)行模型參數(shù)敏感性分析,選擇18個(gè)敏感性強(qiáng)的參數(shù)進(jìn)行率定,結(jié)果見(jiàn)表1。
從表1中可以看出,對(duì)清水河流域徑流模擬結(jié)果影響最大的參數(shù)為決定流域徑流量大小的SCS徑流曲線(xiàn)數(shù)CN2,CN2越大則地表徑流量就越大;其次是主要影響下滲情況的表層土壤容重SOL_BD以及主要表示河岸地下徑流對(duì)于河道流量的補(bǔ)給快慢的主河道調(diào)蓄系數(shù)ALPHA_BNK。另外,其余參數(shù)如降雪氣溫SFTMP、淺層地下水徑流系數(shù)GWQMN、基流消退系數(shù)ALPHA_BF、主河道曼寧系數(shù)CH_N2都對(duì)流域徑流模擬結(jié)果的影響較為敏感。
表1 清水河流域SWAT模型參數(shù)敏感性排名
將收集到的韓府灣水文站實(shí)測(cè)月平均流量數(shù)據(jù)劃分為預(yù)熱期(1976年和1977年)、率定期(1978—1982年)和驗(yàn)證期(1983—1986年)3個(gè)階段。在SWAT-CUP軟件中,使用率定期數(shù)據(jù)對(duì)上文選取的18個(gè)參數(shù)對(duì)進(jìn)行參數(shù)率定,將參數(shù)進(jìn)行優(yōu)化迭代,直至決定系數(shù)R2>0.6且納什效率系數(shù)NS>0.5,則可認(rèn)為模型的模擬結(jié)果可信。將率定期得到的最優(yōu)參數(shù)值回代到驗(yàn)證期后再次運(yùn)行SWAT模型,并根據(jù)驗(yàn)證期的月流量模擬值與實(shí)測(cè)值再次計(jì)算決定系數(shù)R2和納什效率系數(shù)NS。若驗(yàn)證期的決定系數(shù)R2>0.6且納什效率系數(shù)NS>0.5,則認(rèn)為SWAT模型的模擬結(jié)果可信,確定的模型最優(yōu)參數(shù)值可用于后續(xù)的月流量模擬計(jì)算。
韓府灣水文站率定期(1978—1982年)和驗(yàn)證期(1983—1986年)月平均流量實(shí)測(cè)值與模擬值對(duì)比結(jié)果見(jiàn)圖3。
圖3 1978-1986年韓府灣站月平均流量過(guò)程模擬值與實(shí)測(cè)值對(duì)比
經(jīng)計(jì)算,模型率定期的決定系數(shù)R2為0.76、納什效率系數(shù)NS為0.74;驗(yàn)證期的決定系數(shù)R2為0.88,納什效率系數(shù)NS為0.85。因此認(rèn)為清水河流域月流量過(guò)程的模擬結(jié)果可信,SWAT模型在清水河流域具有良好的適用性。
利用SWAT模型可以得到1978—1986年清水河流域入黃口斷面以及流域93個(gè)子流域出口斷面的月流量過(guò)程數(shù)據(jù)以及歷年年徑流深等數(shù)據(jù)。
清水河流域SWAT模型輸出文件中的Sub表中的子流域產(chǎn)水量(WYLD)字段包含1978—1986年93個(gè)子流域的逐月月產(chǎn)水量數(shù)據(jù),以此可以推求出各子流域的多年平均徑流深H;經(jīng)統(tǒng)計(jì)分析,從Sub表中查得各子流域歷年的最大月徑流深hm。將徑流侵蝕功率公式進(jìn)一步變換得到公式(4)—(5)。
hmA′=QmΔt
(4)
(5)
式中:hm為年內(nèi)最大月徑流深(m);A′為子流域面積(km2);Δt為時(shí)段長(zhǎng)(s),此處平均按每月30.4 d計(jì)算各月總秒數(shù),Δt=2.627×106s。
由公式(5)得到1978—1986年清水河各子流域歷年年內(nèi)最大月平均流量模數(shù)Q′m,最后由公式(1)計(jì)算得到相應(yīng)年份各子流域的年徑流侵蝕功率。利用1978—1986年清水河各子流域歷年的年徑流侵蝕功率結(jié)果求得多年平均徑流侵蝕功率,并基于A(yíng)rcGIS平臺(tái)制作清水河流域多年平均徑流侵蝕功率空間分布圖(圖4)。
從圖4可以看出,清水河流域多年平均徑流侵蝕功率呈現(xiàn)“支流大,干流??;東部大,西部小”的空間分布規(guī)律,即各支流區(qū)域的多年平均徑流侵蝕功率大于清水河干流,且流域西部地區(qū)的多年平均徑流侵蝕功率明顯小于東部地區(qū)。
圖4 1976-1986年清水河流域多年平均徑流侵蝕功率空間分布 圖5 1976-1986年清水河流域各子流域出口斷面多年平均徑流侵蝕功率空間分布
為了探究流域徑流匯聚效應(yīng)與流域徑流侵蝕功率的關(guān)系,本節(jié)進(jìn)一步研究了子流域年徑流侵蝕功率與子流域出口斷面集水面積間的相關(guān)關(guān)系。
SWAT模型輸出的Rch表包含各子流域出口斷面的集水面積和月流量(FLOW-OUT)。子流域出口斷面的集水面積指子流域出口斷面以上的集水區(qū)總面積。子流域出口斷面流量(FLOW-OUT)指子流域出口斷面以上匯聚到子流域出口處的流量。以各子流域出口斷面的集水面積A′和多年平均月平均流量(FLOW-OUT),依據(jù)公式(1)—(3)計(jì)算各子流域出口斷面多年平均徑流侵蝕功率。
由圖5可知,在各級(jí)支流的徑流匯聚過(guò)程中,子流域出口斷面的多年平均徑流侵蝕功率都有減少的趨勢(shì),尤其以流域上游地區(qū)較為明顯,但此趨勢(shì)在清水河下游集水面積較大的子流域間不太明顯。因此,本研究將出口斷面集水面積小于4 000 km2的子流域的集水面積與子流域出口斷面多年平均徑流侵蝕功率進(jìn)行擬合分析,試圖找到子流域出口斷面集水面積與子流域出口斷面多年平均徑流侵蝕功率間的相關(guān)關(guān)系,二者關(guān)系散點(diǎn)圖見(jiàn)圖6。由圖6可知,子流域出口斷面多年平均徑流侵蝕功率隨著子流域出口斷面集水面積的增大而迅速減小,直到出口斷面集水面積達(dá)到某一數(shù)值時(shí)子流域出口斷面多年平均徑流侵蝕功率對(duì)其敏感性降低以至于趨于穩(wěn)定。
圖6 子流域出口斷面集水面積與多年平均徑流侵蝕功率擬合分析
經(jīng)回歸分析,當(dāng)子流域出口斷面集水面積小于4 000 km2時(shí),清水河子流域出口斷面多年平均徑流侵蝕功率E與子流域出口斷面集水面積Ac間呈如下顯著的冪函數(shù)關(guān)系
(6)
經(jīng)統(tǒng)計(jì)分析,當(dāng)子流域出口斷面集水面積大于4 000 km2時(shí),出口斷面多年平均徑流侵蝕功率數(shù)值穩(wěn)定在1.56×10-5m4/(s·km2)左右。分析其主要原因可能是清水河流域右岸支流的匯水量過(guò)大,導(dǎo)致在下游地區(qū)干流與其他支流間的多年平均徑流侵蝕功率空間分布規(guī)律不明顯。
為確定子流域出口斷面集水面積的臨界值,將公式(6)進(jìn)行求導(dǎo)得到如下導(dǎo)數(shù)方程:
(7)
|E′|<0.001說(shuō)明子流域出口斷面多年平均徑流侵蝕功率變化速率減小,其值趨于穩(wěn)定[6]。當(dāng)出口斷面集水面積Ac大于84.85 km2時(shí),|E′|<0.001,子流域出口斷面多年平均徑流侵蝕功率不再隨著子流域集水面積的增大而明顯減少,因此影響清水河流域多年平均徑流侵蝕功率大小的子流域出口斷面集水面積空間閾值為84.85 km2。
依據(jù)魯克新等[3]提出的流域次暴雨輸沙模型研究結(jié)果,流域次暴雨徑流侵蝕功率與流域輸沙模數(shù)呈冪函數(shù)關(guān)系,徑流侵蝕功率越大則流域輸沙模數(shù)越大。因此,在進(jìn)行清水河流域生態(tài)治理時(shí),優(yōu)先選擇處于流域上中游區(qū)且出口斷面集水面積小于84.85 km2的小流域,可取得良好的水土流失治理效果。
通過(guò)本研究以及龔珺夫等[5-6]在無(wú)定河及延河的徑流侵蝕功率空間分布的研究發(fā)現(xiàn),不同流域的年徑流侵蝕功率存在區(qū)域性差異的同時(shí),也有共同的特點(diǎn)。從空間角度來(lái)看,流域年徑流侵蝕功率的空間分布情況與流域下墊面條件、氣候情況相關(guān),不同流域會(huì)具有不同的年徑流侵蝕功率空間分布特點(diǎn)。從流域水系的角度來(lái)看,不同流域的多年平均徑流侵蝕功率都存在“支流大,干流小”的空間分布規(guī)律。
為了更好研究子流域出口斷面徑流侵蝕功率與子流域出口斷面集水面積間的相關(guān)關(guān)系,本研究相較于龔珺夫[5-6]的研究選取了程圣東[4]提出的年徑流侵蝕功率,使用年徑流深作為徑流侵蝕能量的參數(shù),以期更好地表達(dá)流域的徑流侵蝕功率情況。研究結(jié)果表明,清水河子流域多年平均徑流侵蝕功率與子流域出口斷面集水面積間呈顯著的冪函數(shù)關(guān)系,與龔珺夫等[5-6]的研究結(jié)果相似。不同流域的徑流侵蝕功率的子流域集水面積閾值并不相同。對(duì)于清水河流域而言,小于流域集水面積閾值的子流域均位為各支流的河源處,即圖5中徑流侵蝕功率較大的區(qū)域,這對(duì)于清水河流域生態(tài)治理具有一定的實(shí)際指導(dǎo)意義。
本研究基于SWAT模型模擬研究了清水河流域的年徑流侵蝕功率的空間分布規(guī)律及其空間尺度效應(yīng),但并未闡明清水河流域的多年平均徑流侵蝕功率存在“支流大,干流小”空間分布規(guī)律的機(jī)理。另外,本研究沒(méi)有進(jìn)行清水河流域產(chǎn)沙輸沙量模擬,因此無(wú)法將模擬得到的清水河流域徑流侵蝕功率用于流域輸沙模數(shù)計(jì)算。以上兩個(gè)方面的不足還有待今后深入開(kāi)展相關(guān)研究,以期為流域侵蝕產(chǎn)沙預(yù)報(bào)與流域生態(tài)治理提供理論指導(dǎo)。
(1) 本文構(gòu)建的清水河流域SWAT模型的率定期(1978—1982年)決定系數(shù)R2為0.76,納什效率系數(shù)NS為0.74;驗(yàn)證期(1983—1986年)決定系數(shù)R2為0.88,納什效率系數(shù)NS為0.85,因此SWAT模型在清水河流域徑流模擬中具有良好的適用性,模擬的徑流結(jié)果可以較好地表示流域的實(shí)際徑流情況。
(2) 清水河流域的多年平均年徑流侵蝕功率具有“支流大,干流?。粬|部大,西部小”的空間分布規(guī)律。
(3) 清水河子流域多年平均徑流侵蝕功率隨著出口斷面集水面積的增大而迅速減小,而當(dāng)子流域出口斷面集水面積達(dá)到臨界閾值時(shí)出口斷面徑流侵蝕功率對(duì)其敏感性降低并趨于穩(wěn)定。當(dāng)子流域的出口斷面集水面積小于4 000 km2時(shí),子流域多年平均徑流侵蝕功率與出口斷面集水面積間呈顯著的冪函數(shù)關(guān)系。清水河上中游流域多年平均徑流侵蝕功率的流域出口斷面集水面積臨界閾值約為84.85 km2。優(yōu)先選擇處于清水河上中游區(qū)且出口斷面集水面積小于84.85 km2的小流域進(jìn)行水土保持措施布設(shè)及生態(tài)治理,可取得良好的水土流失治理效果。