趙建飛，梁忠民，劉金濤，2，李彬權(quán)，段雅楠

(1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；2. 河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098)

中國大部分山丘區(qū)是洪水突發(fā)、易發(fā)區(qū)，由于水文氣象條件及特殊地形的綜合作用，山洪災(zāi)害頻繁，是中國防汛工作的難點和薄弱環(huán)節(jié)。以數(shù)據(jù)為導(dǎo)向的經(jīng)驗或概念性模型難以滿足山丘區(qū)洪水預(yù)報的需求[1]，而分布式水文模型因其可以實現(xiàn)任一位置/斷面的洪水預(yù)報而成為未來趨勢，同時這也是水文模型自身發(fā)展的必然方向。因此，研究面向山丘區(qū)的分布式水文模型具有重要的意義。

最新產(chǎn)匯流理論研究表明，壤中流是構(gòu)建山丘區(qū)分布式水文模型的重點和難點。壤中流不僅直接參與形成徑流過程，也可通過影響土壤水分的時空分布進(jìn)而主導(dǎo)流域的水文響應(yīng)，是大多數(shù)山坡徑流形成的主控因素[2-4]。壤中流存在飽和與非飽和2種狀態(tài)，包含垂向和側(cè)向2種流向，并表現(xiàn)出高度的非線性和空間異質(zhì)性?，F(xiàn)有分布式水文模型對壤中流的模擬方法可分為3類[5]，即Richards模型、山坡蓄泄模型和運動波模型。Richards模型的優(yōu)點在于能精細(xì)刻畫壤中流運動的微觀規(guī)律，但存在算力要求高、計算耗時長的問題[6]，難以滿足快速預(yù)報的需求；山坡蓄泄模型為蓄滿產(chǎn)流模式，理論上僅適用于恒定雨強、規(guī)則山坡等理想情況，對實際復(fù)雜流域需進(jìn)行較多概化[7]，一定程度上喪失了其理論優(yōu)勢；運動波模型是由Lighthill和Whitham[8]在研究河道洪水演進(jìn)問題時提出的，Henderson和Wooding[9]、Beven[10]分別將其引入飽和側(cè)向壤中流與非飽和壤中流的模擬。Cabral等[11]進(jìn)一步發(fā)展了該模型，提出了考慮坡度、土壤分層和各向異性影響的Cabral產(chǎn)流模型；在此基礎(chǔ)上，Ivanov等[12]考慮毛管力和變動地下水位作用，耦合基于非結(jié)構(gòu)三角網(wǎng)格的匯流模型，提出tRIBS模型，隨后得到較多應(yīng)用[13]。研究認(rèn)為，運動波模型能充分反映產(chǎn)流的物理機制，可滿足山丘區(qū)洪水預(yù)報的實時性要求，代表壤中流模擬的一個重要發(fā)展方向[14]，但目前的研究忽略了土壤-風(fēng)化基巖界面的相對不透水層作用，無法同時模擬土壤層間的壤中流、土壤-風(fēng)化基巖界面(土巖界面)的壤中流和地下徑流3種成分，模型在理論上存在一定不足。

為此，本文根據(jù)山坡水文試驗的新認(rèn)識，考慮土壤-風(fēng)化基巖界面在產(chǎn)流中的作用，提出變動產(chǎn)流層概念。在此基礎(chǔ)上，根據(jù)運動波模型理論，推導(dǎo)超滲/飽和地面徑流和壤中流的計算公式，進(jìn)一步耦合植被截留、蒸散發(fā)、二維地表和一維河道匯流模型，構(gòu)建了一個新的分布式流域水文模型——山丘區(qū)變動產(chǎn)流層分布式水文模型(VRGL)，并在典型濕潤山丘區(qū)的屯溪流域進(jìn)行示例應(yīng)用。

1 變動產(chǎn)流層

1.1 已有認(rèn)識

從山坡產(chǎn)流的垂向來看，上邊界為植被冠層的頂端，下邊界為一理想的相對或絕對不透水層。盡管不同山坡的表層關(guān)鍵帶結(jié)構(gòu)存在一定區(qū)別[3，15-16]，但最新研究表明，山坡產(chǎn)流的下邊界一般對應(yīng)不受風(fēng)化作用影響的新鮮基巖界面[15，17]；新鮮基巖以上為風(fēng)化基巖，風(fēng)化基巖以上覆蓋土壤[18]，兩者之間存在一個相對不透水的土壤-風(fēng)化基巖界面；地表至土壤-風(fēng)化基巖界面之間的垂向深度為土壤厚度(圖1)。

圖1 山坡表層關(guān)鍵帶結(jié)構(gòu)示意Fig.1 Schematic diagram of hillslope critical zone structure

山坡水文試驗的最新研究表明，壤中流也會發(fā)生在土壤-風(fēng)化基巖界面[3]。McDonnell[19]在MaiMai實驗流域觀測并系統(tǒng)描述了土壤-風(fēng)化基巖界面上的臨時飽和區(qū)和對應(yīng)的側(cè)向壤中流；Fu等[20]在中國南方的徑流小區(qū)試驗表明，當(dāng)滿足初始土壤含水量和降雨閾值條件時，土壤-風(fēng)化基巖界面會出現(xiàn)壤中流；Zhang等[21]在四川龍溪河實驗流域發(fā)現(xiàn)，當(dāng)降雨超過特定閾值后，土壤-風(fēng)化基巖界面的臨時飽和區(qū)開始向上發(fā)展，壤中流出流驟增；Han等[3]在浙江姜灣流域的試驗則進(jìn)一步提供了新的認(rèn)識，研究區(qū)山坡存在土壤-風(fēng)化基巖界面和土壤層間2種壤中流水流路徑，前者是研究區(qū)壤中流水流路徑的主要形式，后者主要形成于較大雨強期間。

土壤-風(fēng)化基巖界面以下是風(fēng)化基巖層，包含一定的快速輸水通道，是地下徑流的形成區(qū)及河道基流的重要來源，對山坡產(chǎn)流過程也有重要貢獻(xiàn)[4]。Rempe和Dietrich[15]指出在Elder Creek流域的一個研究區(qū)內(nèi)，風(fēng)化基巖層可存儲全年降水量的27%，是山丘區(qū)的一個重要蓄水庫。Bouvier等[4]通過Valescure流域內(nèi)徑流小區(qū)的觀測和試驗發(fā)現(xiàn)，上坡位風(fēng)化基巖層中的水分可在重力和壓力作用下進(jìn)入下坡位土壤層，參與土壤-風(fēng)化基巖界面臨時飽和層的形成。

1.2 概念提出

由上可知，近年來山坡水文試驗特別是對山坡表層關(guān)鍵帶結(jié)構(gòu)的研究表明，土壤-風(fēng)化基巖界面和風(fēng)化基巖層在山坡產(chǎn)流中發(fā)揮著重要作用，但已有的產(chǎn)流模型尚無法充分反映這一特點。為此，本文提出變動產(chǎn)流層概念(認(rèn)知模型)，描述山坡產(chǎn)流機制，其徑流形成及與流量過程的對應(yīng)關(guān)系如圖2(θ為土壤含水量，n為地表垂向深度)所示。

圖2 變動產(chǎn)流層概念示意Fig.2 Schematic diagram of the variable runoff generation layer concept

對1個具體的計算單元，降雨首先進(jìn)入植被層，扣除截留等損失后以穿透雨的形式進(jìn)入地表。穿透雨進(jìn)入地表后，根據(jù)地表下滲強度判斷是否產(chǎn)生超滲地面徑流，地表下滲強度由地表飽和水力傳導(dǎo)度和土壤含水量共同決定。下滲水量補充土壤含水量，形成的濕潤鋒不斷向下運動。一般而言，土壤飽和水力傳導(dǎo)度沿垂向呈衰減特性，當(dāng)濕潤鋒運動到某一深度后，若下滲強度大于該深度的飽和水力傳導(dǎo)度，則在該深度上會有水分堆積，使得該深度的土壤達(dá)到飽和，并向上/向下發(fā)展成一個臨時的飽和土壤層(但該層以上、以下的土壤可能尚未飽和)；同時，在地形坡度作用下，在該土壤層上將形成側(cè)向壤中流。該臨時飽和土壤層的位置受下滲強度和前期土壤含水量的影響而變動，本文稱為“土壤層間變動產(chǎn)流層”，其水位稱為“土壤層間水位”。

隨著土壤水分的運動，對土壤-風(fēng)化基巖界面的補給增加，在該界面上同樣會有水分堆積，使得該界面處的土壤逐漸達(dá)到飽和，并向上發(fā)展，同時產(chǎn)生側(cè)向壤中流，本文將該臨時飽和土壤層稱為“土巖界面變動產(chǎn)流層”，其水位稱為“土巖界面水位”。土壤-風(fēng)化基巖界面上的水分可通過下滲對風(fēng)化基巖層進(jìn)行補給，并參與形成地下徑流；風(fēng)化基巖層的水分補給一方面來源于本單元土壤層的下滲水量，也可來源于鄰近單元的同層側(cè)向補給；當(dāng)整個風(fēng)化基巖層達(dá)到飽和后，多余水量會透過土壤-風(fēng)化基巖界面向土壤層反向補給。隨著土壤層間變動產(chǎn)流層和土巖界面變動產(chǎn)流層的繼續(xù)發(fā)展，當(dāng)土壤層間水位或土巖界面水位到達(dá)地表后，將形成飽和地面徑流。

由上可知，通過引入土壤-風(fēng)化基巖界面，變動產(chǎn)流層概念將對壤中流的描述從單一的土壤層間拓展到土巖界面，體現(xiàn)了山坡水文試驗的新認(rèn)識。同時，在變動產(chǎn)流層框架下，借助對壤中流形成及發(fā)展過程的精細(xì)刻畫，可同時描述飽和地面徑流、超滲地面徑流、壤中流和地下徑流4種徑流成分。對1個具體的計算單元，在降雨過程中其產(chǎn)流模式可根據(jù)土壤濕潤鋒的運動自動確定并動態(tài)轉(zhuǎn)換，避免對產(chǎn)流模式進(jìn)行蓄滿或超滲的人為假定。

2 變動產(chǎn)流層分布式水文模型結(jié)構(gòu)

2.1 模型概述

變動產(chǎn)流層分布式水文模型(VRGL)是一個充分反映土壤-風(fēng)化基巖界面作用的分布式水文模型，通過對壤中流形成及發(fā)展過程的精細(xì)刻畫，實現(xiàn)了蓄滿-超滲產(chǎn)流及其轉(zhuǎn)換機制的統(tǒng)一描述。VRGL以DEM網(wǎng)格為計算單元，在計算單元上構(gòu)建產(chǎn)流模塊，模擬植被截留層、土壤層和風(fēng)化基巖層的產(chǎn)流過程并考慮單元間的水量交換，采用二維擴散波和一維擴散波分別進(jìn)行地表匯流和河道匯流，模型結(jié)構(gòu)如圖3所示。

圖3 VRGL模型結(jié)構(gòu)示意Fig.3 Schematic diagram of the VRGL model structure

2.2 產(chǎn)流模型

2.2.1 植被截留

采用Bucket Model模擬植被截留[22]，其水量平衡方程為

(1)

式中：Mv為截留量，mm；t為時間，h；P為降雨強度，mm/h；Ev為植被截留層的實際蒸散發(fā)率，mm/h；Pv為穿透雨強，mm/h。

2.2.2 土壤層產(chǎn)流

為便于產(chǎn)流計算，根據(jù)穩(wěn)態(tài)土壤水分剖面，將土壤層的θ—n域動態(tài)劃分為域Ⅰ和域Ⅱ 2個部分(圖4，θs為土壤飽和含水量；θr為土壤殘余含水量；Nsd為土壤-風(fēng)化基巖界面深度；Nwt為土巖界面水位深度；Nst為土壤層間水位深度)。域Ⅰ和域Ⅱ上邊界均為地表，域Ⅰ下邊界為土壤-風(fēng)化基巖界面，域Ⅱ下邊界為土巖界面水位。初始時刻，土巖界面水位與土壤-風(fēng)化基巖界面重合(圖4(a))；隨著水分的積累，在土壤-風(fēng)化基巖界面上逐漸形成“土巖界面變動產(chǎn)流層”，土巖界面水位開始抬升(圖4(b)和圖4(c))；隨著降雨的持續(xù)，土巖界面水位不斷抬升直至達(dá)到地表，域Ⅱ消失，整個土層歸為域Ⅰ并達(dá)到飽和(圖4(d))。

圖4 土壤層θ—n域劃分示意Fig.4 Schematic diagram of the soil layer moisture content-depth domain partition

產(chǎn)流計算時，對每個計算時段，分別進(jìn)行域Ⅰ和Ⅱ的計算，兩者共同構(gòu)成整個土壤層的產(chǎn)流量，徑流成分可能包括超滲/飽和地面徑流和壤中流。

域Ⅰ和域Ⅱ接收穿透雨入滲和上游網(wǎng)格側(cè)向壤中流的補給，對域Ⅰ，當(dāng)本計算單元的風(fēng)化基巖層含水量達(dá)到最大蓄量時，多余水量也會補給該域，但若未達(dá)到最大蓄量時，則域Ⅰ的水分將下滲補給風(fēng)化基巖層。域Ⅰ和域Ⅱ形成的側(cè)向壤中流補給下游網(wǎng)格(作為下游網(wǎng)格的側(cè)向入流)。對土壤的蒸發(fā)，本次假定：先由域Ⅱ的水分供給蒸發(fā)，消耗殆盡后，再進(jìn)行域Ⅰ蒸發(fā)計算。當(dāng)穿透雨強度大于域Ⅰ與域Ⅱ的下滲強度之和時，將產(chǎn)生超滲地面徑流。域Ⅰ和域Ⅱ均能形成飽和地面徑流，但機制不同：對域Ⅰ，只有當(dāng)土巖界面水位抬升至地表時(即整個土層達(dá)到飽和，圖4(d))，才能形成飽和地面徑流；對域Ⅱ，當(dāng)土壤層間水位抬升至地表時，濕潤鋒以上直至地表均達(dá)到飽和(但此時濕潤鋒以下的土壤可能并未達(dá)到飽和)，將形成飽和地面徑流(圖4(c))。具體計算過程如下：

(1) 域Ⅰ的產(chǎn)流計算

域Ⅰ的水量平衡方程為

(2)

式中：Mg為域Ⅰ的土壤含水總量，mm；t為時間，h；Ig為域Ⅰ的地表下滲強度，mm/h；qg，in、qg，out分別為域Ⅰ上游網(wǎng)格側(cè)向壤中流總?cè)肓鲝姸群捅揪W(wǎng)格域Ⅰ側(cè)向壤中流出流強度，mm3/h；A為計算單元面積，mm2；qwr是風(fēng)化基巖層達(dá)到最大蓄量后對域Ⅰ的補給強度，mm/h；qsr為域Ⅰ對風(fēng)化基巖層的下滲補給強度，mm/h；qd，g為域Ⅰ形成的飽和地面徑流，mm/h；Eg為域Ⅰ的蒸發(fā)強度，mm/h。

Mg通過對域Ⅰ的土壤水分剖面沿n方向積分得到，公式為

(3)

式中：Rg為域Ⅰ的等價恒定雨強[12]，mm/h；θ(Rg，n)為域Ⅰ的未飽和土壤水分剖面分布，mm3/mm3，根據(jù)Cabral產(chǎn)流模型[11]的基本假設(shè)，則有：

(4)

式中：K0n為n方向的地表飽和水力傳導(dǎo)度，mm/h；ε為土壤孔隙分布指數(shù)；f為飽和水力傳導(dǎo)度隨深度衰減系數(shù)，mm-1。

隨著域Ⅰ土壤水分的累積，土壤-風(fēng)化基巖界面處的土壤含水量逐漸變大，當(dāng)該深度的土壤含水量達(dá)到飽和臨界狀態(tài)時(即θn=Nsd=θs)，將形成土巖界面變動產(chǎn)流層，根據(jù)式(3)和式(4)可推導(dǎo)出對應(yīng)的土壤含水總量臨界值的計算公式為

(5)

式中：Mg，*為土壤含水總量臨界值，mm。

按照公式(2) 所述的水量平衡方程進(jìn)行產(chǎn)流計算時，根據(jù)是否形成土巖界面變動產(chǎn)流層(Mg是否超過Mg，*)，分成2種情況，具體如下：

① 未形成土巖界面變動產(chǎn)流層，即Mg

等價恒定雨強公式為

(6)

地表下滲強度公式為

Ig=Rgcosα

(7)

壤中流強度公式為

qg=NsdRg(ar-1)sinαL

(8)

式中：qg為域Ⅰ的壤中流強度，mm3/h ；ar為各向異性比率；α為地表坡度；L為壤中流過流寬度，mm。

由于此時土巖界面水位未抬升到地表，域Ⅰ不產(chǎn)生飽和地面徑流，即qd，g=0。

域Ⅰ對風(fēng)化基巖層的下滲強度見后文公式(31)和公式(32)；風(fēng)化基巖層對域Ⅰ的補給強度(qwr)見后文公式(33)。

公式(2)中的蒸發(fā)強度公式為

Eg=min[max(Ep-Ev-Eu， 0)， (Mg-Mg，i)/dt]

(9)

式中：Ep為蒸散發(fā)能力，mm/h；Eu為域Ⅱ的蒸發(fā)強度，mm/h，見后文公式(22)；Mg，i為Mg的最小取值，mm，可通過在公式(4)中取Rg的最小值并代入公式(3)計算，本文Rg的最小值取為0.01 mm/h。

② 已形成土巖界面變動產(chǎn)流層，即Mg≥Mg，*，可推導(dǎo)出如下計算公式。

此時，土巖界面水位深度的土壤含水量達(dá)到飽和，即θ(Rg，Nwt)=θs，根據(jù)公式(3) 和(4) ，可推導(dǎo)出域Ⅰ的土壤含水總量公式為

(10)

公式(10)為Nwt的非線性方程，在Mg和其余參數(shù)已知時，可根據(jù)該式反向求解，公式為

(11)

式中：c1=c4(Mg-θsNsd)+c3；c2=c3c4；c3=θs-θr；c4=-f/ε；LambertW函數(shù)為f(w)=w·exp(w)函數(shù)的反函數(shù)。

等價恒定雨強公式為

Rg=K0nexp(-fNwt)

(12)

地表下滲強度公式為

(13)

壤中流強度公式為

(14)

飽和地面徑流強度公式為

qd，g=max[(Mg-θsNsd)/dt， 0]

(15)

同樣，此時域Ⅰ對風(fēng)化基巖層的下滲強度見后文公式(31)和公式(32)；風(fēng)化基巖層對域Ⅰ補給強度見后文公式(33);蒸發(fā)強度計算同公式(9)。

(2) 域Ⅱ的產(chǎn)流計算

域Ⅱ的水量平衡公式為

(16)

式中：Mu為域Ⅱ的土壤含水總量，mm；Iu為域Ⅱ的地表下滲強度，mm/h；qu，in、qu，out分別為域Ⅱ上游網(wǎng)格側(cè)向壤中流總?cè)肓鲝姸群捅揪W(wǎng)格域Ⅱ側(cè)向壤中流出流強度，mm3/h；qd，u為飽和地面徑流，mm/h；域Ⅱ的蒸發(fā)強度記為Eu，mm/h。

當(dāng)穿透雨強滿足域Ⅰ的下滲補給時，剩余部分開始通過下滲進(jìn)入域Ⅱ，形成濕潤鋒并不斷向下發(fā)展，記濕潤鋒深度為Nf(單位為mm)。當(dāng)濕潤鋒到達(dá)某一深度時，等價恒定雨強超過該深度的飽和水力傳導(dǎo)度，此時土壤層間變動產(chǎn)流層形成，該深度稱為臨界深度，根據(jù)Cabral產(chǎn)流模型[11]，臨界深度的計算公式為

(17)

式中：Nf，*為濕潤鋒臨界深度，mm；Ru為域Ⅱ的等價恒定雨強，mm/h。

按照公式(16)進(jìn)行產(chǎn)流計算時，根據(jù)土壤層間變動產(chǎn)流層是否形成(Nf是否超過Nf，*)，分成2種情況，具體如下：

① 未形成土壤層間變動產(chǎn)流層，即Nf

等價恒定雨強公式為

(18)

濕潤鋒運動公式為

(19)

土壤層間變動產(chǎn)流層未形成時，可認(rèn)為土壤層間水位深度與Nf相同，即Nst=Nf。

地表下滲強度公式為

Iu=(Ru-Rg)cosα

(20)

壤中流強度公式為

qu=Nf(Ru-Rg)(ar-1)sinαL

(21)

式中：qu為域Ⅱ的側(cè)向壤中流強度，mm3/h。

由于此時土壤層間水位未抬升到地表，域Ⅱ不產(chǎn)生飽和地面徑流，即qd，u=0。

蒸發(fā)強度公式為

Eu=min[max(Ep-Ev， 0)，Mu/dt]

(22)

② 已形成土壤層間變動產(chǎn)流層，即Nf≥Nf，*，可推導(dǎo)出如下計算公式。

等價恒定雨強公式為

Ru=K0nexp(-fNst)

(23)

濕潤鋒運動公式為

(24)

土壤層間水位深度公式為

(25)

地表下滲強度公式為

(26)

壤中流強度公式為

(27)

飽和地面徑流強度公式為

(28)

蒸發(fā)強度計算同公式(22)。

從上可以看出，整個土壤層的產(chǎn)流計算包括域Ⅰ和域Ⅱ 2個部分，形成的徑流成分可能包括超滲/飽和地面徑流以及壤中流。對于域Ⅰ，根據(jù)土巖界面產(chǎn)流層是否形成，使用式(7)或式(13)計算Ig，使用式(8)或式(14)計其qg，當(dāng)土巖界面水位抬升到地表后，域Ⅰ形成飽和地面徑流對應(yīng)為公式(15)。對于域Ⅱ，根據(jù)土壤層間變動產(chǎn)流層是否形成，使用式(20)或式(26)計算其Iu，使用式(21)或式(27)計算其qu，當(dāng)土壤層間水位抬升到地表后，域Ⅱ形成飽和地面徑流對應(yīng)為公式(28)。對于整個土壤層，其飽和地面徑流(qd)為域Ⅰ和域Ⅱ形成的飽和地面徑流之和，即qd=qd，u+qd，g，單位為mm/h；當(dāng)穿透雨強超過域Ⅰ和域Ⅱ的地表下滲強度之和時，才會形成超滲地面徑流(qh)，即qh=max(Pv-Iu-Ig，0)，單位為mm/h。

2.2.3 風(fēng)化基巖層產(chǎn)流

將風(fēng)化基巖層的蓄泄過程概化為非線性水庫，以計算地下徑流，其水量平衡方程為

(29)

式中：Mw為非線性水庫的蓄量，mm；qw，in、qw，out分別為上游網(wǎng)格流入的地下徑流總強度和本網(wǎng)格流向下游網(wǎng)格的地下徑流強度，mm/h。

假設(shè)新鮮基巖界面坡度與地表坡度一致，則可使用考慮蓄滿程度的非線性運動波方程計算非線性水庫地下徑流出流強度，為

qw=kw(Mw/Mw，max)btanα

(30)

式中：qw為非線性水庫地下徑流出流強度，mm/h；kw為風(fēng)化基巖層出流系數(shù)，mm/h；Mw，max為非線性水庫的最大蓄量，mm；b為形狀參數(shù)?？梢圆捎檬?30)計算上游網(wǎng)格匯入本網(wǎng)格的地下徑流強度，其和即為qw，in；而本網(wǎng)格按此公式計算的地下徑流出流強度，即為qw，out。

土壤層域Ⅰ對風(fēng)化基巖層的下滲補給強度由風(fēng)化基巖層的剩余蓄水量和土壤層底部Nsd處的n方向水流強度決定，為

qsr=max[qn，sd， (Mw，max-Mw)/dt]

(31)

式中：qn，sd為土壤層底部Nsd處的n方向水流強度，mm/h。根據(jù)是否形成土巖界面變動產(chǎn)流層，qn，sd采用不同公式：

(32)

當(dāng)風(fēng)化基巖層非線性水庫蓄滿時，超出的水量補給上層土壤，公式為

qwr=max(Mw-Mw，max， 0)/dt

(33)

2.3 匯流模型

分別采用一維擴散波[22]、二維擴散波[23]進(jìn)行河道匯流、地表匯流計算，數(shù)值求解方法為有限體積法，時間推進(jìn)使用三階TVD Runge Kutta法[24]。

3 模型應(yīng)用

3.1 研究區(qū)域及資料

本文選取屯溪流域作為研究區(qū)域。屯溪流域位于新安江水系的源頭，是典型的濕潤山丘區(qū)流域，面積約為2 670 km2，多年平均降水量約為1 750 mm。流域?qū)贀P子地層區(qū)的江南地層分區(qū)，褶皺、斷裂等構(gòu)造發(fā)育，物理風(fēng)化強烈，主要地貌類型有山地、高低丘陵和山間盆地，地形周高中低、山高坡陡。流域土壤分布有垂直地帶性，典型土壤剖面分層明顯，森林覆蓋率高，壤中流發(fā)育。暴雨類型主要分為鋒面型、低壓型、臺風(fēng)外圍型和對流單體型，受特殊地形的強對流作用影響，降水強度大，洪水陡漲陡落。

使用的數(shù)據(jù)可分為站點數(shù)據(jù)和柵格數(shù)據(jù)2類。站點數(shù)據(jù)主要為水文氣象資料，包含屯溪水文站流量數(shù)據(jù)、21個雨量站的降水?dāng)?shù)據(jù)、日尺度的蒸發(fā)數(shù)據(jù)和河道斷面數(shù)據(jù)，數(shù)據(jù)來源為水文年鑒、安徽省水信息網(wǎng)和中國地面氣候資料日值數(shù)據(jù)集(V3.0)，資料年限為2010—2019年。場次洪水期間的降水流量資料主要來自降摘表與洪摘表，同時也收集了場次洪水前6個月的日尺度降水流量資料用于模型預(yù)熱。屯溪流域及站點分布如圖5所示。柵格數(shù)據(jù)主要為DEM、土壤類型、植被覆蓋等，數(shù)據(jù)及來源見表1。

圖5 屯溪流域及站點分布Fig.5 Tunxi watershed and its observation station distribution

表1 屯溪流域柵格數(shù)據(jù)及來源表

3.2 結(jié)果及分析

根據(jù)上述資料構(gòu)建模型，計算單元采用DEM網(wǎng)格，分辨率為1 000 m×1 000 m，網(wǎng)格數(shù)為5 950。假定網(wǎng)格內(nèi)降水和地形地貌等下墊面特征分布一致，不同網(wǎng)格間存在區(qū)別。降水、蒸發(fā)、流量數(shù)據(jù)的時間分辨率均統(tǒng)一處理為1 h，降水?dāng)?shù)據(jù)通過普通克里金法插值到計算網(wǎng)格，蒸發(fā)數(shù)據(jù)使用實測資料作為流域統(tǒng)一值。通過投影和重采樣計算，得到每個網(wǎng)格的植被、土壤、高程等數(shù)據(jù)。模型使用Python語言編寫，采用模塊化設(shè)計，有較好的跨平臺性和可拓展性。

選取了24場洪水進(jìn)行模擬研究，其中19場用于率定，5場用于驗證；評價指標(biāo)選取洪峰相對誤差、洪量相對誤差、納什確定性系數(shù)(ENS)和Kling-Gupta 效率系數(shù)(EKG)[25]共4個指標(biāo)。本次采用人工率定的方法確定參數(shù)，主要產(chǎn)流參數(shù)共8個，需要率定的參數(shù)有3個，其中飽和水力傳導(dǎo)度隨深度衰減系數(shù)根據(jù)場次干旱且雨強較大的場次洪水初步確定、土壤孔隙分布指數(shù)根據(jù)土壤參數(shù)數(shù)據(jù)和土壤類型數(shù)據(jù)估算初值、土壤各向異性比率參照已有研究取固定值，其余參數(shù)均可以根據(jù)柵格數(shù)據(jù)計算得到；匯流參數(shù)包括地表曼寧系數(shù)、河道曼寧系數(shù)和堰流公式流量系數(shù)等主要參數(shù)，本次采用相關(guān)文獻(xiàn)的推薦方法確定[22]。模型的主要參數(shù)及取值見表2。

表2 VRGL主要參數(shù)及取值表

模型精度評定結(jié)果如圖6所示。從圖6(a)可以看出，不管是率定期還是驗證期，洪峰流量和洪量的相對誤差均在±20%范圍內(nèi)，其中，洪峰流量相對誤差的絕對值平均為5.5%(率定期為5.2%、驗證期為6.6%)，洪量相對誤差的絕對值平均為7.4%(率定期為7.3%、驗證期為7.7%)。在洪水過程的評定方面，從圖6(b)可以看出，率定期和驗證期ENS均值分別為0.83和0.85，EKG均值分別為0.83和0.82，ENS和EKG2個指標(biāo)的評價結(jié)果較為一致；另外，除20140513場次ENS為0.56外，其余場次ENS均大于0.7，率定期和驗證期的ENS最大值可達(dá)0.95和0.91。上述結(jié)果表明，VRGL對洪峰流量、洪量以及洪水過程均有較好的模擬精度。作為示例，圖7提供了2場洪水的模擬與實測對比結(jié)果。

圖6 屯溪流域場次洪水模擬精度評定結(jié)果Fig.6 Evaluation results of flood simulation accuracy in Tunxi watershed

圖7 典型場次洪水模擬流量過程Fig.7 Simulated hydrograph of typical flood events

4 結(jié) 論

(1) 本文提出了變動產(chǎn)流層概念，將對壤中流的模擬從單一的土壤層間拓展到土壤-風(fēng)化基巖界面，充分反映了山坡水文試驗的最新認(rèn)識。通過對壤中流形成及發(fā)展過程的精細(xì)刻畫，實現(xiàn)了蓄滿-超滲產(chǎn)流及其轉(zhuǎn)換機制的統(tǒng)一描述。

(2) 基于變動產(chǎn)流層概念建立了分布式水文模型(VRGL)。該模型以DEM網(wǎng)格作為計算單元，通過運動波模擬土壤水濕潤鋒的運動，以描述變動產(chǎn)流層的形成和發(fā)展，并進(jìn)行產(chǎn)流計算；考慮單元間的水量交換，采用二維擴散波進(jìn)行地表匯流計算，河道匯流采用一維擴散波。VRGL共有11個主要的產(chǎn)匯流參數(shù)，且大多物理意義明確，易于通過植被、土壤等公開數(shù)據(jù)產(chǎn)品進(jìn)行確定。在屯溪流域的應(yīng)用結(jié)果表明，該模型對洪峰流量、洪量及洪水過程均具有較好的模擬精度。

本文主旨是提出新的分布式水文模型(VRGL)，限于篇幅，只著重介紹了其產(chǎn)流部分，對匯流等內(nèi)容從簡介紹。本次研究使用屯溪流域的資料進(jìn)行了初步應(yīng)用驗證，后續(xù)將結(jié)合更多流域、更多資料對模型精度提升、模型連續(xù)模擬及流域內(nèi)部狀態(tài)驗證進(jìn)行進(jìn)一步的研究。