周方錄，黃金柏，王斌，檜谷治

（1.東北農(nóng)業(yè)大學(xué)水利與建筑學(xué)院，哈爾濱150030；2.大慶市松嫩工程管理處，黑龍江大慶163300；3.揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院，江蘇揚(yáng)州225008；4.鳥(niǎo)取大學(xué)大學(xué)院工學(xué)研究科，日本鳥(niǎo)取680-8552）

兩種地形條件下小尺度流域降雨-徑流數(shù)值模擬

周方錄1,2，黃金柏3*，王斌1，檜谷治4

（1.東北農(nóng)業(yè)大學(xué)水利與建筑學(xué)院，哈爾濱150030；2.大慶市松嫩工程管理處，黑龍江大慶163300；
3.揚(yáng)州大學(xué)水利與能源動(dòng)力工程學(xué)院，江蘇揚(yáng)州225008；4.鳥(niǎo)取大學(xué)大學(xué)院工學(xué)研究科，日本鳥(niǎo)取680-8552）

分別以位于黃土高原北部溝壑區(qū)的六道溝流域和日本鳥(niǎo)取縣境內(nèi)的山區(qū)流域-Bukuro流域?yàn)檠芯繀^(qū)，對(duì)兩個(gè)研究區(qū)域土壤垂直剖面模型化，以運(yùn)動(dòng)波基礎(chǔ)方程式結(jié)合GIS-ArcMap構(gòu)建分布式流域降雨-徑流數(shù)值模型；以對(duì)觀測(cè)徑流過(guò)程的數(shù)值模擬檢驗(yàn)?zāi)Ｐ蛯?shí)用性。結(jié)果表明，模型計(jì)算精度較高，數(shù)值模擬結(jié)果誤差＜3%，觀測(cè)流量過(guò)程和數(shù)值計(jì)算結(jié)果之間的皮爾森相關(guān)系數(shù)在0.90以上，為地表徑流準(zhǔn)確推算提供科學(xué)計(jì)算方法。

流域；地形；分布式水文模型；降雨-徑流；數(shù)值模擬

分布式水文模型具有流域物理基礎(chǔ)和分布式參數(shù)[1-2]，模擬功能強(qiáng)大，發(fā)展迅速[3-4]。隨著計(jì)算機(jī)技術(shù)、地理信息系統(tǒng)（GIS）和衛(wèi)星遙感技術(shù)（RS）發(fā)展，獲取和描述流域下墊面空間分布信息技術(shù)日益完善，分布式水文模型構(gòu)建技術(shù)日漸成熟。與數(shù)字高程模型（DEM：digital elevation model）相結(jié)合的數(shù)字分布式水文模型成為國(guó)際水文科學(xué)研究領(lǐng)域熱點(diǎn)[5-6]，分布式水文模型發(fā)展表現(xiàn)在向不同尺度流域、多功能、多學(xué)科、模塊化模型系統(tǒng)發(fā)展[7-8]。

我國(guó)以分布式水文模型圍繞中小尺度流域開(kāi)展研究，如唐麗莉等利用分布式水文模型DHSVM（Distributed hydrology soil vegetation model）對(duì)蘭江流域徑流變化的模擬試驗(yàn)[9]；張利平等將分布式水文模型VIC（Variable infiltration capacity）和SWAT（Soil and water assessment tool）模型應(yīng)用于白蓮河流域，探討模型在中小流域的適用性[10]；張會(huì)蘭等利用BPCC（Basic pollution calculated center）分布式水文模型，研究岷江鎮(zhèn)江關(guān)流域氣候波動(dòng)和覆被變化對(duì)流域徑流影響[11]；陳瑩等以長(zhǎng)江三角洲地區(qū)太湖上游西苕溪流域?yàn)檠芯繀^(qū)，借助分布式水文模型SWAT對(duì)流域降雨徑流過(guò)程進(jìn)行模擬[12]；楊?yuàn)檴櫟葘WAT模型應(yīng)用于臥虎山水庫(kù)流域徑流模擬[13]；李致家等將物理基礎(chǔ)分布式流域水文模型TOPKAPI與新安江模型應(yīng)用于呈村流域洪水模擬計(jì)算[14]。

由以上分析可知，我國(guó)近期圍繞小流域分布式水文模型研究多為現(xiàn)有國(guó)外模型的應(yīng)用性研究，依托我國(guó)小流域物理性條件自開(kāi)發(fā)分布式水文模型的研究相對(duì)缺乏。

針對(duì)當(dāng)前小尺度流域分布式水文模型研究現(xiàn)狀，本研究選取兩個(gè)位于不同地形條件下的小流域?yàn)檠芯繀^(qū)，基于GIS依托研究流域的地形和基礎(chǔ)水文地質(zhì)條件自開(kāi)發(fā)分布式降雨-徑流模型，以期為小尺度流域分布式水文模型深入研究提供技術(shù)參考，為推進(jìn)中小尺度流域分布式水文數(shù)值模型平臺(tái)構(gòu)建提供基礎(chǔ)數(shù)據(jù)。

1 材料與方法

1.1 小尺度流域概念

對(duì)于小尺度流域，尚無(wú)嚴(yán)格定義。小流域通常指二、三級(jí)支流以下以分水嶺和下游河道出口斷面為界，集水面積在100 km2以下相對(duì)獨(dú)立和封閉的自然匯水區(qū)域，水利上通常指面積＜1 000 km2或河道基本上在一個(gè)縣屬范圍內(nèi)的流域。按照全國(guó)科學(xué)技術(shù)名詞審定委員會(huì)的定義，我國(guó)小流域面積均以＜100 km2為限，以5～30 km2占多數(shù)[15-16]。參考國(guó)內(nèi)及國(guó)際上水利類相關(guān)文獻(xiàn)，小尺度流域一般面積數(shù)量級(jí)為50km2[17-19]。參照上述有關(guān)研究或相關(guān)部門對(duì)小流域尺度定義，根據(jù)多年從事分布式水文模型研究過(guò)程中對(duì)流域尺度與流域基礎(chǔ)水文地質(zhì)條件之間關(guān)系的探究，本研究所指小尺度流域，流域面積數(shù)量級(jí)≤50km2，流域水文地質(zhì)條件相對(duì)穩(wěn)定，即在同一流域內(nèi)部，土壤垂直剖面含水層及土層分布，在數(shù)量上保持一致，只是尺度（含水層及土層厚度）上不同。

1.2 研究區(qū)調(diào)查

研究分別選取位于陜西省神木縣的北部黃土高原溝壑區(qū)六道溝流域（東經(jīng)110°21'～110°23'，北緯38°46'～38°51'）及位于日本鳥(niǎo)取縣境內(nèi)的山區(qū)流域-Bukuro流域（東經(jīng)134°12'～134°26'，北緯35° 25'～35°29'）為研究區(qū)。

在各研究流域水文觀測(cè)及土壤垂直剖面分層情況（基礎(chǔ)水文地質(zhì)條件）的實(shí)際調(diào)查，降雨觀測(cè)采用雨量計(jì)[型號(hào)：7852M-L10，尺寸：φ165×240H（mm）]。地表徑流（河道水位）采用自動(dòng)記數(shù)水位計(jì)觀測(cè)（六道溝流域：水位計(jì)型號(hào)，KADEC21-MZPT，制造商，KONA System Co.,Ltd.；Bukuro流域，水位計(jì)型號(hào)，HM-910-02-309，制造商：Sensez Co.,Ltd.）。六道溝流域地表徑流觀測(cè)點(diǎn)位于流域上游的一條支溝內(nèi)（控制面積0.1 km2），Bukuro流域地表徑流觀測(cè)點(diǎn)控制流域面積為31.8 km2。在六道溝流域水文數(shù)據(jù)觀測(cè)期間內(nèi)，進(jìn)行氣象數(shù)據(jù)同期觀測(cè)。

1.3 模型構(gòu)建

1.3.1 土壤垂直剖面模型化

研究需要對(duì)各研究流域土壤垂直剖面分層情況模型化，構(gòu)建自地面開(kāi)始至地下某一含水層土壤垂直剖面模型，承載雨水降落到地面后在垂直方向的運(yùn)動(dòng)過(guò)程，使水的運(yùn)動(dòng)方式受下墊面實(shí)際物理?xiàng)l件約束。根據(jù)對(duì)研究流域地形及基礎(chǔ)水文地質(zhì)條件調(diào)查結(jié)果分析，篩選土壤垂直剖面結(jié)構(gòu)特征參數(shù)，構(gòu)建各研究流域土壤垂直剖面模型如圖1和2所示。

圖1 六道溝流域土壤垂直剖面模型化及局部地形寫真Fig.1Modeling of soil vertical profile of Liudaogou Basin and topographic photo(partial view)

圖2Bukuro流域土壤垂直剖面模型化及地形寫真Fig.2Modeling of soil vertical profile of Bukuro Basin and topographic photo(partial view)

六道溝流域土壤垂直剖面模型由坡面和河道構(gòu)成，因長(zhǎng)期受水蝕和風(fēng)蝕交錯(cuò)作用，表層土壤（厚度約10 cm）沙化嚴(yán)重，與其土壤孔隙度、滲透系數(shù)等物理性參數(shù)存在較明顯差異，在雨季某些降雨條件下，表層土壤底部和下層土壤過(guò)渡區(qū)域，會(huì)有短歷時(shí)的壤中流產(chǎn)生，因此，坡面區(qū)間自地表至砂巖層表面被認(rèn)為由兩層構(gòu)成。河道區(qū)間自河床表面至砂巖層表面被開(kāi)發(fā)成一層（見(jiàn)圖1a）。

根據(jù)土壤垂直剖面和含水層分布情況實(shí)際調(diào)查結(jié)果，Bukuro流域河道區(qū)間和坡面區(qū)間自地表至砂巖層表面均被開(kāi)發(fā)成兩層，坡面和河道第一層厚度不同，隨流域內(nèi)地理位置改變存在空間分異性，坡面和河道第二層厚度較均勻，全流域內(nèi)變化不大，厚度約2 m。該流域地表徑流和各含水層地下水（包括潛水和承壓含水層）流向均由坡面區(qū)間流向河道區(qū)間（見(jiàn)圖2a）。

1.3.2 河網(wǎng)及流域劃分

利用GIS-ArcMap，構(gòu)建各研究流域DEM及河網(wǎng)，因?yàn)辄S土高原六道溝流域溝壑密度較大，為清楚表達(dá)溝壑地形，流域DEM柵格尺寸采用5 m× 5 m，Bukuro河流域DEM柵格尺寸為50 m×50 m?；诟髁饔駾EM及河網(wǎng)，按照河網(wǎng)上各小流域空間連接關(guān)系（水流入關(guān)系）對(duì)流域分級(jí)和編號(hào)，生成分布式流域模型如圖3和4所示。

圖3 六道溝流域計(jì)算區(qū)域河網(wǎng)和各小流域的連接關(guān)系Fig.3River channel network of the calculation area of Liudaogou Basin and spatial connection of each small catchment

圖4Bukuro流域計(jì)算區(qū)域河網(wǎng)及各小流域的連接關(guān)系Fig.4River channel network of the calculation area of Bukuro Basin and spatial connection of each small catchment

1.3.3 基礎(chǔ)方程式

運(yùn)動(dòng)波模型被廣泛應(yīng)用于降雨-徑流過(guò)程的計(jì)算[20-21]，該模型以物理性參數(shù)描述流域及水流運(yùn)動(dòng)過(guò)程，可對(duì)地表徑流進(jìn)行準(zhǔn)確計(jì)算，且可通過(guò)達(dá)西公式描述和計(jì)算滲透及地下水徑流過(guò)程[22-23]，所以，本研究采用運(yùn)動(dòng)波理論基礎(chǔ)方程式構(gòu)建降雨-徑流計(jì)算方法，以坡面區(qū)間為例，給出計(jì)算公式如下：

地表徑流連續(xù)方程式：

式中，dt-計(jì)算的時(shí)間步長(zhǎng)（s）；dx-水流方向上計(jì)算的空間步長(zhǎng)（m）；h-水深（m）；q-單寬流量（m·s-1）；r-降雨（m·s-1）；f1-第一層土壤平均滲透速度（m·s-1）；R-水力半徑（m）；v-流速（m）；n-糙率（曼寧粗度系數(shù)）（s·m-1/3）；I-河道坡度；Q-流量（m3·s-1）；A-過(guò)水?dāng)嗝婷娣e（m2）。

滲透流（地下水）連續(xù)方程式（以第一層為例）：

式中，λ-有效孔隙率；hˉ-地下水（潛水）水深（m）；qˉ-地下水單寬流量（m2·s-1）；ET-蒸散發(fā)（m·s-1）；f2-第二層平均滲透速度（m·s-1）；k-橫向透水系數(shù)（m·s-1）；Iˉ-潛水含水層坡降，其他因子與前述相同。

第二層用于地下水計(jì)算的公式與第一層計(jì)算采用的基本公式相同，只是層號(hào)不同導(dǎo)致部分參數(shù)腳標(biāo)發(fā)生變化；河道與坡面采用同一組方程式計(jì)算，個(gè)別水流入或流出成分變化導(dǎo)致公式在形式上略有差別，相應(yīng)內(nèi)容在此略去。

1.3.4 初始條件和邊界條件

1.3.4.1 初始條件

因降雨-徑流計(jì)算從流域最末級(jí)流域的源點(diǎn)開(kāi)始，在此位置上，流域以分水線為邊界開(kāi)始產(chǎn)匯流，該處地表徑流水深（流量）為0。所以，計(jì)算開(kāi)始時(shí)刻，在各研究區(qū)最末級(jí)小流域（如圖3a，4a所示河網(wǎng)上各分布式小流域）源點(diǎn)處，地表徑流水深（流量）被設(shè)為0。

對(duì)于地下水初始水深的確定，在六道溝流域，河道區(qū)間或坡面區(qū)間，地下水（潛水）埋深較大（見(jiàn)圖1a），坡面區(qū)間自地面開(kāi)始至潛水自由表面距離在20 m以上，該研究區(qū)位于典型半干旱區(qū)，即使是降雨相對(duì)集中的雨季，也很少發(fā)生長(zhǎng)歷時(shí)強(qiáng)降雨事件，經(jīng)對(duì)2006～2010年實(shí)測(cè)降雨事件的特征分析，該流域發(fā)生的降雨事件多為短歷時(shí)降雨事件，歷時(shí)＜1 h次降雨事件達(dá)64%[24]，而短歷時(shí)降雨事件在產(chǎn)流期間尚未下滲到潛水含水層，即潛水含水層變動(dòng)不會(huì)對(duì)降雨徑流產(chǎn)生影響。因此在六道溝流域，計(jì)算開(kāi)始時(shí)刻潛水層初始水深，由溝道內(nèi)通過(guò)人工結(jié)合簡(jiǎn)單機(jī)械方法，通過(guò)鉆孔調(diào)查近似確定，即在溝道內(nèi)用土鉆開(kāi)挖小口徑水井（直徑5 cm），其深度至潛水含水層下部弱透水層表面，分別記下鉆孔至潛水層表面和底部深度，確定含水層深度。而在雨季某些降雨條件下，其松散表土和表土下硬黃土之間，會(huì)產(chǎn)生短歷時(shí)的壤中流，對(duì)降雨-徑流過(guò)程影響較大，確定表土中初始水深方法為，先利用環(huán)刀法對(duì)表層土取樣，然后烘干土樣，推算表土含水量，而后利用層厚和含水量關(guān)系，換算為表土中初始（計(jì)算開(kāi)始時(shí)刻）水深。

Bukuro流域用于數(shù)值計(jì)算的潛水（地下水）初始水深，通過(guò)鉆孔調(diào)查結(jié)合對(duì)計(jì)算區(qū)域內(nèi)現(xiàn)有潛水井水深測(cè)量方法確定。

1.3.4.2 邊界條件

在兩個(gè)研究流域，水流在垂直方向運(yùn)動(dòng)的邊界條件由地表徑流和第一層水深約束，以研究區(qū)坡面區(qū)間為例，自地表向第一層滲透速度由地表徑流深度h與滲透速度f(wàn)1關(guān)系通過(guò)計(jì)算確定，如當(dāng)?shù)乇韽搅魉頷i與計(jì)算時(shí)間步長(zhǎng)dt比值hi/dt≥f1時(shí)，實(shí)際滲透到第一層的速度為f1；當(dāng)?shù)乇韽搅魉頷i與計(jì)算時(shí)間步長(zhǎng)dt比值hi/dt＜f1時(shí)，實(shí)際滲透到第一層速度為hi/dt（小于f1）；地表徑流水深hi為0時(shí)，自地表向第一層的滲透速度為0。

對(duì)于地下水（滲透流），當(dāng)?shù)谝粚铀钆c時(shí)間步長(zhǎng)比值hˉ1/dt≥f2時(shí)，自第一層向第二層的入滲速度為f2；當(dāng)hˉ1/dt＜f2時(shí)，自第一層向第二層的入滲速度為hˉ1/dt；當(dāng)hˉ1=0時(shí)，由第一層向第二層的入滲速度為0。

1.3.5 參數(shù)率定

利用分布式水文模型對(duì)流域降雨-徑流過(guò)程計(jì)算時(shí)，需要大量的物理性參數(shù)，如各級(jí)流域坡面長(zhǎng)度、坡度，河道長(zhǎng)度、河道寬度和坡度，糙率，各層土壤滲透系數(shù)、橫向透水系數(shù)等。以上參數(shù)中，坡面與河道尺度參數(shù)及坡度，利用各流域數(shù)字高程模型（DEM）確定；各流域土壤垂直剖面分層情況如土層和含水層（潛水含水層）厚度及分布主要通過(guò)多點(diǎn)鉆孔調(diào)查方法率定。土壤水力學(xué)特性，特別是表層土壤水力學(xué)特性參數(shù)如滲透系數(shù)、透水系數(shù)、有效孔隙率等對(duì)產(chǎn)流和入滲過(guò)程有很大影響[25]。表層土壤水力學(xué)特性參數(shù)通過(guò)野外簡(jiǎn)易滲透試驗(yàn)確定，以表層土壤垂向滲透系數(shù)為例，其確定方法和過(guò)程如下：

以自制簡(jiǎn)易滲透試驗(yàn)裝置進(jìn)行野外試驗(yàn)（見(jiàn)圖5a）。試驗(yàn)裝置主要由支架（圖5a中未繪出）和標(biāo)記刻度的貯水器兩部分組成，下部設(shè)有模擬點(diǎn)降雨出水孔，根據(jù)水量平衡原理，自滲透試驗(yàn)開(kāi)始至有地表徑流產(chǎn)生的臨界狀態(tài)為止，表層土壤已達(dá)到或接近飽和狀態(tài)，表層土入滲量和土壤入滲能力達(dá)到平衡，此時(shí)降雨量全部入滲到土壤中。

以貯水器在地面投影為參照面積，根據(jù)貯水器直徑、水面下降高度及試驗(yàn)時(shí)間可計(jì)算出降雨強(qiáng)度。同一研究區(qū)研究結(jié)果表明，該地區(qū)表土飽和滲透系數(shù)數(shù)量級(jí)為10-6m·s-1[26]，表土初期（不飽和）滲透系數(shù)遠(yuǎn)大于飽和滲透系數(shù)。因此試驗(yàn)將雨量控制在13.27 cm3·min-1（相當(dāng)參照面積范圍內(nèi)降雨強(qiáng)度為4 mm·min-1）。試驗(yàn)結(jié)束后滲透輪廓如圖4b、4c所示，其中，水平面滲透輪廓近似為圓形，其直徑可用刻度尺沿著兩個(gè)垂直方向多次測(cè)量后取平均值確定；土壤垂直剖面滲透輪廓近似為弧形區(qū)域，其最大弦長(zhǎng)為地面滲透輪廓直徑，垂向最大滲透深度為地表至滲透輪廓線底部中心部位距離，每試驗(yàn)點(diǎn)測(cè)量3次取平均值，作為滲透深度。

圖5 簡(jiǎn)易滲透裝置示意圖及滲透輪廓Fig.5Sketch of simple infiltration device and infiltration profile

基于穩(wěn)定滲透理論和水量平衡原理建立計(jì)算模型如圖5a和圖6a所示，當(dāng)水降落到地面后，運(yùn)動(dòng)方式可以分解為水平方向的均勻擴(kuò)散及豎直方向的入滲兩部分。設(shè)到某時(shí)刻為止，降落到地面上的總水量為Vt，在單位時(shí)間Δt內(nèi)，水沿圓徑向擴(kuò)散的距離為Δrc，自圓心開(kāi)始沿徑向每增加Δrc對(duì)應(yīng)的圓環(huán)面積增量為ΔAc，則在水平方向上滲透的圓形區(qū)域?yàn)橐唤M環(huán)形（見(jiàn)圖6a），以ΔVi表示流入和流出第i個(gè)環(huán)形區(qū)域的水量差，則該差值為第i個(gè)環(huán)形區(qū)域內(nèi)的滲透量。根據(jù)水量平衡及穩(wěn)定滲透原理，可以建立入滲量、時(shí)間以及水平擴(kuò)散距離的函數(shù)關(guān)系式；在豎直方向上，單位時(shí)間內(nèi)的滲透深度變化Δh是入滲量ΔV、時(shí)間t以及土壤滲透能力的函數(shù)。圖6b為模型結(jié)構(gòu)下滲透過(guò)程中表土含水量變化示意圖。

各點(diǎn)滲透試驗(yàn)開(kāi)始時(shí)初始含水率和試驗(yàn)終止時(shí)飽和含水率分別設(shè)為θint和θsat，在穩(wěn)定滲透理論下，各滲透試驗(yàn)點(diǎn)初始和飽和含水量被設(shè)定為兩個(gè)不等的常數(shù)。滲透試驗(yàn)開(kāi)始后，表土層含水率在表土層內(nèi)（各滲透試驗(yàn)點(diǎn)最大深度h）隨滲透時(shí)間增加自上而下由初始含水率θint開(kāi)始逐漸增加至飽和含水率θsat（見(jiàn)圖6b，t3＞t2＞t1），可以根據(jù)上述水量和滲透關(guān)系建立式（6）～（9）所示的滲透計(jì)算方法。

①質(zhì)量（水量）平衡連續(xù)方程式

式中，t-時(shí)間因子（s）；rc-徑向空間因子（m）；f-滲透系數(shù)（m·s-1）；V-水體積，其變化量dV-單位時(shí)間內(nèi)滲透量（m3·s-1）；dAci-環(huán)形區(qū)域面積（m2）；rci-第i個(gè)環(huán)形區(qū)域內(nèi)徑（m）；λ-土壤有效孔隙率；h-到某一時(shí)刻為止的滲透深（m）；θ-土壤體積含水率（cm3·cm-3）；K-引入?yún)?shù)，用來(lái)表示與土壤含水率、飽和度及孔隙率關(guān)系，是與土壤實(shí)際滲透能力有關(guān)的參數(shù)。

通過(guò)上述試驗(yàn)和計(jì)算方法，可近似推算研究區(qū)域內(nèi)各種植被條件下表層土壤垂向滲透系數(shù)。

圖6 滲透模型Fig.6Infiltration model

兩個(gè)研究流域土壤垂直剖面模型（圖1，圖2）所示的第二層土壤垂向滲透系數(shù)、有效孔隙率等難以通過(guò)試驗(yàn)、實(shí)際調(diào)查以及DEM等方法確定的參數(shù)，首先為其賦予合理初值，通過(guò)模型海量計(jì)算考查數(shù)值模擬誤差的方法確定[26-27]。因?yàn)橥悈?shù)隨河網(wǎng)上各分布式小流域的不同而不同，土壤水力學(xué)特性參數(shù)呈現(xiàn)時(shí)空變動(dòng)特性。因此數(shù)值計(jì)算過(guò)程中，各計(jì)算參數(shù)隨計(jì)算進(jìn)程呈現(xiàn)動(dòng)態(tài)變化，因計(jì)算參數(shù)較多，相應(yīng)內(nèi)容在此略去。

1.3.6 蒸散發(fā)推求

蒸散發(fā)一般主要由坡面第一層土壤產(chǎn)生，因?yàn)榱罍狭饔蛴?jì)算區(qū)域面積較?。?.1 km2），計(jì)算區(qū)域主要為草地，故以草地蒸散發(fā)作為評(píng)價(jià)降雨-徑流計(jì)算時(shí)段內(nèi)由于蒸散發(fā)而導(dǎo)致的降雨損失（式4），利用觀測(cè)氣象數(shù)據(jù)以彭曼-蒙蒂斯（Penmen-Monteith）公式（10）計(jì)算六道溝流域草地蒸散發(fā)[28]。

式中，lET-土壤潛熱通量（MJ·m-2·h-1）；Δ-不同溫度下飽和水汽壓曲線上的斜率（kPa·℃-1）；Rn-凈輻射量（MJ·m-2·h-1）；γ-溫度表常數(shù)（kPa·℃-1）；G-土壤熱通量（MJ·m-2·h-1）；es-不同溫度下飽和水汽壓（kPa）；ea-實(shí)際水汽壓（kPa）；l-汽化潛熱（MJ·kg-1）；ET-蒸散發(fā)（mm·h-1）；cp-空氣比熱（1.0×10-3MJ·kg-1·℃-1）；ρ-空氣密度（kg·m-3）；ra-空氣阻抗（s·m-1）；rs-表面阻抗（s·m-1）。

基于實(shí)測(cè)氣象數(shù)據(jù)，率定式（10）中用于計(jì)算散發(fā)的各參數(shù)，進(jìn)而推算草地蒸散發(fā)[26,29]，相應(yīng)內(nèi)容在此略去。

在Bukuro流域，植被多樣且分布比較復(fù)雜，如在雨季，地面幾乎被綠色植物覆蓋，而植被蒸散發(fā)量也因生長(zhǎng)期內(nèi)氣象條件和植物本身生物條件而有所不同，蒸散發(fā)準(zhǔn)確計(jì)算難于實(shí)現(xiàn)，在降雨-徑流過(guò)程計(jì)算時(shí)，引入一個(gè)損失系數(shù)近似地代替蒸散發(fā)因子，該損失系數(shù)確定方法為首先賦予合理初值，通過(guò)模型海量計(jì)算考查數(shù)值模擬誤差方法確定。

2 結(jié)果與分析

2.1 模型驗(yàn)證

利用計(jì)算機(jī)軟件Fortran開(kāi)發(fā)計(jì)算程序（數(shù)值模型），通過(guò)對(duì)觀測(cè)地表徑流數(shù)值模擬檢驗(yàn)?zāi)Ｐ蛯?shí)用性和計(jì)算精度，六道溝流域和Bukuro流域?qū)崪y(cè)徑流過(guò)程數(shù)值模擬結(jié)果分別如圖7、8所示。

六道溝流域位于黃土高平原水蝕風(fēng)蝕交錯(cuò)區(qū)，年間降雨蒸發(fā)比為0.20～0.50，屬于典型半干旱區(qū)，又因計(jì)算區(qū)域尺度較小，平時(shí)河道內(nèi)無(wú)徑流存在，只有在雨季集中降雨發(fā)生期間才能產(chǎn)生短歷時(shí)徑流，且徑流存在時(shí)間較短，降雨結(jié)束后很快消失，通過(guò)觀測(cè)得到徑流數(shù)據(jù)相對(duì)較少。圖7所示為六道溝流域兩次降雨（圖7a為2005年8月15日，圖7b為2006年6月23日）-徑流事件的數(shù)值模擬結(jié)果。在Bukuro流域，即使在無(wú)雨期間，河道內(nèi)仍有基流存在，且有較長(zhǎng)時(shí)間（1999年5月31～7月8日）序列的實(shí)測(cè)徑流數(shù)據(jù)，故對(duì)該流域徑流過(guò)程進(jìn)行較長(zhǎng)時(shí)段的模擬（見(jiàn)圖8）。由圖7和8可知，數(shù)值模擬結(jié)果可較好再現(xiàn)觀測(cè)徑流產(chǎn)生過(guò)程。

圖7 六道溝流域降雨-徑流事件數(shù)值模擬結(jié)果Fig.7Simulation results of rainfall-runoff events for Liudaogou Basin

圖8Bukuro流域降雨-徑流數(shù)值模擬結(jié)果(1999年5月31～7月8日)Fig.8Simulation result of rainfall-runoff for Bukuro Basin

以皮爾森相關(guān)系數(shù)（Pearson correlation coefficient）評(píng)價(jià)觀測(cè)徑流序列和數(shù)值模擬序列間契合程度，公式如下：

式中，rp-皮爾森相關(guān)系數(shù)；n-計(jì)算次數(shù)或觀測(cè)數(shù)據(jù)個(gè)數(shù)；i-編號(hào)；Q-觀測(cè)流量（m3·s-1）；Qs-模擬值（m3·s-1）；Qˉ-觀測(cè)徑流序列；Qˉs-模擬值序列均值。

分別計(jì)算圖7所示的六道溝流域觀測(cè)徑流序列和對(duì)應(yīng)模擬序列之間的皮爾森相關(guān)系數(shù)，結(jié)果分別為0.98和0.99，圖8所示Bukuro流域觀測(cè)徑流和模擬序列皮爾森相關(guān)系數(shù)為0.91，說(shuō)明兩個(gè)研究流域觀測(cè)徑流與數(shù)值模擬結(jié)果序列之間相關(guān)性很高，針對(duì)兩個(gè)研究流域分別構(gòu)建的分布式降雨-徑流數(shù)值模型在降雨-徑流計(jì)算中具有實(shí)用性。

2.2 誤差分析

以式（12）作為觀測(cè)流量與數(shù)值模擬結(jié)果之間誤差判別基準(zhǔn)，該基準(zhǔn)要求誤差＜3%。

式中，Er-誤差；Qmax-計(jì)算時(shí)段內(nèi)最大觀測(cè)流量（m3·s-1）；其他因子同前述對(duì)應(yīng)相同。

對(duì)圖7模擬結(jié)果的峰值流量誤差進(jìn)行計(jì)算，結(jié)果分別為0.010和0.012，計(jì)算序列誤差分別為0.016和0.018；圖8所示模擬結(jié)果的峰值流量誤差為0.035，觀測(cè)徑流和模擬值兩序列之間誤差為0.015。對(duì)兩個(gè)研究流域其他時(shí)期降雨-徑流數(shù)值模擬結(jié)果隨機(jī)誤差分析表明，次降雨產(chǎn)生徑流以及長(zhǎng)時(shí)間序列（計(jì)算時(shí)段內(nèi)有多次降雨發(fā)生的情況）徑流模擬結(jié)果誤差均在誤差基準(zhǔn)允許范圍內(nèi)。而個(gè)別峰值流量與其模擬值之間誤差雖超出判斷基準(zhǔn)允許范圍，但由于計(jì)算次數(shù)（數(shù)據(jù)個(gè)數(shù)）n只是峰值產(chǎn)生時(shí)刻的單值，不符合誤差判斷基準(zhǔn)對(duì)數(shù)據(jù)序列數(shù)量要求。綜合數(shù)值模擬誤差分析結(jié)果可知，研究構(gòu)建的降雨-徑流數(shù)值模型具有很高模擬（計(jì)算）精度。

另一方面，六道溝流域降雨-徑流數(shù)值模型構(gòu)建過(guò)程中，以彭曼-蒙蒂斯公式計(jì)算草地蒸散發(fā)作為蒸散發(fā)（ET）因子用于降雨-徑流過(guò)程計(jì)算，該蒸散發(fā)因子通過(guò)觀測(cè)氣象數(shù)據(jù)計(jì)算得來(lái)，該方法在六道溝流域草地蒸散發(fā)計(jì)算中實(shí)用性已得到驗(yàn)證，模型具有較強(qiáng)物理性。而在Bukuro流域，因?yàn)橹脖欢鄻有詫?dǎo)致蒸散發(fā)難于準(zhǔn)確計(jì)算，引入損失系數(shù)近似評(píng)價(jià)蒸散發(fā)因子用于該流域降雨-徑流過(guò)程計(jì)算，并得到精度較高的模擬結(jié)果，但該損失系數(shù)在一定程度上降低了Bukuro流域降雨-徑流數(shù)值模型物理性。

3 結(jié)論

本研究中，基于六道溝流域和Bukuro流域的實(shí)際地形和水文地質(zhì)條件，利用運(yùn)動(dòng)波理論基礎(chǔ)方程式結(jié)合GIS-ArcMap分別構(gòu)建兩個(gè)研究流域分布式水文模型，對(duì)兩個(gè)流域降雨-徑流過(guò)程進(jìn)行數(shù)值模擬并對(duì)結(jié)果進(jìn)行分析，結(jié)論如下：

a.基于數(shù)值模擬檢驗(yàn)?zāi)Ｐ驮趦蓚€(gè)研究流域的適用性，結(jié)果表明，基于不同地形條件開(kāi)發(fā)的小尺度流域分布式水文模型，適用于各研究流域降雨-徑流過(guò)程計(jì)算。

b.通過(guò)對(duì)降雨-徑流過(guò)程數(shù)值模擬結(jié)果檢驗(yàn)和誤差分析，檢驗(yàn)兩個(gè)研究流域降雨-徑流過(guò)程數(shù)值模擬結(jié)果契合度和模型計(jì)算精度，結(jié)果表明，觀測(cè)流量序列和數(shù)值模擬結(jié)果之間相關(guān)性較高，皮爾森相關(guān)系數(shù)在0.90以上，模型計(jì)算精度可達(dá)到誤差＜3%。

[1]陳仁升,康爾泗,楊建平,等.水文模型研究綜述[J].中國(guó)沙漠, 2003,23(3)：222-229.

[2]金鑫,郝振純,張金良.水文模型研究進(jìn)展及發(fā)展方向[J].水土保持研究,2006,13(4)：197-199.

[3]閆紅飛,王船海,文鵬.分布式水文模型研究綜述[J].水電能源科學(xué),2008,26(6)：1-4.

[4]劉昌明,鄭紅星,王中根,等.流域水循環(huán)分布式模擬[M].鄭州：黃河水利出版社,2006.

[5]Abbott M B,Refsgaara J F.Distributed hydrological modelling [M].Boston：Kluwer Academic Publishers,1996.

[6]任立良.流域數(shù)字水文模型研究[J].河海大學(xué)學(xué)報(bào),2004,28 (4)：1-7.

[7]Jain M K,Kothyari U C,Ranga R K G.A GIS based distributed rainfall-runoff model[J].Journal of Hydrology,2004,299(1-2)：107-135.

[8]Karvonena T,Koivusaloa H,Jauhiainena M.A hydrological model for predicting runoff from different land use areas[J].Journal of Hydrology,1999,217(3-4)：253-265.

[9]唐麗莉,王守榮,顧駿強(qiáng).分布式水文模型DHSVM對(duì)蘭江流域徑流變化的模擬試驗(yàn)[J].熱帶氣象學(xué)報(bào),2008,24(2)：176-182.

[10]張利平,陳小鳳,張曉琳,等.VIC模型與SWAT模型在中小流域徑流模擬中的對(duì)比研究[J].長(zhǎng)江流域資源與環(huán)境,2009,18 (8)：745-752.

[11]張會(huì)蘭,李丹勛,王興奎.基于BPCC模型的徑流對(duì)氣候波動(dòng)和覆被變化的響應(yīng)[J].2010,50(3)：376-379.

[12]陳瑩,徐有鵬,陳興偉.長(zhǎng)江三角洲地區(qū)中小流域未來(lái)城鎮(zhèn)化的水文效應(yīng)[J].資源科學(xué),2011,33(1)：64-69.

[13]楊?yuàn)檴?徐征和,孔珂,等.基于SWAT模型的臥虎山水庫(kù)流域徑流模擬[J].中國(guó)農(nóng)村水利水電,2013(5)：11-14,19.

[14]李致家,王秀慶,呂雁翔,等.TOPKAPI模型的應(yīng)用及與新安江模型的比較研究[J].水力發(fā)電,2013,39(11)：6-10.

[15]小流域[DB/OL].http：//baike.baidu.com/link?url=IN0NZhP8QFJ5S3Un4Ddvj_DdAVQq44EvNIm1iK_Ey1uz7MseF7YHJsCoVG 7m_PXCelQopEZuqeFdne09ENsY_#1.2015-09-15.

[16]黃凱,劉永,郭懷成,等.小流域水環(huán)境規(guī)劃方法框架及應(yīng)用[J].環(huán)境科學(xué)研究,2006,19(5)：136-141.

[17]朱麗,秦富倉(cāng),姚云峰,等.SWAT模型靈敏性分析模塊在中尺度流域的應(yīng)用—以密云縣紅門川流域?yàn)槔齕J].水土保持研究, 2001,18(1)：161-166.

[18]郭建中,胡和平,翁文斌.中小流域防洪規(guī)劃決策支持系統(tǒng)—Ⅰ系統(tǒng)研究[J].水科學(xué)進(jìn)展,2001,12(2)：222-226.

[19]胡和平,郭建中,翁文斌.中小流域防洪規(guī)劃決策支持系統(tǒng)—Ⅱ個(gè)例分析[J].水科學(xué)進(jìn)展,2001,12(2)：227-231.

[20]Yomoto A,Islam M N.Kinematic analysis of flood runoff for a small-scale upland field[J].Journal of Hydrology,1992,137 (1-4)：311-326.

[21]Chua Lloyd H C,Wong Tommy S W,Sriramula L K.Comparison between kinematic wave and artificial neural network models in event based runoff simulation for an overland plane[J].Journal of Hydrology,2008,357(3-4)：337-348.

[22]Cabral M C,Garrote L,Bras R L,et al.A kinematic model of infiltration and runoff generation in layered and sloped soils[J].Advances in Water Resources,1992,15(5)：311-324.

[23]Sarkar R,Dutta S.Field investigation and modeling of rapid subsurface storm flow through preferential pathways in a vegetated hillslope of northeast India[J].Journal of Hydrologic Engineering, 2012,17(2)：333-341.

[24]Huang J B,Li J,Yasuda H,et al.Rainfall characteristics of the Liudaogou Catchment on the northern Loess Plateau of China[J]. Journal of Rangeland Science,2013,3(3)：252-264.

[25]Huang J B,Wen J W,Wang B,et al.Numerical analysis of the combined rainfall-runoff process and snowmelt for the Alun River Basin,Heilongjiang,China[J].Environmental Earth Sciences, DOI：10.1007/s12665-015-4694-y.

[26]Huang J B.Study on runoff and water balance in the northern Loess Plateau,China[D].Tottori：Tottori University,Japan,2010.

[27]Huang J B,Hinokidani O,Yasuda H,et al.Study on characteristics of the surface flow of the upstream region of Loess Plateau[J]. Annual Journal of Hydraulic Engineering,JSCE,2008,52：1-6.

[28]Kimura R,Fan J,Zhang X C,et al.Evapotranspiration over the grassland field in the Liudaogou basin of the Loess Plateau,China [J].Acta Oecologica,2005,29(1)：45-53.

Numerical simulation of rainfall-runoff for small-scale basin of two dif- ferent topographic conditions/

ZHOU Fanglu1,2,HUANG Jinbai3,WANG Bin1,HINOKIDANI
Osamu4
(1.School of Water Conservancy and Architecture,Northeast Agricultural University,Harbin 150030,China;2.Daqing Songnen Project Management Office,Daqing Heilongjiang,163300,China; 3.School of Hydraulic Energy and Power Engineering,Yangzhou University,Yangzhou Jiangsu 225009,China;4.Faculty of Engineering of the Graduate School of Tottori University,Tottori 680-8552, Japan)

The Liudaogou Basin which is located at hilly-gully area of the northern Loess Plateau and Bukuro Basin which is located at mountain area of the Tottori County of Japan were adopted as the study locations.Modeling of soil vertical profile of the study locations was carried out;the numerical model for rainfall-runoff was developed based on the equations of kinematic wave theory combined with GIS-ArcMap; Applicability of the model was validated through numerical simulation of the observed runoff process.The results indicated that the simulation results were with relatively high precision,the error was less than 3%. And Pearson correlation coefficient between the curve of the observed flow discharge and the calculated one was more than 0.90.A scientific method for surface runoff estimation was therefore provided for the study locations.

basin;topography;distributed hydrological model;rainfall-runoff;numerical simulation

TV121.1；P333.9

A

1005-9369（2015）09-0093-09

時(shí)間2015-9-23 10：43：27[URL]http：//www.cnki.net/kcms/detail/23.1391.S.20150923.1043.026.html

周方錄,黃金柏,王斌,等.兩種地形條件下小尺度流域降雨-徑流數(shù)值模擬[J].東北農(nóng)業(yè)大學(xué)學(xué)報(bào),2015,46(9):93-101.

Zhou Fanglu,Huang Jinbai,Wang Bin,et al.Numerical simulation of rainfall-runoff for small-scale basin of two different topographic conditions[J].Journal of Northeast Agricultural University,2015,46(9):93-101.(in Chinese with English abstract)

2015-03-17

國(guó)家自然科學(xué)基金項(xiàng)目（41271046）；揚(yáng)州大學(xué)新世紀(jì)人才工程項(xiàng)目-優(yōu)秀青年教師（5015-137050209）；揚(yáng)州市“綠揚(yáng)金鳳計(jì)劃”優(yōu)秀博士人才項(xiàng)目（yzlyjfjh2013YB105）

周方錄（1972-），男，高級(jí)工程師，碩士，研究方向?yàn)榱饔蛩呐c水資源、水資源優(yōu)化利用。E-mail：zhoufanglu@163.com

*通訊作者：黃金柏，副教授，研究方向?yàn)樗倪^(guò)程數(shù)值模擬、數(shù)字流域、河流工學(xué)。E-mail：huangjinbai@aliyun.com