WRF-Hydro模型與新安江模型在陳河流域的應(yīng)用對比*

孫明坤 ，李致家，劉志雨,，侯愛中，霍文博，溫婭惠,3，孔祥意，戴金旺,4，梁世強,5

(1：河海大學(xué)水文水資源學(xué)院，南京210098) (2：水利部信息中心(水利部水文水資源監(jiān)測預(yù)報中心)，北京100053) (3：黃河水利委員會河南水文水資源局，鄭州450003) (4：寧波市鎮(zhèn)海區(qū)建筑(交通)工程安全質(zhì)量管理站，寧波315202) (5：北京水文總站，北京100089)

暴雨與洪水災(zāi)害嚴(yán)重威脅人民的生命財產(chǎn)安全，如何系統(tǒng)地研究兩者間的聯(lián)系并及時準(zhǔn)確做出水文氣象預(yù)報一直是研究的熱點難點問題[1]. 一般地，暴雨事件可通過大氣數(shù)值模式模擬[2-3]，如WRF模式[3-5]；而洪水一般通過水文模型預(yù)報[6-7]，如新安江模型[8]和TOPMODEL模型[9]等. 事實上，大氣模式和水文模型尺度不盡相同，前者一般尺度較大，后者網(wǎng)格較精細(xì). 近些年來，出現(xiàn)了將兩者相耦合的建模方法(即陸氣耦合模型)，該法可顯著延長洪水預(yù)報的預(yù)見期，為暴雨洪水的模擬和預(yù)報提供有力工具[10-12]. 例如，高冰等[11]將WRF模式與分布式水文模型GBHM進(jìn)行耦合，較好地模擬三峽庫區(qū)入庫洪水過程. 但多數(shù)研究均采用單向的陸氣耦合方式，即僅將大氣模式的輸出作為水文模型的輸入，而未考慮后者對前者的反饋機制. 雙向耦合的水文氣象模型可以較好地考慮兩個過程間的交互作用，更合理地模擬水文氣象狀態(tài)，提供更可靠的暴雨洪水預(yù)報.

WRF-Hydro模型[13]是有物理基礎(chǔ)的、多尺度的、多參數(shù)化方案的、可與大氣模式進(jìn)行單向或雙向耦合的、新一代的分布式水文氣象預(yù)報系統(tǒng)，能定量研究流域內(nèi)大氣-陸面水熱交換過程[14]. WRF-Hydro模型以較大尺度的陸面模式[2]為基礎(chǔ)，加入基于精細(xì)網(wǎng)格的匯流模塊，以解決某些陸面模式中未考慮側(cè)向流的問題. 同時可以更好地模擬土壤水的分布情況，為大氣模式提供更合理的下邊界條件. 此外，該模型在不同的水文氣象過程中能提供多個參數(shù)化方案以供選擇，如河道演算模塊提供馬斯京根法、馬斯京根-康吉法與擴散波演算法[13]. 近4年來，國外學(xué)者對該系統(tǒng)進(jìn)行了深入研究. 例如Yucel等[15]探究了該模型的參數(shù)率定和移植方法，發(fā)現(xiàn)在無資料流域利用移植的參數(shù)和WRF同化數(shù)據(jù)所得的模擬結(jié)果誤差較小. Senatore等[16]發(fā)現(xiàn)WRF/WRF-Hydro雙向耦合系統(tǒng)的模擬結(jié)果略優(yōu)于僅使用WRF模式的結(jié)果. Naabil等[17]認(rèn)為WRF-Hydro模型在降水模擬、洪水預(yù)報、水文水資源管理等各個領(lǐng)域都有應(yīng)用潛力. 而在實際應(yīng)用之前，從水文學(xué)角度研究WRF-Hydro模型在我國匯流時間較短的中小流域中的表現(xiàn)及適用性很有必要.

為了評價該模型在中小流域的表現(xiàn)和適用性，本文利用在我國濕潤半濕潤地區(qū)應(yīng)用較好的新安江模型進(jìn)行對比分析. 新安江模型是經(jīng)典的概念性半分布式水文模型，也經(jīng)常作為工具來評價其他模型在水文方面的表現(xiàn)[11-12,18]. 例如謝帆等[19]在淮河的息縣流域?qū)Π敕植渴剿哪Ｐ蚑OPMODEL模型和新安江模型進(jìn)行對比，發(fā)現(xiàn)兩者均可達(dá)到精度要求. 李致家等[20]在新安江的呈村流域分別利用分布式水文模型TOPKAPI模型[21]和新安江模型進(jìn)行洪水模擬，發(fā)現(xiàn)兩者精度相差不大，均可用于水文過程的模擬與評價.

本文分別利用陸氣耦合模型WRF-Hydro模型和新安江模型對陳河流域4年共16場洪水進(jìn)行模擬分析，討論WRF-Hydro模型在研究流域的表現(xiàn)及適用性，以及兩個模型在結(jié)構(gòu)、參數(shù)、輸入輸出和模擬結(jié)果方面的差異性. 考慮到尺度問題，引入土層厚度乘子ZSOILFAC對WRF-Hydro模型的次表面層的厚度進(jìn)行調(diào)節(jié)，以合理地提高模擬精度. 同時利用新安江模型反推的包氣帶厚度與調(diào)節(jié)后的次表面層厚度進(jìn)行對比分析，進(jìn)一步論證ZSOILFAC的必要性與合理性. 此外，探究WRF-Hydro模型在陳河流域中合適的計算步長，并分析可能的原因，為之后利用單向或雙向耦合的WRF-Hydro模型在中小流域進(jìn)行水文模擬和實時洪水預(yù)報提供參考.

1 研究流域與方案

1.1 陳河流域概況

選取陜西省渭河右岸一級支流黑河的陳河水文站以上集水區(qū)域為研究流域(簡稱陳河流域). 陳河流域地處陜西省中南部，流域面積約為1380 km2，流域內(nèi)黑河干流河長約為105 km. 該流域恰好位于“秦嶺-淮河”的北部，為典型的半濕潤流域，流域徑流深100～500 mm. 多年平均降水量為800 mm，且多集中在夏、秋兩季，易形成洪水. 夏季降雨歷時較短，雨強較大；秋季降雨歷時較長，雨強較小. 該流域地處秦嶺褶皺帶，地形以山地為主，其中秦嶺主峰太白山位于流域西北部. 海拔高程約為600～3500 m，落差較大. 流域平均坡度為26.73°，其中坡度為10°～50°的面積約占總面積的95.3%，而20°～40°的面積約占68.2%. 根據(jù)美國地質(zhì)勘探局(USGS)資料[22]，陳河流域的表層土壤均為壤土，主要土地利用類型有：稀樹草原(60.9%)、草地(12.3%)、混合林(8.9%)和落葉闊葉林(7.9%). 流域內(nèi)有8個雨量站和1個水文站(圖1).

圖1 陳河流域水系及水文站的分布Fig.1 Geographic location of rain/flow stations and main rivers in Chenhe Basin

1.2 資料與研究方案

流域內(nèi)可用資料為2003-2012年的日流量資料和以小時為尺度的部分洪水資料，其中洪水事件共計24場. 而2003、2006、2010、2011年內(nèi)有記錄的洪水共16場，占總場次的66.7%，其余年份僅有1～2場. 考慮計算成本和效率問題，選擇這4年6-10月為研究時段，詳細(xì)信息見表1. 注意到本研究中所有洪水事件均以1 h時段資料為基礎(chǔ). 根據(jù)WRF-Hydro模型的相關(guān)研究[15-16]，將2003和2006兩年共6場洪水作為率定期，2010和2011兩年共10場洪水作為檢驗期. WRF-Hydro模型需要預(yù)熱以獲得與實際較為相符的流域初始狀態(tài). 對于每一年，預(yù)熱期的跨度為汛期中有記錄的第一場洪水的前一個月左右. 利用全球再分析資料

表1 預(yù)熱期、模擬期和研究洪水事件信息

FNL(Final Operational Global Analysis)驅(qū)動WRF模式，所獲得的初始狀態(tài)場視為WRF-Hydro模型預(yù)熱期的初始狀態(tài)場. 其中FNL資料來源于https://rda.ucar.edu/. 相對于假設(shè)的模型狀態(tài)場，這種方法可以獲得相對準(zhǔn)確的初始場，適當(dāng)?shù)乜s短預(yù)熱期長度，提高計算效率. WRF-Hydro模型所需的土地利用、土壤類型和反射率等均源于WRF模式前處理系統(tǒng)WPS的數(shù)據(jù)庫http://www2.mmm.ucar.edu/wrf/. 此外該模型需要的高分辨率數(shù)字高程資料源自http://www.gscloud.cn/.

由于WRF-Hydro模型是基于陸面模式的水文模型系統(tǒng)，需要大量的氣象強迫數(shù)據(jù)，包括近地面的氣溫、氣壓、比濕、風(fēng)速、向下長波輻射、向下短波輻射和降水率(雨強). 一方面，陳河流域目前沒有可以獲得以上全部資料的觀測站點，這為驅(qū)動模型帶來很大困難. 而全球陸面數(shù)據(jù)同化系統(tǒng)(GLDAS)提供全球范圍內(nèi)較大尺度的氣象再分析資料，為資料匱乏地區(qū)提供寶貴的強迫數(shù)據(jù). 本文利用分辨率為3 h，0.25°的GLDAS資料，時間上經(jīng)線性插值，空間上經(jīng)雙線性插值[23]，得到陳河流域1 h、1 km的強迫數(shù)據(jù)場. 另一方面，在洪水模擬中，降雨的時空分布往往起很大的作用. 因此，利用反距離權(quán)重法(以下簡稱IDW法)對雨量站觀測數(shù)據(jù)進(jìn)行空間插值，得到1 h、1 km的降雨時空分布場. 將處理后的強迫數(shù)據(jù)場中的非降雨場與利用IDW法插值后的降雨場結(jié)合，得到最終的模型強迫數(shù)據(jù).

新安江模型的主要輸入是站點的降雨數(shù)據(jù)和蒸發(fā)皿觀測數(shù)據(jù). 蒸發(fā)皿蒸發(fā)數(shù)據(jù)源于流域外黑峪口站(位于流流域出口東北方向約9 km處)，經(jīng)換算得到流域蒸發(fā)能力. 對于雨量數(shù)據(jù)，將整個流域按照每個自然子流域有且僅有一個觀測站的原則進(jìn)行劃分. 在每個自然子流域內(nèi)，對IDW法插值后的降雨場求面平均值作為模型輸入，以保證兩個模型降雨的一致性.

1.3 評價指標(biāo)

本文選取徑流深誤差(AErd)、洪峰相對誤差(REfp)、峰現(xiàn)時間誤差(AEpt)及納什系數(shù)(NSE)作為評價模型表現(xiàn)的主要指標(biāo)，具體計算方法見文獻(xiàn)[24]. 在模擬和實測過程之間，AErd反映兩者洪量的差異，REfp與AEpt則綜合評價了模型在洪峰方面的表現(xiàn)，NSE表示模擬與實測過程在線型上的擬合程度. 此外計算合格洪水?dāng)?shù)與總洪水?dāng)?shù)之比(即合格率)以整體衡量模型的表現(xiàn). 當(dāng)某場洪水的徑流深誤差小于徑流深許可誤差(PAErd)時，則徑流深預(yù)報合格；當(dāng)洪峰相對誤差在20%以內(nèi)時，則洪峰預(yù)報合格[24]. 其中PAErd是實測徑流深的函數(shù)：當(dāng)實測徑流深小于15 mm時，PAErd取3 mm；當(dāng)實測徑流深大于15 mm小于100 mm時，PAErd取實測徑流深的20%；當(dāng)實測徑流深大于100 mm時，PAErd取20 mm[25]. 此外，在陳河流域內(nèi)，若相鄰時刻流量增量大于10 m3/s，且后續(xù)時刻至洪水極大值處流量持續(xù)增加，則將第一個時刻定義為洪水起漲時刻. 模擬過程與實測過程的洪水起漲時刻的差值定義為洪水起漲點誤差，單位為h.

2 模型與參數(shù)

2.1 WRF-Hydro模型

2.1.1 模型簡介 WRF-Hydro模型是美國國家氣象研究中心(NCAR)開發(fā)的陸氣耦合模型系統(tǒng)，主要包括不同尺度的陸面模式和匯流模塊[13]. 陸面模式的相關(guān)過程均在尺度較大的網(wǎng)格(如1 km，簡稱陸面模式網(wǎng)格)進(jìn)行，如下滲等；水文匯流的相關(guān)過程均在尺度較小的網(wǎng)格(如100 m，簡稱匯流網(wǎng)格)進(jìn)行，如地表徑流等. 在模型運行過程中，變量或參數(shù)需要在兩個不同尺度間相互轉(zhuǎn)換，具體方法見文獻(xiàn)[13]和[23]. 模型的本質(zhì)是在陸面模式的基礎(chǔ)上添加匯流模塊以更精細(xì)地考慮水文過程. 本文中WRF-Hydro模型各個部分的選取方案見表2.

表2 WRF-Hydro模型系統(tǒng)的主要物理過程配置方案*

*指該過程目前無其他方案可供選擇.

Noah陸面模式主要模擬大氣下層、植被冠層和淺層土層之間的水量和能量的交互過程[2,27]. 冠層截留[26]、超滲地面徑流、土壤含水量和蒸散發(fā)均在陸面模式網(wǎng)格中計算. Noah模式的淺層土層被稱為次表面層，主要計算水量能量的分布情況. 其總厚度為2.0 m，共4層，每層默認(rèn)厚度由上至下分別為0.1、0.3、0.6、和1.0 m. 當(dāng)雨強大于最大下滲率(式1)時，地表出現(xiàn)超滲地面徑流.

(1)

式中，Imax為最大下滲率(m/s)，Pd為經(jīng)冠層截留后的雨強(m/s)，Dx為次表面層總?cè)彼?m/s)，Ks為飽和水力傳導(dǎo)度(m/s)，Kref=2.0×10-6m/s，kdtref為可調(diào)參數(shù)，δi為與時段轉(zhuǎn)換有關(guān)的系數(shù)，δi與Dx的計算詳見文獻(xiàn)[26]. 次表面層的濕度和溫度分別由Richard方程和熱擴散方程計算. Noah模式的總蒸發(fā)分為地表直接蒸發(fā)、冠層截留量蒸發(fā)和植物根系散發(fā)3部分，蒸發(fā)能力由Penman公式計算[27].

匯流網(wǎng)格主要計算次表面徑流(即壤中流)、坡面流、地下徑流和河網(wǎng)匯流等. 利用達(dá)西定律計算側(cè)向飽和壤中流：

(2)

式中，Qi為飽和壤中流的流量(m3/s)，αx和αy分別為x和y方向的水面坡度，d為土層厚度(m)，w為網(wǎng)格邊長(m)，h為水面深度(m)，n為可調(diào)參數(shù)，默認(rèn)為1. 在計算步長Δt內(nèi)，根據(jù)水量平衡計算當(dāng)前網(wǎng)格的水深變化Δh. 模型認(rèn)為每個網(wǎng)格的地表積水由3部分構(gòu)成：超滲地面徑流、飽和地面徑流及與其他網(wǎng)格的水量交換. 當(dāng)且僅當(dāng)網(wǎng)格水深hs超過參數(shù)地表積水深hret時，出現(xiàn)可自由流動的地表徑流. 地表徑流由淺水波的擴散波方程計算，詳見文獻(xiàn)[6,13]. 其中，單寬流量與糙率的關(guān)系由曼寧公式計算：

(3)

式中，q為單寬流量(m2/s)，nov為坡面的曼寧糙率，Sf為摩阻坡度. 使用概念性水箱模型描述流域?qū)Υ伪砻鎸右韵禄?即地下徑流)的調(diào)蓄作用，該水箱的水深并不代表真實含水層中的水深[13]. 水箱的入流量為次表面層最下層向下的滲透量，出流量采用經(jīng)驗蓄泄關(guān)系計算，詳見文獻(xiàn)[13]. 河道形狀均采用等腰梯形，其底寬、邊坡系數(shù)、糙率為河道等級的函數(shù)，在相關(guān)參數(shù)文件中有預(yù)定值，也可在隨后過程中率定. 河道演算采用一維擴散波方程：

(4)

式中，A為過水?dāng)嗝婷娣e(m2)，Q為河道流量(m3/s)，nch為河道曼寧糙率，H為河道水面高程(m)，R為水力半徑(m)，sgn(·)為符號函數(shù).

2.1.2 參數(shù)率定 WRF-Hydro模型可調(diào)參數(shù)較多，主要分為陸面網(wǎng)格參數(shù)(如土層的田間持水量SMCREF)和匯流網(wǎng)格參數(shù)(地表糙率乘子參數(shù)). 模型中多利用乘子參數(shù)在整個空間范圍內(nèi)對某些屬性值(如地表糙率)進(jìn)行合理地等比縮放. 考慮到次表面層是針對大尺度的陸面或天氣模式而設(shè)計的，其厚度可能不適用于中小流域的洪水模擬. 因此，本文借鑒其他乘子參數(shù)，引入土層厚度乘子ZSOILFAC. 利用該參數(shù)在合適范圍內(nèi)等比例地對整個次表面層厚度進(jìn)行調(diào)節(jié). 相比于文獻(xiàn)[16]將次表面層厚度直接定為0.05、0.2、0.45和0.8 m，經(jīng)過ZSOILFAC修改的厚度更加合理. 綜上，本文選取5個主要參數(shù)進(jìn)行率定，分別是土層厚度乘子ZSOILFAC、滲透參數(shù)REFKDT(kdtref)、地表截流深乘子RETDEPRTFAC(hret的乘子參數(shù))、地表糙率乘子OVROUGHRTFAC(nov的乘子參數(shù))和河道曼寧糙率乘子MANNFAC(nch的乘子參數(shù)). 其中，前3個參數(shù)主要控制徑流分配和洪水水量，后兩個參數(shù)主要控制洪峰和洪水過程線的形狀.

調(diào)參的目的是使率定期內(nèi)盡可能多的模擬洪水與實測過程充分接近，使AErd和REfp較小，同時使NSE較大. 考慮到計算成本，對主要參數(shù)使用手動逐步逼近法調(diào)參，以避免消耗大量時間[15]. 手動逐步逼近法是指對某一參數(shù)在合理的范圍內(nèi)以適當(dāng)?shù)牟介L分別驅(qū)動模型，結(jié)合指標(biāo)定量評價模擬結(jié)果，從而確定最佳的參數(shù). WRF-Hydro模型在陳河流域的參數(shù)率定結(jié)果見表3.

表3 WRF-Hydro模型系統(tǒng)的主要參數(shù)及率定數(shù)值

2.2 新安江模型

三水源新安江模型是經(jīng)典的概念性半分布式水文模型，在我國濕潤和半濕潤地區(qū)有較好的應(yīng)用. 模型結(jié)構(gòu)包含三層蒸散發(fā)、蓄滿產(chǎn)流、三水源劃分和流域匯流4個模塊[6,8]. 鑒于新安江模型在陳河流域有較好的應(yīng)用，本文使用該模型對16場洪水進(jìn)行模擬，并與WRF-Hydro模型的結(jié)果比較，評價后者的表現(xiàn)與適用性. 利用自然子流域分割法將陳河流域分為9塊，每塊有且只有1個雨量/水文站，以此考慮降水的空間差異性. 新安江模型的主要參數(shù)率定結(jié)果見表4.

表4 新安江模型的主要參數(shù)及率定數(shù)值

2.3 模型結(jié)構(gòu)對比

WRF-Hydro模型是基于物理的分布式水文模型系統(tǒng)，可定量考慮下墊面屬性的空間差異，并與WRF模式有很好的相容性. 新安江模型雖然是概念性半分布式水文模型，但對實際過程進(jìn)行合理概化，使得其在土壤含水量較高的地區(qū)得以廣泛應(yīng)用. 兩個模型的結(jié)構(gòu)對比結(jié)果見表5.

表5 WRF-Hydro模型和新安江模型結(jié)構(gòu)對比*

*次表面徑流與壤中流的概念相似，基流與地下徑流的概念相似.

3 結(jié)果分析

鑒于在陳河流域內(nèi)對WRF-Hydro模型與新安江模型模擬結(jié)果進(jìn)行對比的目的，本部分以率定期和檢驗期為基礎(chǔ)分別對它們的模擬結(jié)果和評價指標(biāo)進(jìn)行分析. 此外發(fā)現(xiàn)文獻(xiàn)[13]推薦的WRF-Hydro模型的計算步長不適用于陳河流域，需要重新尋找使模擬結(jié)果穩(wěn)定的計算步長，并分析造成這一現(xiàn)象的可能原因. 本研究中WRF-Hydro模型為3.0版本，陸面模式網(wǎng)格分辨率為1 km，匯流網(wǎng)格為100 m. 如無特殊聲明，本文所涉及單位及符號采用國際單位制(SI).

3.1 模型計算步長選取對計算結(jié)果的影響

由于需使用差分格式或一階牛頓迭代法對坡面淺水波方程組與河道的擴散波方程組求近似解，所以WRF-Hydro模型的匯流網(wǎng)格的計算步長選取是影響模型計算穩(wěn)定性的重要因素. 該計算步長的選取往往與網(wǎng)格空間分辨率有關(guān)，且對模擬結(jié)果的穩(wěn)定性影響較大. 尤其對水深較深的河道洪水演算而言，模擬結(jié)果受時間計算步長影響更大[13]. 為了保證模型的計算穩(wěn)定性，WRF-Hydro模型用戶手冊[13]給出了不同匯流網(wǎng)格分辨率對應(yīng)的計算步長. 如100 m的匯流網(wǎng)格對應(yīng)的時間計算步長為6 s. 然而，經(jīng)數(shù)值實驗發(fā)現(xiàn)該時間步長并不適用于陳河流域，如圖2中紅色虛線所示. 因此，需要在建議的基礎(chǔ)上進(jìn)一步縮小計算步長. 以兩場模擬洪水060926和110804為例，圖2描述這兩場洪水的計算步長分別取6、3和1 s時，WRF-Hydro模型的模擬結(jié)果. 當(dāng)Δt=6 s時，兩場洪水在洪峰段與退水段均出現(xiàn)數(shù)值波動現(xiàn)象；Δt=3 s時退水段模擬結(jié)果顯著提升，但在洪峰段仍偶爾存在波動；Δt=1 s時洪水過程表現(xiàn)平滑. 注意到110804洪水在Δt取3 s和6 s時相對于1 s洪峰整體向后坦化，隨后又于8月5日20:00時陡落，但洪量整體保持一致. 這說明若Δt選取不合適，甚至可能對洪峰模擬造成較大偏差，增加模型不確定性.

地形起伏較大、坡度和比降較陡，可能是默認(rèn)時間空間步長關(guān)系在陳河流域不適用的主要原因. 當(dāng)坡度或河道比降較大時，導(dǎo)致重力沿坡面的分力對洪水波運動影響較大. 水量相同的情況下導(dǎo)致流速變大，洪水波傳播速度c可能變快. 根據(jù)Courant數(shù)的定義Cn=cΔt/Δx，如果波速增大且保持空間分辨率Δx不變，則計算步長Δt需要適當(dāng)減小以保證Courant數(shù)數(shù)值不變. 進(jìn)一步說明在陳河流域適當(dāng)減小計算步長的必要性. 綜上，在陳河流域中利用當(dāng)前版本的WRF-Hydro模型進(jìn)行洪水模擬時，若要與100 m的匯流網(wǎng)格相對應(yīng)，計算步長須取1 s才能保證模型計算的穩(wěn)定性.

圖2 模型計算步長Δt的選取(模擬結(jié)果輸出步長：1 h)Fig.2 The determination of values of calculation step Δt (the time step of simulation results: 1 h)

3.2 率定期結(jié)果分析

本研究中率定期共由6場洪水構(gòu)成，即2003年的3場和2006年的3場，兩個模型的評價指標(biāo)結(jié)果見表6. WRF-Hydro模型的徑流深合格率和洪峰合格率均為66.7%，而新安江模型均為83.3%. 在2003年WRF-Hydro模型的模擬效果整體較好，NSE均值為0.91，而2006年模擬效果欠佳，NSE均值僅為0.03. 總體而言，率定期內(nèi)WRF-Hydro模型各場洪水均方根誤差均值為49.3 m3/s，新安江模型為53.3 m3/s. 在66.7%(6場中的4場)的洪水中，前者的均方根誤差小于后者. 說明前者的模擬結(jié)果誤差較小，且NSE表現(xiàn)尚可；后者模擬效果相對穩(wěn)定. 選取030828、030916、060827、060925為代表洪水，其過程線如圖3所示.

表6 WRF-Hydro模型和新安江模型在率定期的模擬結(jié)果*

*W代表WRF-Hydro模型，X代表新安江模型.

圖3 兩個模型在率定期內(nèi)4場代表洪水的對比(輸入輸出數(shù)據(jù)時間步長：1 h)Fig.3 The simulation results of four events in calibration period of two models (the time step of data: 1 h)

030828和030916號洪水的徑流深和洪峰均合格，而且兩個模型的模擬效果較好. WRF-Hydro模型模擬洪峰誤差略大于新安江模型，但AErd和NSE整體優(yōu)于新安江模型. 對于多峰洪水030828(圖3a)，由起漲段至第一個洪峰處，WRF-Hydro模型與觀測過程較為貼合，相關(guān)系數(shù)為0.9578；而新安江模型模擬過程出現(xiàn)較早的抬升，相關(guān)系數(shù)為0.9409. 兩個模型模擬洪峰偏大15%～20%，表現(xiàn)較接近. 最后的退水段，新安江模型退水幅度過大，模擬水量偏小. 對于單峰洪水030916(圖3b)，由起漲段至洪峰處，WRF-Hydro模型仍然在起漲時間和趨勢上模擬較好，相關(guān)系數(shù)為0.9834；而新安江模型再次出現(xiàn)較早抬升，相關(guān)系數(shù)為0.9656. 洪峰處WRF-Hydro模型模擬洪峰偏大7.6%；而新安江模型偏小7.3%. 退水階段與前場洪水相似，新安江模型低估退水段水量. 綜合兩場洪水發(fā)現(xiàn)WRF-Hydro模型對洪水過程的起漲過程和退水過程的模擬較好，其中洪水起漲點與實測過程幾乎重合；而新安江模型的洪水起漲點比實測過程提前8 h，且有高估漲洪段水量、低估退水段水量的特點.

對于洪水060827，結(jié)合圖3c發(fā)現(xiàn)，兩個模型均未成功模擬洪水過程. 根據(jù)洪水起漲點誤差和退水段的相關(guān)系數(shù)，WRF-Hydro模型在洪水起漲時刻和退水過程與實測過程相當(dāng)接近，但模擬洪峰大幅偏大. 雖然新安江模型對洪峰高估程度略小于WRF-Hydro模型，但其對整個洪水過程的模擬效果很不理想. 對于洪水060925(圖3d)，WRF-Hydro模型的初始流量為56.8 m3/s，是實測流量(18.0 m3/s)的3倍；新安江模型的初始流量為26.6 m3/s，是實測流量的近1.5倍. 在整個洪水過程中，WRF-Hydro模型模擬的洪量和洪峰顯著偏大；新安江模型在起漲段和洪峰處偏大，在退水段偏小，使徑流深誤差較小，僅為0.7 mm. 這兩場洪水模擬效果偏差的原因?qū)⒃诘?.3節(jié)討論.

綜合率定期所有洪水，發(fā)現(xiàn)雖然WRF-Hydro模型的徑流深合格率低于新安江模型，但在4場徑流深模擬較好的洪水中(030828、030903、030916和060903)，有3場NSE高于新安江模型. 而且該模型善于刻畫洪水的細(xì)節(jié)，例如率定期洪水起漲點誤差均在3 h之內(nèi)，而新安江模型在4～10 h間. 新安江模型整體表現(xiàn)較穩(wěn)定，有較高的徑流深合格率和洪峰合格率，66.7%(6場中的4場)的NSE大于0.8. 事實上，結(jié)合圖3可知，新安江模型在洪水起漲段和退水段均存在偏差，起漲段為正偏差，退水段為負(fù)偏差. 因此，在計算整體的徑流深誤差時部分偏差被抵消，可能使其在洪量方面表現(xiàn)較好.

3.3 檢驗期結(jié)果分析

檢驗期由2010年和2011年共10場洪水組成，兩個模型的評價指標(biāo)結(jié)果見表7. WRF-Hydro模型的徑流深合格率和洪峰合格率均分別為60.0%和50.0%，而新安江模型均為50.0%. 率定期內(nèi)WRF-Hydro模型的均方根誤差均值為115.7 m3/s，新安江模型為80.0 m3/s. 其中60%(10場中的6場)的WRF-Hydro模型的均方根誤差小于新安江模型，且這些洪水多為兩個指標(biāo)均合格的洪水. 在兩個指標(biāo)均合格的洪水中，WRF-Hydro模型的均方根誤差均值為54.8 m3/s，新安江模型為68.6 m3/s，且前者的AErd、REfp和NSE整體優(yōu)于后者. 其他洪水多為不合格洪水，其中前者的均方根誤差的均值為207.1 m3/s，后者為97.1 m3/s. 說明WRF-Hydro模型在洪峰洪量表現(xiàn)好的洪水中誤差較小，NSE較高；在其他洪水中，偏離程度較明顯. 新安江模型對峰現(xiàn)時間的刻畫較好，表現(xiàn)相對穩(wěn)定；對洪峰洪量誤差較大的洪水模擬效果也較差，但優(yōu)于WRF-Hydro模型. 選取100812、100823、110728、110909為代表洪水，其過程線如圖4所示.

表7 WRF-Hydro模型和新安江模型在檢驗期的模擬結(jié)果*

*W代表WRF-Hydro模型，X代表新安江模型.

兩個模型對100823和110909號洪水模擬結(jié)果較好(圖4a，圖4b). 在洪水100823的起漲段，WRF-Hydro模型與觀測過程較為貼合，相關(guān)系數(shù)為0.9876；新安江模型高估洪水的初始狀態(tài)，從而模擬的起漲過程水量偏大，相關(guān)系數(shù)為0.9285. 兩個模型對洪峰模擬效果都很出色，WRF-Hydro模型的REfp為2.1%，新安江模型的REfp為-1.9%. 退水段，WRF-Hydro模型表現(xiàn)出色，相關(guān)系數(shù)為0.9981；新安江模型在退水前期和后期略微高估洪量，在退水中期低估洪量，相關(guān)系數(shù)為0.9872. 對于洪水110909，WRF-Hydro模型和新安江模型的NSE分別達(dá)到了0.96和0.92，模擬結(jié)果相當(dāng)出色. 在洪水起漲段，WRF-Hydro模型受前場洪水退水過程影響，模擬流量過程偏大，相關(guān)系數(shù)為0.9819；新安江模型對洪量模擬較好，相關(guān)系數(shù)為0.9468. WRF-Hydro模型模擬洪峰較準(zhǔn)確，REfp和AEpt分別為9.0%和0 h. 在退水段，WRF-Hydro模型仍與觀測過程較接近，相關(guān)系數(shù)為0.9965；而新安江模型低估退水中后期的洪量. 此外，WRF-Hydro模型甚至能較準(zhǔn)確模擬出9月14日2:00時的小洪水波，表現(xiàn)該模型在退水段的優(yōu)秀性能，這可能是該模型可利用的信息較豐富導(dǎo)致的.

兩個模型對另外兩場洪水的模擬效果欠佳. 100812和110728號洪水與060827號洪水相似，WRF-Hydro模型模擬的洪量和洪峰仍然大幅偏大；新安江模型也大幅偏大，但偏離程度小于WRF-Hydro模型. 這兩場洪水模擬效果偏差的原因也將在第4.3節(jié)討論.

圖4 兩個模型在檢驗期內(nèi)4場代表洪水的對比(輸入輸出數(shù)據(jù)時間步長：1 h)Fig.4 The simulation results of four events in verification period of two models (the time step of data: 1 h)

綜合檢驗期所有洪水，發(fā)現(xiàn)WRF-Hydro模型的徑流深合格率和洪峰合格率平于甚至優(yōu)于新安江模型. 在多數(shù)兩個指標(biāo)均合格的洪水中(100820、100822和110909)，前者的均方根誤差小于后者；在其他洪水中，前者模擬過程的誤差很大. WRF-Hydro模型檢驗期起漲時刻誤差均值為1.3 h，而新安江模型為5 h；在70%(10場中的7場)的洪水中，前者模擬的退水段的相關(guān)系數(shù)大于后者. 說明WRF-Hydro模型善于模擬洪水的細(xì)節(jié)；新安江模型在整體上的表現(xiàn)非常穩(wěn)定，很少出現(xiàn)大幅度的偏差.

3.4 模擬結(jié)果對比

在所有洪水事件中，WRF-Hydro模型的徑流深合格率和洪峰合格率分別為62.5%和56.3%，均方根誤差均值為90.8 m3/s；新安江模型徑流深合格率和洪峰合格率均為62.5%，均方根誤差均值為70.0 m3/s. 這說明WRF-Hydro模型的徑流深和洪峰合格率平于或略低于新安江模型，且均方根誤差均值稍大. 在徑流深和洪峰指標(biāo)均合格的洪水中(030828、030902、030916、060903、100820、100822和110909；共7場)，WRF-Hydro模型均方根誤差均值為43.5 m3/s，NSE的均值為0.87，其中在57.1%(7場中的4場)的事件中，NSE達(dá)到 0.9 以上；新安江模型均方根誤差均值為55.4 m3/s，NSE的均值為0.84，其中在42.9%(7場中的3場)的事件中，NSE達(dá)到0.9以上. 在其他洪水中，WRF-Hydro模型均方根誤差均值為127.6 m3/s；新安江模型均方根誤差均值為81.3 m3/s. 由此發(fā)現(xiàn)WRF-Hydro模型在兩個指標(biāo)均合格的洪水中均方根誤差均值比新安江模型小21.5%，而在其他洪水中，前者均方根誤差均值比后者大56.2%.

在兩個指標(biāo)均合格的洪水事件中，可能WRF-Hydro模型需要的資料信息較多，相對于新安江模型更加善于模擬洪水過程線的細(xì)節(jié)，使模擬洪水的誤差較小，NSE較高. 但是，在其他表現(xiàn)效果較差的事件中，大量信息中的不確定性可能對WRF-Hydro模型的模擬過程產(chǎn)生較大影響，使模擬洪水的偏離幅度較大. 新安江模型反而因為未使用大量的資料，模擬過程的偏離幅度小于WRF-Hydro模型. 綜上，WRF-Hydro模型可以較好地模擬洪水過程，刻畫洪水細(xì)節(jié)，但可能對資料質(zhì)量較為敏感. 新安江模型誤差整體較小，模擬效果較穩(wěn)定. 此外，WRF-Hydro模型在洪水起漲處模擬較好，其漲洪點誤差均值為-1 h，而新安江模型為-5.7 h. 這可能是由于WRF-Hydro模型利用地形、糙率等信息獲得了更為精細(xì)的模擬結(jié)果.

4 對比討論

本部分主要從模型參數(shù)和輸入輸出信息對WRF-Hydro模型和新安江模型進(jìn)行對比討論，最后對模擬效果較差洪水進(jìn)行分析.

4.1 參數(shù)對比

WRF-Hydro模型可調(diào)參數(shù)較多，主要包含植被參數(shù)(如冠層持水能力等)、土壤特性參數(shù)(如水力傳導(dǎo)度)、下墊面參數(shù)(如地表糙率等)及匯流參數(shù)(如河道糙率等)[15, 29]. 三水源新安江模型約有17個參數(shù)，主要分為產(chǎn)流參數(shù)(如流域平均張力水容量等)、蒸散發(fā)參數(shù)(如流域蒸散發(fā)折算系數(shù)等)、分水源參數(shù)(如表層自由水蓄水容量等)和匯流參數(shù)(如河網(wǎng)消退系數(shù)等)[6,8,24].

4.1.1 產(chǎn)流參數(shù) 兩個模型在概化流域土壤時有較大差異. 對于大尺度的天氣或陸面模式，通常只考慮包氣帶中對天氣過程影響較大的部分. 如在Noah陸面模式中，認(rèn)為地表至地下2 m之間的土層對天氣過程影響較顯著，而2 m以下的土層對模擬的影響可忽略不計. 而WRF-Hydro模型是基于Noah陸面模式發(fā)展而來，繼承了這一特點，即只考慮近地面2 m深的土層內(nèi)水量分配情況，其余土層內(nèi)的水視為地下徑流. 但是對于小尺度的水文模擬而言，默認(rèn)的次表面層厚度可能不適用于中小流域. 需要經(jīng)過土層厚度乘子ZSOILFAC對其進(jìn)行調(diào)節(jié)，以控制模擬過程中土層的含水量. 此外，幾乎每層土層內(nèi)的土壤水都會受滲透率REFKDT和蒸散發(fā)影響，各層水量平衡利用Richard方程計算[26-27]. 而對于新安江模型，其著重研究包氣帶中的水量平衡，并將土壤水劃分為張力水和表層自由水，分別由張力水蓄水容量WM和自由水蓄水容量SM控制. WM影響土壤含水量和產(chǎn)流計算，SM影響地表徑流、壤中流、地下徑流3種徑流成分的比例[6,8].

此外，研究發(fā)現(xiàn)經(jīng)過ZSOILFAC調(diào)節(jié)后的次表面土層厚度和新安江模型反推的包氣帶厚度有較好的關(guān)系. 后者的計算過程如下. 由于張力水蓄水容量存在空間差異，故對于任意匯流網(wǎng)格，假設(shè)其張力水蓄水容量為WM(i,j)，且可由式(5)估計[6]：

WM=(θfc-θwp)La

(5)

式中，θfc和θwp分別為田間持水量和凋萎含水量，均為體積含水量，La為當(dāng)前網(wǎng)格的包氣帶厚度(mm). 根據(jù)(5)式，利用WM的最大值和最小值即可反推出包氣帶厚度的分布范圍.WM的最大值為流域最大張力水蓄水容量Wmax，可由式(6)[6,24]計算：

Wmax=(1+B)WM

(6)

式中，B為張力水蓄水容量曲線的指數(shù)，取0.3；WM取180 mm(表4). 由式(6)可知，Wmax=234 mm. 根據(jù)(5)式及壤土的屬性(θfc=0.3,θwp=0.1)，與Wmax相對應(yīng)的La的最大值為1.17 m. 根據(jù)流域變動產(chǎn)流理論[6,9]，在濕潤地區(qū)河道附近包氣帶較薄，河流上游坡地的包氣帶相對較厚[30-31]. 故假設(shè)河道附近網(wǎng)格幾乎飽和，相應(yīng)的包氣帶厚度為0. 綜上，由新安江模型反推的流域包氣帶厚度分布在0～1.17 m內(nèi)，均值為0.8 m. 而次表面層經(jīng)ZSOILFAC調(diào)節(jié)后，厚度為0.4 m，與包氣帶厚度均值接近. 由此說明，經(jīng)ZSOILFAC作用后的次表面層厚度與新安江模型反推所得的包氣帶厚度有較好的一致性，進(jìn)一步說明在中小流域中利用ZSOILFAC調(diào)節(jié)次表面層厚度的合理性與必要性.

4.1.2 匯流參數(shù) WRF-Hydro模型主要利用擴散波方程，并結(jié)合OVROUGHTFAC和MANNFAC兩個參數(shù)考慮地表和河道的調(diào)蓄作用；次表面徑流主要可以由側(cè)向水力傳導(dǎo)度LKSATFAC調(diào)節(jié)[16]；地下徑流調(diào)蓄作用可利用概念水箱模型[13]調(diào)節(jié). 新安江模型在流域匯流方面主要利用線性水庫考慮流域調(diào)蓄作用，如壤中流消退系數(shù)(CI)和地下水消退系數(shù)(CG)；子流域河網(wǎng)匯流利用滯后演算法，如CS和滯時L[6,24]. WRF-Hydro模型的OVROUGHTFAC和MANNFAC與新安江模型的CS均可以調(diào)整洪水形狀，有相似的作用.

4.2 輸入輸出信息對比

兩個模型在各個過程的輸入輸出信息見表8. 作為基于物理的分布式水文模型，WRF-Hydro模型需要大量的氣象強迫和下墊面屬性數(shù)據(jù)，對資料依賴性和計算機性能要求較高. 基于高質(zhì)量的強迫數(shù)據(jù)及高性能服務(wù)器，該模型可以很好地描述各個水文氣象變量的空間不均勻性，獲得較精細(xì)的輸出結(jié)果. 作為概念性半分布式水文模型，新安江模型合理地簡化實際產(chǎn)匯流情況，如分別利用自然子流域分割法和蓄水容量曲線考慮降水和土壤水的空間分布不均勻性[6,8,24]. 該模型能在保證一定模擬精度的基礎(chǔ)上降低計算量. 但正是由于對實際過程簡化，無法獲得精細(xì)土壤水空間分布及各個河道的流量. 兩個模型輸入信息量和計算量的巨大差別導(dǎo)致計算精度存在一定差別.

表8 WRF-Hydro模型和新安江模型輸入輸出對比

4.3 表現(xiàn)較差洪水分析

表9對模擬結(jié)果較差的洪水進(jìn)行歸納. 根據(jù)成因不同，可將洪水分為兩類：060827、100812和110728為第1類，060925和110731為第2類. 對于第1類洪水，兩個模型模擬的洪峰和洪量均偏大，而實測徑流量與降雨量的比值較小，均小于0.2，說明大部分降雨被流域保持. 可能有以下兩點原因：1)因為這類洪水事件前期缺少降水，洼地、庫塘閘壩未蓄滿，使這類洪水的實測值較小. 2)由于缺少有效降水，可能導(dǎo)致土層相對干燥、腐殖質(zhì)層松軟、地下水位偏低、包氣帶偏厚. 進(jìn)行參數(shù)率定時未考慮這些洪水事件產(chǎn)流的差異性，使模型低估了實際的下滲率和土壤調(diào)蓄能力，導(dǎo)致模擬過程偏大. 此外，該類型洪水在率定期內(nèi)只包含060827，針對該問題考慮較少. 綜上，實測值較小且模擬過程偏大，導(dǎo)致兩個模型模擬的徑流量與洪峰偏大. WRF-Hydro模型利用更多的降水空間分布的下墊面信息，使其更易受影響.

表9 WRF-Hydro模型和新安江模型主要不合格洪水統(tǒng)計*

*W代表WRF-Hydro模型，X代表新安江模型，√代表合格，×代表不合格.

第2類洪水的實測徑流量與降雨量的比值均大于0.35. 雖然兩個模型模擬的洪水過程偏大，但新安江模型的表現(xiàn)整體較好. 對于洪水060925，WRF-Hydro模型的參數(shù)可能低估了洪水開始前的包氣帶厚度、下滲率和糙率，導(dǎo)致洪水初始時刻模擬流量偏大近3倍. 當(dāng)降雨達(dá)到極值后，仍以較小的雨量持續(xù)近5天，進(jìn)一步促使WRF-Hydro模型模擬結(jié)果偏高. 而新安江模型的初始狀態(tài)源于日模型，對前期水量和徑流成分模擬較好. 洪水110731與洪水060925的原因相似.

總體而言，上述洪水模擬效果較差的可能原因是流域初始狀態(tài)和其他洪水有較大差異，如包氣帶厚度、土壤含水量、下滲率、地表糙率和河道糙率等. 這種情況下仍使用相同的參數(shù)，可能導(dǎo)致兩個模型模擬的徑流量與洪峰偏大. 此外，WRF-Hydro模型需要大量的資料輸入，這些資料的誤差可能在模擬流域狀態(tài)和洪量過程時傳遞，增加不確定性. 新安江模型需要資料較少，由資料誤差引起不確定性的概率稍低. 改善這類問題主要有3點：1) 利用高質(zhì)量高空間分辨率的資料或融合數(shù)據(jù)驅(qū)動WRF-Hydro模型，提高模擬結(jié)果的精度與可靠性；2) 盡可能地延長預(yù)熱期，把握好計算成本和預(yù)熱期長度的平衡；3) 增加率定期洪水場次數(shù)，在計算成本和率定期洪水場次數(shù)之間尋找平衡.

5 結(jié)論與展望

在半濕潤區(qū)中小流域陳河流域中，利用GLDAS非降水?dāng)?shù)據(jù)場與IDW降雨場的結(jié)合數(shù)據(jù)驅(qū)動WRF-Hydro模型，利用站點的降雨數(shù)據(jù)和蒸發(fā)皿觀測數(shù)據(jù)驅(qū)動新安江模型. 探究WRF-Hydro模型在陳河流域的表現(xiàn)及適用性，對比分析兩個模型的結(jié)構(gòu)、主要參數(shù)、輸入輸出和模擬結(jié)果，并對模擬效果較差洪水進(jìn)行討論. 研究表明：

1)在陳河流域中，需要在建議的基礎(chǔ)上進(jìn)一步縮小計算步長的值，以滿足穩(wěn)定性的要求. 本文中匯流網(wǎng)格的空間分辨率為100 m，對應(yīng)的計算時間步長應(yīng)取1 s. 原因可能是陳河流域地形起伏較大、坡度和比降較陡.

2)借鑒WRF-Hydro模型的其他乘子參數(shù)，在中小流域中引入土層厚度乘子ZSOILFAC，對整個次表面層的厚度等比縮放，使次表面層與包氣帶建立聯(lián)系. 經(jīng)過ZSOILFAC調(diào)節(jié)后的次表面土層厚度和新安江模型反推的包氣帶厚度空間分布場有較好的關(guān)系，進(jìn)一步說明在中小流域中利用ZSOILFAC調(diào)節(jié)次表面層厚度的合理性與必要性.

3)WRF-Hydro模型能夠較好地模擬實測洪水過程，捕捉洪水的細(xì)節(jié). 其漲洪點誤差均值為-1 h，比新安江模型小82.5%. 這可能是WRF-Hydro模型可利用大量信息，提高洪水模擬效果，但同時可能帶來較大的不確定性. 該模型與新安江模型性能相近，表現(xiàn)出在半濕潤地區(qū)中小流域內(nèi)實際應(yīng)用的潛力.

4)對于所有洪水，WRF-Hydro模型的徑流深合格率和洪峰合格率分別為62.5%和56.3%，均方根誤差均值為90.8 m3/s；新安江模型均為62.5%，均方根誤差均值為70.0 m3/s. 說明WRF-Hydro模型的徑流深和洪峰合格率平于或略低于新安江模型，且均方根誤差均值稍大；新安江模型整體誤差小于WRF-Hydro模型，模擬效果相當(dāng)穩(wěn)定.

5)在徑流深和洪峰指標(biāo)均合格的洪水事件中，WRF-Hydro模型均方根誤差均值比新安江模型小21.5%；但在其他事件中，前者均方根誤差均值比后者大56.2%. 可能是WRF-Hydro模型對信息質(zhì)量較為敏感，信息質(zhì)量較好時善于模擬洪水細(xì)節(jié)；反之，資料的不確定性以及參數(shù)與某些洪水實際情況不匹配將對模擬過程產(chǎn)生較大影響.

本研究使用GLDAS資料的分辨率相對較粗，下一步將對多種資料或組合資料進(jìn)行對比研究，如WRF模擬資料、降水融合資料等，以尋找更適用于WRF-Hydro模型的資料，為之后的實時預(yù)報提供理論基礎(chǔ). 研究洪水場次較少，后期研究需要適當(dāng)增加洪水場次. 而針對討論部分提出的率定參數(shù)可能與某些洪水實際狀態(tài)不匹配的問題，需要在預(yù)熱期長度、率定期洪水場次數(shù)和計算成本之前尋找平衡. WRF-Hydro模型作為新一代陸氣耦合的分布式水文模型系統(tǒng)，在水文模擬、洪水預(yù)報和水資源評估等方面有較大應(yīng)用的潛力，同時也為水文氣象研究提供新的思路.

致謝：感謝NCAR的喻煒博士對WRF-Hydro模型參數(shù)調(diào)整的指導(dǎo)與建議.