董力軒, 常順利,*, 張毓?jié)?/p>

1 新疆大學資源與環(huán)境科學學院綠洲生態(tài)教育部重點實驗室,烏魯木齊 830046 2 新疆林科院森林生態(tài)研究所,烏魯木齊 830063

在森林生態(tài)系統(tǒng)中,降水量被分為林冠截留量、穿透降雨量和樹干莖流量三部分,其中林冠截留量在該系統(tǒng)的水文循環(huán)和水量平衡中占有重要地位,約占降水量的20%—40%[1—3]。因此,在應用SWAT模型[4—5]模擬森林生態(tài)水文過程時,林冠截留對地表徑流的影響不可忽視。

對林冠截留水文過程的模型化主要可以分為三類,即經(jīng)驗模型、半理論模型和理論模型[6—8]。雖然經(jīng)驗模型可以根據(jù)實測數(shù)據(jù)運用統(tǒng)計學方法直接推求林冠截留量和降雨量之間線性或非線性關系,但不能刻畫水文物理過程。因運算繁瑣,不易求解,如結(jié)合光線與林冠分配降雨規(guī)律形成的理論模型在實例中運行困難[6]。以Gash模型為代表的半理論模型相較于前兩者工作量小,參數(shù)易確定且運算求解方便,顯現(xiàn)了一定的優(yōu)勢[9]。Gash模型的適用性較強,能夠廣泛使用在不同氣候類型或林分類型的研究中,如海洋性氣候森林[9],熱帶與亞熱帶森林[10]、干旱半干旱區(qū)灌木林[11]、地中海氣候森林等[12]。在國內(nèi),Gash模型應用在大興安嶺地區(qū)落葉松林[13]、熱帶季雨林[14]、華北油松人工林[15]、貢嘎山暗針葉林[16]和黃土區(qū)人工刺槐林[17]等區(qū)域并獲得了比較合理的模擬效果。SWAT模型與Gash模型在計算平均蒸發(fā)速率時都選擇了Penman-Monteith公式,這為兩種模型的耦合提供了便利條件[4, 9]。因此,將林冠截留的物理過程引入SWAT模型的流域模擬中意義重大。

本文以天山北坡中段雪嶺云杉(Piceaschrenkiana)林為研究對象,結(jié)合氣象觀測、林下穿透雨實測數(shù)據(jù)和野外調(diào)查獲得模型所需參數(shù)進行建模。考慮到研究區(qū)潛在蒸發(fā)遠遠大于實際蒸發(fā)[18],林冠截留作用明顯,符合Gash模型對每次降雨事件發(fā)生時初始條件的要求,進行基于物理過程的林冠截留模擬可以提高模型精度,為水資源管理提供更可靠的依據(jù)。

1 研究區(qū)概況

研究區(qū)位于烏魯木齊縣境內(nèi)的烏魯木齊河流域(86°45′—87°56′E,43°00′—44°07′N),山區(qū)流域面積1070 km2,屬于中型流域,見圖1。河源位于天山烏魯木齊河源1號冰川,自南流向東北,出山口位于英雄橋。根據(jù)世界土壤數(shù)據(jù)庫(Harmonized World Soil Database,HWSD),土壤類型有簡育灰色土、鈣積黑鈣土、粘化栗鈣土、簡育栗鈣土、粘化鈣積土、永凍薄層土、松軟薄層土和冰川。研究區(qū)屬于溫帶大陸性氣候[18],春秋短暫,夏季晝夜溫差大,冬季漫長且多雪。日照時間長且日照輻射大,年平均氣溫為2—3℃,極端最高氣溫30℃(7月),極端最低氣溫-38℃(1月)。降水量少,年降水量400—600 mm,年內(nèi)降水多集中在6—8月；蒸發(fā)量大,年蒸發(fā)量980—1150 mm[19]；年均相對濕度64%—66%;干燥度1.2—1.6。該區(qū)域植被類型是以雪嶺云杉純林為主的溫帶針葉林,占流域面積的10%左右,研究區(qū)內(nèi)云杉分布廣泛且林冠截留作用明顯,流域特征符合研究需求,見圖1。

2 研究方法

2.1 模型的輸入數(shù)據(jù)

2.1.1氣象數(shù)據(jù)

選擇研究流域范圍內(nèi)的大西溝、英雄橋和小渠子三站1993—2012的日降水、日最高溫、日最低溫、平均風速等氣象數(shù)據(jù)構(gòu)建天氣發(fā)生器,整理后輸入氣象數(shù)據(jù)庫。

2.1.2數(shù)字高程數(shù)據(jù)

由地理空間數(shù)據(jù)云http://www.gscloud.cn/下載研究區(qū)DEM數(shù)字高程數(shù)據(jù),選擇30 m 分辨率數(shù)字高程數(shù)據(jù)并且轉(zhuǎn)為WGS_1984_UTM_Zone_45N地理投影(圖1)。

2.1.3土地利用數(shù)據(jù)

選擇對地物分類更好的哨兵數(shù)據(jù),經(jīng)最大似然分類法得到精度為10 m的土地利用數(shù)據(jù)(圖1)。

2.1.4土壤數(shù)據(jù)

本研究采用世界土壤數(shù)據(jù)庫,裁剪出研究區(qū),根據(jù)SAWT內(nèi)置土壤類型進行重分類和重采樣處理以滿足模型的精度要求,使用SPAW軟件計算建模所需參數(shù)建立土壤數(shù)據(jù)庫。

2.1.5林冠截留觀測區(qū)布設

依托天山森林生態(tài)系統(tǒng)定位站于2010年7—8月設立林冠截留觀測區(qū)(東經(jīng)87°27′,北緯43°25′)。該觀測區(qū)內(nèi)主要植被為雪嶺云杉,無灌木分布。在林間空地布設自記數(shù)字雨量儀測定總降水量和降雨強度,在林下設置接雨蓬鏈接集水桶收集穿透雨量后換算為降雨量深度(mm),具體數(shù)據(jù)見表1。

表1 林冠層對降雨的分配

圖1 烏魯木齊河流域示意圖Fig.1 Sketch map of Urumqi River basin

2.2 SWAT林冠截留模型

當使用SCS曲線法計算地表徑流時,不會對林冠截留進行具體物理過程的模擬,估計為產(chǎn)流前地表蓄水量和入滲量的20%。當使用Green & Ampt入滲方程計算地表徑流和入滲時,必須額外計算林冠截留量[20]。SWAT允許林冠截留的最大水量每天在變化,該方法根據(jù)葉面積指數(shù)函數(shù)來計算其值。公式如下：

(1)

式中,canday為當日林冠截留的最大水量(mm),canmx為樹冠充分發(fā)育時的林冠截留量(mm),該參數(shù)需要使用者手動輸入。LAI為該植被當日的葉面積指數(shù),LAImx為該植被的最大葉面積指數(shù)。

在SWAT模型中,得到給定日的降雨數(shù)據(jù),模型會通過降雨量的大小來判斷是否會有穿透雨形成。公式如下：

(2)

(3)

2.3 Gash模型

Gash模型描述的是一系列彼此分離的降雨事件,每個降雨事件都包含林冠加濕、林冠飽和、以及降雨停止后林冠干燥的過程,且假定每次降雨事件之間有足夠的時間讓林冠完全恢復到降雨前的干燥程度[9, 21]。模型將整個林冠在降雨過程中各個階段的截留損失相加得到總的林冠截留量。因為SWAT模型對于林冠截留的計算采用了單次降雨事件單獨計算的模式,并且根據(jù)降雨量的大小分為林冠截留量未飽和與林冠截留量飽和兩種情況,故將Gash模型拆解為分段函數(shù)以更好的與SWAT模型擬合。單一降雨事件發(fā)生時公式如下：

(4)

(5)

(6)

式中,S為林冠枝葉持水能力(mm)。

2.4 Gash模型參數(shù)的確定

圖2 降雨量與穿透雨量的回歸關系Fig.2 The relationship between throughfall and precipitation

自由穿透雨系數(shù)p,為林內(nèi)測得的平均可見天空率[23],經(jīng)過野外實地調(diào)查郁閉度后確定為0.2。

2.5 對SWAT模型子程序的修改

2.6 模型校準和驗證

模型的校準和驗證均在英雄橋出水口進行,本文通過決定系數(shù)R2、效率系數(shù)NSE、相對偏差系數(shù)PBIAS和均方根誤差RMSD來對模擬值和實測值的擬合度做出評價。決定系數(shù)R2的表達式[25]為：

(7)

Nash-Sutcliffe模型效率系數(shù)(NSE)的表達式[26]為：

(8)

偏差系數(shù)(PBIAS)的表達式為：

(9)

均方根誤差(RMSD)的表達式為：

(10)

標準差與皮爾遜相關系數(shù)的表達式為：

(11)

(12)

(13)

2.7 不確定性分析

SWAT模擬的流量不確定性主要是由于輸入數(shù)據(jù)的偏差,建模者的專業(yè)知識不足或改進模型的結(jié)構(gòu)不合理等原因?qū)е隆１疚闹惺褂昧朔治粩?shù)回歸方法,通過計算不同分位數(shù)的條件分布來進行非參數(shù)不確定性分析[27—29]。確定置信區(qū)間后,使用P因子和R因子來量化不確定性。P因子代表了95%的預測不確定性(95PPU),是落在95PPU范圍內(nèi)數(shù)據(jù)的數(shù)量與總數(shù)據(jù)數(shù)量的比值。評價校準質(zhì)量的另一個參數(shù)是R因子,表示95PPU范圍的平均寬度與觀測數(shù)據(jù)的標準差的比值。在理想情況下,P因子的取值在0—100%之間,R因子介于0—∞,當R<1.5時認為可接受,在P因子為1且R因子為0是模擬值與觀測值完全一致。

3 結(jié)果與分析

3.1 參數(shù)的敏感性分析

SWAT模型和SWAT-Gash模型均使用SWAT-CUP內(nèi)SUFI- 2算法進行模擬值的敏感性分析。首先根據(jù)類似研究區(qū)常用的25個參數(shù)[30—31]對SWAT模型進行1000次模擬,發(fā)現(xiàn)其中13個參數(shù)對SWAT模型模擬結(jié)果影響明顯(表2),再將其進行500次迭代得到最優(yōu)校準。同理,對SWAT-Gash模型進行參數(shù)的敏感性分析后選擇以下16個參數(shù)。經(jīng)分析發(fā)現(xiàn)兩種模型均對SCS徑流曲線系數(shù)、土壤飽和水力傳導系數(shù)、土壤可利用有效水、淺層含水層水位、融雪基溫和平均坡長等參數(shù)體現(xiàn)出更高的敏感性。

3.2 SWAT模型與SWAT-Gash模型在徑流模擬中的性能對比

考慮到氣象數(shù)據(jù)的年際變化,設置1993—1996年為模型的預熱期,以減少由于氣象數(shù)據(jù)而導致的潛在不確定性。1997—2004為校準期,2005—2012為驗證期。如圖3和圖4是SWAT模型與SWAT-Gash模型在校準期和驗證期的徑流模擬過程線,SWAT-Gash模型的模擬值優(yōu)于SWAT模型的模擬值,體現(xiàn)在SWAT-Gash模型更好的重現(xiàn)了高流量與低流量,SWAT模型會低估高流量。在校準期與驗證期,SWAT-Gash模型與SWAT模型在出山口位置的月平均徑流量分別為8.45—8.50 m3/s和6.55—6.86 m3/s之間;NSE值分別在0.63—0.85和0.58—0.82之間；R2值分別在0.65—0.86和0.59—0.83之間；PBIAS值分別在9.7%—11.5%和9.2%—17.1%之間(表3)；表明在該研究區(qū)SWAT-Gash模型的適用性高于SWAT模型。

校準期和驗證期兩種模型的均方根誤差(RMSD)、標準差(SD)和相關系數(shù)等均在泰勒圖中體現(xiàn)(圖5)。圖中各點與x軸上實測值點之間的徑向距離表示RMSD的大小,反映了各模型的模擬能力。SWAT-Gash模型得出的RMSD更低,在校準期和驗證期的值分別為3.22 m3/s和4.68 m3/s。使用SWAT模型時,RMSD值分別為3.49和7.80 m3/s,說明與SWAT-Gash模型相比SWAT模型在干旱半干旱區(qū)的徑流模擬效果欠佳。各點到原點之間的距離表示各組模擬值與觀測值的標準差。在校準期SWAT模型和SWAT-Gash模型的標準差分別為7.59和7.33,驗證期標準差分別為6.45和6.54。SWAT-Gash模型在驗證期標準差較大,可能是因為敏

表2 SWAT模型與SWAT-Gash模型敏感性參數(shù)

SWAT：水土評估工具 Soil and water assessment tool;CH_K2：主河道有效滲透系數(shù) Effective hydraulic conductivity in main channel;CN2：徑流曲線數(shù) Curve number of moisture condition Ⅱ;SOL_K：飽和滲透系數(shù) Hydraulic conductivity of the soil;SOL_AWC：土壤有效含水率 Available water capacity of soil layer;ALPHA_BF：基流α因子 Baseflow alpha factor;REVAPMN：發(fā)生再蒸發(fā)的淺層含水層水位閾值 Threshold depth of water in the shallow aquifer to “revap” to occur;GWQMN：發(fā)生回歸流的淺層含水層水位閾值 Threshold depth of water in the shallow aquifer required for return flow to occur;GW_DELAY：地下水時間延遲 Groundwater delay;SMTMP：融雪基溫 Snow melt base temperature;SFTMP：降雪基溫 Snowfall temperature;SURLAG：地表徑流滯后系數(shù) Surface runoff lag time;TLAPS：氣溫垂直變率 Temperature laps rate;SLSUBBSN：平均坡長 Average slope length; ALPHA_BNK：河岸調(diào)蓄基流α因子 Baseflow alpha factor for blank storage;CH_N2：主河道曼寧系數(shù) Manning′s "n" value for the main channel;GW_REVAP：地下水再蒸發(fā)系數(shù) Groundwater “revap” coefficient;TIMP：積雪溫度滯后因子 Snow pack temperature lag factor;PLAPS：降水垂直變率 Precipitation laps rate;CH_K1：支流有效滲透系數(shù) Effective hydraulic conductivity in tributary channel alluvium;CH_N1：支流曼寧系數(shù) Manning′s "n" value for the tributary channels;DEP_IMP：土壤不透水層深度 Depth to impervious layer in soil profile;EVPOT：壺穴政法系數(shù) pothole evaporation coefficient；EPCO：植物吸收補償因子 Plant uptake compensation factor

表3 模型徑流模擬評價指標

圖3 SWAT模型與SWAT-Gash模型校準期徑流模擬預測對比Fig.3 Comparison of runoff simulation prediction between SWAT model and SWAT gash model in calibration period

圖4 SWAT模型與SWAT-Gash模型驗證期徑流模擬預測對比Fig.4 Comparison of runoff simulation prediction between SWAT model and SWAT gash model in validation period

感性參數(shù)更多,但總體上該模型仍有更好的表現(xiàn)。相關系數(shù)表示出山口實測徑流數(shù)據(jù)和模型模擬值之間的關系,表明地表徑流在時間序列上的相似性。SWAT-Gash模型在校準期與驗證期的相關系數(shù)分別為0.93和0.81,高于SWAT模型的0.91和0.77。相較于SWAT模型,改進模型在校準期和驗證期模擬高流量的能力更加突出。

圖5 SWAT模型和SWAT-Gash模型在校準期和驗證期的泰勒圖Fig.5 Taylor diagram of calibration period and validation period for SWAT model and SWAT-Gash model

3.3 基于分位數(shù)回歸的不確定性分析

本文中定義了95%的置信區(qū)間β∈[0.025,0.975],使用兩種模型的驗證期徑流量與觀測流量進行基于分位數(shù)回歸的不確定性分析。如圖6所示,在驗證期SWAT模型的95%置信區(qū)間的平均厚度為13.50 m3/s,SWAT-Gash模型為12.86 m3/s,表明SWAT-Gash模型在本研究區(qū)的模擬表現(xiàn)更佳。SWAT-Gash模型的P因子為0.96,即96%的觀測值可以被該模型95%置信區(qū)間所包括,略高于SWAT模型(0.93)。此外,SWAT-Gash模型的R因子(1.19)低于SWAT模型(1.26),表明前者的不確定度范圍相對更窄。

圖6 SWAT模型和SWAT-Gash模型月徑流模擬的不確定性分布Fig.6 Uncertainty distribution of SWAT model and SWAT-Gash model in monthly runoff simulation

4 討論與結(jié)論

4.1 SWAT模型與SWAT-Gash模型在烏魯木齊河流域的應用

本文將Gash林冠截留模型與SWAT模型耦合后的模擬結(jié)果與原始SWAT模型的模擬結(jié)果進行對比分析,經(jīng)過校準和驗證后發(fā)現(xiàn)SWAT模型與SWAT-Gash模型在該流域的應用均取得了較好的效果,該結(jié)論與前人在烏魯木齊河流域的建模結(jié)果一致[30—32],但SWAT-Gash模型表現(xiàn)出更高的模擬精度。當使用SWAT模型在烏魯木齊河流域模擬時,校準期徑流重現(xiàn)表現(xiàn)較好,驗證期表現(xiàn)不佳,造成這種現(xiàn)象的原因可能是模型結(jié)構(gòu)復雜,敏感性參數(shù)對模擬結(jié)果的影響較大[33—35]。該流域內(nèi)坡面匯流過程、地下水與河道的補給排泄過程和融雪過程[36—37]是影響山區(qū)地表徑流的主要因素,僅通過調(diào)整描述水文過程的參數(shù)無法進一步提升模型在該區(qū)域的適用性。與SWAT模型相比,SWAT-Gash模型會高估該流域的低流量,Gash模型輸入數(shù)據(jù)的誤差也會影響SWAT-Gash模型地表徑流的預測。Gash模型計算時要求參數(shù)較多,雖然這些參數(shù)的確定方法經(jīng)過國內(nèi)外前人研究的驗證[9,12,17,22—23],可以得到較好的結(jié)果,但王艷萍等人[17]研究表明本文方法得出的參數(shù)模擬值偏高。Gash模型與SWAT模型耦合時,對林冠截留物理過程的描述更加詳細,結(jié)合雪嶺云杉林下穿透雨實測數(shù)據(jù)得出滿足Gash模型的參數(shù),使SWAT-Gash模型在該研究區(qū)的模擬表現(xiàn)更佳。

4.2 SWAT模型與SWAT-Gash模型的不確定性分析

本文采用了基于分位數(shù)的不確定性分析。與SWAT模型相比SWAT-Gash模型在基于分位數(shù)回歸的不確定性分析中的表現(xiàn)更佳,也表明Gash林冠截留模型對SWAT模型在烏魯木齊河流域的模擬改進效果明顯。通常在SWAT模型的建模中,將許多復雜水文過程用簡單的經(jīng)驗系數(shù)進行參數(shù)化,增加了模型的不確定性[27]。而本研究提出的SWAT-Gash模型中,將林冠截留過程和最大冠層截留量參數(shù)替換為基于物理過程的方程,從而提高了模型的適用性。SWAT-Gash模型的改進方法可以提高天山北坡中段中小型流域的水文模擬精度,為該區(qū)域水資源管理提供更可靠的依據(jù)。此外,本研究的降雨量分配數(shù)據(jù)代表性有所欠缺,還需要在不同地形條件下收集時間尺度更長的降雨分配數(shù)據(jù),使用更多樣的方法驗證Gash模型所需參數(shù),從而在該流域進行更精細化的建模。

4.3 結(jié)論

Gash模型參數(shù)易于獲取,能夠通過降雨量、穿透雨量和林冠截留量等數(shù)據(jù)推求,可以較好的在天山林區(qū)進行林冠截留過程的模擬。本文基于SWAT模型,將Gash模型方程代替經(jīng)驗公式提出新的SWAT-Gash模型。此方法可以提高烏魯木齊河流域水文模擬的精度。同時使用基于分位數(shù)的不確定分析發(fā)現(xiàn)改進后模型與出山口處實測地表徑流數(shù)據(jù)的相似性更高,適用性高于原始SWAT模型。