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

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

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

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

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

1 研究區(qū)概況

研究區(qū)位于烏魯木齊縣境內(nèi)的烏魯木齊河流域(86°45′—87°56′E,43°00′—44°07′N),山區(qū)流域面積1070 km2,屬于中型流域,見圖1。河源位于天山烏魯木齊河源1號(hào)冰川,自南流向東北,出山口位于英雄橋。根據(jù)世界土壤數(shù)據(jù)庫(Harmonized World Soil Database,HWSD),土壤類型有簡育灰色土、鈣積黑鈣土、粘化栗鈣土、簡育栗鈣土、粘化鈣積土、永凍薄層土、松軟薄層土和冰川。研究區(qū)屬于溫帶大陸性氣候[18],春秋短暫,夏季晝夜溫差大,冬季漫長且多雪。日照時(shí)間長且日照輻射大,年平均氣溫為2—3℃,極端最高氣溫30℃(7月),極端最低氣溫-38℃(1月)。降水量少,年降水量400—600 mm,年內(nèi)降水多集中在6—8月；蒸發(fā)量大,年蒸發(fā)量980—1150 mm[19]；年均相對(duì)濕度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的日降水、日最高溫、日最低溫、平均風(fēng)速等氣象數(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ù)

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

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

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

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

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

表1 林冠層對(duì)降雨的分配

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

2.2 SWAT林冠截留模型

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

(1)

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

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

(2)

(3)

2.3 Gash模型

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

(4)

(5)

(6)

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

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

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

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

2.5 對(duì)SWAT模型子程序的修改

2.6 模型校準(zhǔn)和驗(yàn)證

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

(7)

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

(8)

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

(9)

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

(10)

標(biāo)準(zhǔn)差與皮爾遜相關(guān)系數(shù)的表達(dá)式為：

(11)

(12)

(13)

2.7 不確定性分析

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

3 結(jié)果與分析

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

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

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

考慮到氣象數(shù)據(jù)的年際變化,設(shè)置1993—1996年為模型的預(yù)熱期,以減少由于氣象數(shù)據(jù)而導(dǎo)致的潛在不確定性。1997—2004為校準(zhǔn)期,2005—2012為驗(yàn)證期。如圖3和圖4是SWAT模型與SWAT-Gash模型在校準(zhǔn)期和驗(yàn)證期的徑流模擬過程線,SWAT-Gash模型的模擬值優(yōu)于SWAT模型的模擬值,體現(xiàn)在SWAT-Gash模型更好的重現(xiàn)了高流量與低流量,SWAT模型會(huì)低估高流量。在校準(zhǔn)期與驗(yàn)證期,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模型。

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

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

SWAT：水土評(píng)估工具 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：地下水時(shí)間延遲 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：植物吸收補(bǔ)償因子 Plant uptake compensation factor

表3 模型徑流模擬評(píng)價(jià)指標(biāo)

圖3 SWAT模型與SWAT-Gash模型校準(zhǔn)期徑流模擬預(yù)測對(duì)比Fig.3 Comparison of runoff simulation prediction between SWAT model and SWAT gash model in calibration period

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

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

圖5 SWAT模型和SWAT-Gash模型在校準(zhǔn)期和驗(yàn)證期的泰勒?qǐng)DFig.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],使用兩種模型的驗(yàn)證期徑流量與觀測流量進(jìn)行基于分位數(shù)回歸的不確定性分析。如圖6所示,在驗(yàn)證期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),表明前者的不確定度范圍相對(duì)更窄。

圖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模型在烏魯木齊河流域的應(yīng)用

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

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

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

4.3 結(jié)論

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