基于SWAT模型的沂蒙山水土流失治理區(qū)地表水與地下水轉(zhuǎn)化規(guī)律研究

譚秀翠，莊會波，季 妤

(1.山東農(nóng)業(yè)大學水利土木工程學院，山東 泰安 271018；2.山東省水文局，濟南 250002)

地表水與地下水的水量轉(zhuǎn)化是水文循環(huán)的重要過程，地表水與地下水的相互作用及轉(zhuǎn)化關系是水文地質(zhì)等領域研究的熱點和難點，分析和掌握其規(guī)律特征，對進行水資源分析評價、物質(zhì)遷移和能量變化的相關研究具有重要意義[1]。目前，對地表水與地下水轉(zhuǎn)化的研究方法有：水文測量法[2]、基流分割法[3]、水量平衡法[4]、水文地球化學法[5,6]、溫度示蹤法[7,8]、數(shù)值模型法[9,10]等。

植被作為陸地生態(tài)環(huán)境中的重要組成部分，對水分具有吸收、駐留、緩釋等作用，是水文循環(huán)的重要環(huán)節(jié)。氣候變化和土地利用/覆被變化引起的各種時間尺度下流域下墊面植被特性的改變將對流域水文循環(huán)產(chǎn)生顯著的影響[11]。Loch[12]對比分析澳大利亞昆士蘭州南部植被覆蓋度分別為0、23%、37%、47%以及100%區(qū)域的產(chǎn)流和地表入滲，結(jié)果顯示植被覆蓋度增加將大大增加地表入滲，從而減小地表產(chǎn)流。石峰[13]利用LCM模型對兩種不同情景下次降雨徑流進行模擬，表明林草植被增加都明顯減小了流域產(chǎn)流量。史曉亮等[14]研究發(fā)現(xiàn)諾敏河流域植被覆蓋變化對徑流影響較小。成向榮等[15]構(gòu)建6種不同植被覆蓋情景，研究該植被覆蓋類型變化對徑流的影響。目前，關于植被覆蓋的水文效應主要集中在地表產(chǎn)流、入滲，對于地表水與地下水轉(zhuǎn)化規(guī)律的研究較少。

山東省沂蒙山區(qū)地處山東省中南部、淮河流域北部，是淮河流域重要的水源區(qū)和生態(tài)屏障，也是淮河流域水土流失最嚴重的地區(qū)之一。2006年，水利部公告了42個國家級水土流失重點防治區(qū)，山東省的沂蒙山區(qū)被列為重點治理區(qū)，其空間分布位置見圖1。經(jīng)過多年治理，2016年，沂蒙山治理區(qū)的監(jiān)測區(qū)域中蒙陰縣和沂水縣，水土流失面積1 577.61 km2，占土地總面積的39.34%，與2013年調(diào)查結(jié)果相比，水土流失面積減少11.29 km2(中國水土保持公報，2016)。本文以山東省沂蒙山區(qū)蒙陰縣岸堤水庫控制流域為研究對象，分析氣候與植被覆蓋時空變化對地表水與地下水轉(zhuǎn)化的影響規(guī)律，研究成果可為沂蒙山區(qū)水資源的合理開發(fā)利用、水土流失治理提供理論依據(jù)和技術支撐。

圖1 沂蒙山水土流失治理區(qū)位置分布圖Fig.1 The location map of Yimeng mongtain soil and water loss control area

1 研究區(qū)域

蒙陰縣岸堤水庫控制流域地處沂蒙山區(qū)腹地，蒙山北麓，東汶河上游。流域地貌受構(gòu)造、巖性、氣候、河流等內(nèi)外營力作用的控制和影響，整個地形南北高，中間低，由北西向東南傾斜。流域地處中緯度地帶，屬暖溫帶半濕潤大陸性季風型氣候，具有四季分明，光照充足，雨量充沛，無霜期長等特點，多年平均降水量781.8 mm，降水量多集中在6-8月，占全年的65%，多年平均溫度13.4 ℃(蒙陰氣象站，1959-2018年數(shù)據(jù))。汛期徑流量占年徑流量的82.3%，易形成春旱夏洪，約有2/3徑流量形成洪水流失。蒙陰縣降水量和徑流量總的空間分布特征為：南多北少，東多西少，山區(qū)大，平地少，從東南向西北遞減[16]。

蒙陰縣城東部東汶河北部，屬層間巖溶裂隙水，水位埋深2 m左右。縣城南部東汶河南側(cè)，淺層屬松散堆積物孔隙水、沖積粗砂礫石層潛水，水位埋深3 m左右。山南莊村東部屬層間巖溶裂隙水。東汶河入云蒙湖口西南部屬裂隙巖溶水。淺層地下水主要靠大氣降水補給，其水位深度隨季節(jié)變化幅度較大。

流域內(nèi)分布有3個水文測站(1984-2016年逐日數(shù)據(jù))，19個雨量站(1984-2016年逐日數(shù)據(jù))，1個氣象站(1959-2018年逐日數(shù)據(jù))，各站點空間分布見圖2。

圖2 研究流域分布圖Fig.2 The location map of research watershed

歸一化植被指數(shù)NDVI(normalized difference vegetation index)對植被的生長勢和生長量非常敏感，因此NDVI不僅可以反映植被、土地利用變化特征，也可以反映氣候?qū)χ脖挥绊懀ǔ１挥糜谕恋乩?覆蓋變化、生態(tài)環(huán)境監(jiān)測以及植被動態(tài)監(jiān)測等方面[17]。文中采用的NDVI數(shù)據(jù)來源于中國年度、月度植被指數(shù)(NDVI)空間分布數(shù)據(jù)集[18]，該數(shù)據(jù)集是基于連續(xù)時間序列的SPOT/VEGETATION NDVI衛(wèi)星遙感數(shù)據(jù)，采用最大值合成法生成的植被指數(shù)數(shù)據(jù)集。2000年流域平均NDVI為0.68，至2016年流域平均NDVI為0.76，NDVI空間分布圖見圖3。

圖3 流域NDVI空間分布Fig.3 Spatial distribution of NDVI in the watershed

2 研究方法

2.1 Mann-Kendall突變檢驗

Mann-Kendall突變檢驗是一種非參數(shù)統(tǒng)計檢驗方法[19]，其優(yōu)點是不需要樣本遵從一定的分布，也不受少數(shù)異常值的干擾，更適用于類型變量和順序變量，計算比較簡便。

對于有n個樣本的時間序列x，構(gòu)造一個秩序列：

(1)

可見，秩序列Sk是第i時刻數(shù)值大于j時刻數(shù)值個數(shù)的累計數(shù)。在時間序列隨機獨立的假定下，定義統(tǒng)計量：

(2)

式中：UF1=0，E(sk)和var(sk)是sk的均值和方差，在x1,x2,…,xn相互獨立，且有相同連續(xù)分布時，它們可由下式算出：

(3)

(4)

UFi為標準正態(tài)分布，它是按時間序列x順序x1,x2,…,xn計算出的統(tǒng)計量序列，給定顯著性水平α，若|UFi|>Uα，則表明序列存在明顯的趨勢變化。按時間序列的逆序xn,xn-1,…,x1，再重復上述過程，同時使UBk=-UFk(k=n,n-1,…,1)，UB1=0。

2.2 NDVI空間變化趨勢分析

基于一元線性回歸的方法，構(gòu)建研究區(qū)域每個柵格單元NDVI值與時間的回歸方程，采用回歸方程的斜率表示植被指數(shù)的變化趨勢[20]，計算公式如下所示：

(5)

式中：k為柵格單元回歸方程的斜率；pi為第i年NDVI均值；n為時間序列，取值為17(2000-2016年)。

k>0，表示該柵格單元NDVI值呈增加趨勢，植被覆蓋情況改善；k<0，表明該柵格單元NDVI值呈減小趨勢，植被覆蓋情況退化。

2.3 SWAT數(shù)值模擬

SWAT模型是由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心Jeff Arnold博士研發(fā)的基于物理過程的流域尺度模型，計算包括降水、地表徑流、蒸散發(fā)、壤中流、滲透、地下水回流和河道運移損失等水量的水平衡關系，水量平衡方程如下[21]：

(6)

式中：SWt為土壤最終含水量，mm；SW0為第i天土壤初始含水量，mm；t為時間步長，d；Rday為第i天降水量，mm；Qsurf為第i天的地表徑流量，mm；Ea為第i天的蒸散發(fā)量，mm；Wseep為第i天從土壤剖面進入包氣帶的水量，mm；Qgw為第i天回歸流的水量，mm。SWAT模型產(chǎn)流計算流程，如圖4所示。

圖4 SWAT模型計算流程圖Fig.4 The flow chart of SWAT model

SWAT模型需要的主要輸入數(shù)據(jù)包括：DEM數(shù)據(jù)、土地利用數(shù)據(jù)、土壤類型數(shù)據(jù)、氣象數(shù)據(jù)、水文數(shù)據(jù)及土壤屬性數(shù)據(jù)，各數(shù)據(jù)來源如表1所示。

表1 SWAT 模型所需數(shù)據(jù)類型Tab.1 Data types for SWAT model

借助ArcGIS10.2，建立SWAT2012模型數(shù)據(jù)庫，基于研究流域DEM數(shù)據(jù)，劃分38個子流域。將土地利用類型重分類，對土地利用類型進行編碼賦值以匹配土地利用類型數(shù)據(jù)庫及作物數(shù)據(jù)庫，土地利用類型分布見圖5。

圖5 土地利用類型分布圖Fig.5 Distribution of land use types

土壤數(shù)據(jù)包括：土壤類型分布圖(見圖6)與土壤屬性數(shù)據(jù)，土壤基礎屬性數(shù)據(jù)主要來源世界土壤數(shù)據(jù)庫(HWSD)，根據(jù)土壤的基礎參數(shù)，如黏土(Clay)、砂土(Sand)、有機質(zhì)(Organic Matter)等含量運用SPAW軟件計算土壤有效持水量、濕密度、飽和導水率。

圖6 土壤類型分布圖Fig.6 Distribution of soil types

土壤侵蝕力因子(USLE_K)采用Williams提出的方程計算：

KUSLE=fcsandfcl-siforgcfhisand

(7)

式中：fcsand、fcl-si、fhisand分別為粗糙沙土地質(zhì)、黏壤土、高沙質(zhì)土壤侵蝕因子；forgc為土壤有機質(zhì)因子，可通過土壤各粒徑組成和有機碳含量計算得到。一般情況下，USLE_K值在0.02～0.75之間。

影響土壤產(chǎn)流能力的屬性是指那些影響土壤在完全濕潤并且不凍的條件下的最小下滲率屬性[22]，根據(jù)最小下滲率確定土壤水文分組，見表2。

表2 土壤水文分組Tab.2 Soil hydrologic groups

氣象資料，采用蒙陰站1984-2016年逐日溫度、風速、相對濕度、蒸發(fā)等氣象數(shù)據(jù)及19個雨量站1984-2016年逐日降水數(shù)據(jù)。

水文資料，采用蒙陰、水明崖、岸堤水庫3個水文1984-2016年的逐日徑流數(shù)據(jù)。水庫數(shù)據(jù)，采用岸堤水庫1984-2016年逐月出流數(shù)據(jù)。

對于模型適用性的驗證，采用相關系數(shù)(R2)和納什效率系數(shù)(ENS)來判斷，但目前并沒有公認的統(tǒng)一標準，綜合近年來的研究成果，Motovilov等[23]認為模擬徑流量時ENS>0.75結(jié)果較好，在0.36～0.75之間結(jié)果較為滿意；郝芳華等[22]、Braemort[24]認為R2>0.60、ENS>0.50時結(jié)果較為滿意。

3 結(jié)果分析

3.1 NDVI特征分析

2006年，沂蒙山被劃為國家級水土流失重點治理區(qū)，經(jīng)過多年的治理，植被覆蓋條件改善，由圖7(a)可以看出，2000-2016年流域NDVI呈增大趨勢，2006年前后，流域平均NDVI由0.7增大到0.74。由Mann-Kendall法突變分析結(jié)果可以看出，見圖7(b)，除2002年外，NDVI統(tǒng)計量UFk均大于0，表明NDVI呈增大趨勢，其中，2007-2009年，UFk明顯大于顯著性水平線，NDVI顯著增大。在0.05置信度線內(nèi)，UFk與UBk存在一個交點，該交點為NDVI的突變點。

圖7 流域NDVI特征分析Fig.7 Characteristics analysis of NDVI in the watershed

由2000-2016年年度NDVI資料，計算每個柵格NDVI的一元回歸方程的斜率，計算結(jié)果如圖8所示。k取值范圍在-0.019～0.011之間，k值大于零的面積占流域面積84%，表明流域內(nèi)大部分面積NDVI都呈現(xiàn)增加趨勢，植被覆蓋條件改善。k值小于零的區(qū)域主要集中在蒙陰縣城周邊，主要受城鎮(zhèn)發(fā)展，土地利用方式變化的影響。

圖8 2000-2016年NDVI變化趨勢空間分布圖Fig.8 Spatial distribution of NDVI trend from 2000 to 2016

3.2 SWAT模擬結(jié)果分析

采用SWAT-CUP對參數(shù)進行調(diào)試、優(yōu)化，1984-1989年作為模型預熱期，1990-2016年作為校準期和驗證期。根據(jù)蒙陰、水明崖和岸堤水庫3個水文站的實測流量數(shù)據(jù)驗證模擬結(jié)果，如圖9所示，3個水文站實測流量與模擬流量的過程線，變化趨勢基本一致，峰值相對應，模擬效果良好，3個水文站徑流量模擬結(jié)果與實測結(jié)果的R2及ENS均>0.5，表明調(diào)試后的模型適用于研究流域。

圖9 實測流量與模擬流量過程線Fig.9 Process line of measured runoff and simulated runoff

SWAT計算結(jié)果顯示，1990-2016年，流域平均降水量611.96 mm，實際蒸散發(fā)量383.51 mm，徑流總量225.96 mm，地表徑流量126.07 mm，基流量72.07 mm，壤中流25.02 mm，潛水補給量119.52 mm。根據(jù)子流域輸出文件(output.sub)，繪制各子流域水文要素空間分布圖見圖10。

圖10 水文要素空間分布圖Fig.10 Spatial distribution of hydrological elements

SWAT模型的預熱期為1984-1989年，數(shù)據(jù)分析采用1990-2016年。由圖11可以看出，1990-2016年流域水文要素的時間變化，從水循環(huán)三要素角度考慮，降水、蒸散發(fā)、徑流呈周期性變化，其中徑流的變化趨勢及峰值位置與降水基本一致，兩者相關度較高(R2=0.92)，且線性趨勢變化平穩(wěn)，蒸散發(fā)量變幅較小，與徑流量相關度較低，如圖11(a)。徑流量由3種水源匯聚而成，包括地表水、土壤水、地下水，其分布見圖10(b)，其中占比最大的是地表水，占55%，呈線性減小趨勢，在降水量變化平穩(wěn)的情況下，主要受下墊面因素影響；其次是地下水，占27%，基流量年際變化比較明顯，呈線性增加趨勢，多年平均基流系數(shù)為0.27；土壤水對徑流量的貢獻最小，占18%。潛水的補給主要來源于降水和地表水，由圖10(c)可以看出，潛水補給呈下降趨勢，地表水對其貢獻率較小，僅占5%左右，主要補給水源為降水。

圖11 水文要素時間分布圖Fig.11 Temporal distribution of hydrological elements

為分析植被覆蓋變化對地表水與地下水轉(zhuǎn)換關系的影響，以2006年為分界點，采用SWAT數(shù)值模擬結(jié)果，分析各水文要素的變化量，統(tǒng)計結(jié)果見表3，可以看出，2006年前后，地表徑流、基流變幅分別為-22%、30%，除受降水影響外，主要受下墊面因素的影響。2006年前后，流域的NDVI值增大，表明植被覆蓋條件改善，促進了地表水流入滲，因而地表徑流減少，地表水補給地下水量減少，基流量增大。綜上，在流域降水量減少、植被覆蓋度增大的情況下，地表水向地下水轉(zhuǎn)化量減少，地下水向地表水轉(zhuǎn)化量增多。

表3 2006年前后水文要素統(tǒng)計Tab.3 Statistics of hydrological elements around 2006

4 結(jié) 論

本文對沂蒙山國家水土流失治理區(qū)典型小流域的NDVI進行特征分析，采用SWAT模型開展氣候與植被覆蓋影響下的地表水與地下水的轉(zhuǎn)化規(guī)律研究，得到如下結(jié)論。

(1)2000-2016年，岸堤水庫控制流域NDVI呈增大趨勢，Mann-Kendall突變分析表明，2007-2009年，NDVI顯著增大。流域內(nèi)80%以上面積NDVI都呈現(xiàn)增加趨勢。

(2)通過驗證，SWAT模型適用于研究流域。1990-2016年，流域平均降水量611.96 mm，實際蒸散發(fā)量383.51 mm，徑流總量225.96 mm，地表徑流量126.07 mm，基流量72.07 mm，壤中流25.02 mm，潛水補給量119.52 mm。

(3)流域徑流與降水變化趨勢基本一致，兩者相關度較高(R2=0.92)。徑流量中，地表徑流占55%，基流量占27%，壤中流占18%。潛水的補給中，地表水對其貢獻率較小，僅占5%左右，主要補給水源為降水。

(4)2006年前后，岸堤水庫流域降水量減少、植被覆蓋度增大，地表水向地下水轉(zhuǎn)化量減少，地下水向地表水轉(zhuǎn)化量增多。

□