山丘區(qū)變動(dòng)產(chǎn)流層分布式水文模型

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

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

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

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

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

1 變動(dòng)產(chǎn)流層

1.1 已有認(rèn)識(shí)

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

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

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

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

1.2 概念提出

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

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

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

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

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

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

2.1 模型概述

變動(dòng)產(chǎn)流層分布式水文模型(VRGL)是一個(gè)充分反映土壤-風(fēng)化基巖界面作用的分布式水文模型，通過對壤中流形成及發(fā)展過程的精細(xì)刻畫，實(shí)現(xiàn)了蓄滿-超滲產(chǎn)流及其轉(zhuǎn)換機(jī)制的統(tǒng)一描述。VRGL以DEM網(wǎng)格為計(jì)算單元，在計(jì)算單元上構(gòu)建產(chǎn)流模塊，模擬植被截留層、土壤層和風(fēng)化基巖層的產(chǎn)流過程并考慮單元間的水量交換，采用二維擴(kuò)散波和一維擴(kuò)散波分別進(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為時(shí)間，h；P為降雨強(qiáng)度，mm/h；Ev為植被截留層的實(shí)際蒸散發(fā)率，mm/h；Pv為穿透雨強(qiáng)，mm/h。

2.2.2 土壤層產(chǎn)流

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

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

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

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

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

域Ⅰ的水量平衡方程為

(2)

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

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

(3)

式中：Rg為域Ⅰ的等價(jià)恒定雨強(qiáng)[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)時(shí)(即θn=Nsd=θs)，將形成土巖界面變動(dòng)產(chǎn)流層，根據(jù)式(3)和式(4)可推導(dǎo)出對應(yīng)的土壤含水總量臨界值的計(jì)算公式為

(5)

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

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

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

等價(jià)恒定雨強(qiáng)公式為

(6)

地表下滲強(qiáng)度公式為

Ig=Rgcosα

(7)

壤中流強(qiáng)度公式為

qg=NsdRg(ar-1)sinαL

(8)

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

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

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

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

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

(9)

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

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

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

(10)

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

(11)

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

等價(jià)恒定雨強(qiáng)公式為

Rg=K0nexp(-fNwt)

(12)

地表下滲強(qiáng)度公式為

(13)

壤中流強(qiáng)度公式為

(14)

飽和地面徑流強(qiáng)度公式為

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

(15)

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

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

域Ⅱ的水量平衡公式為

(16)

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

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

(17)

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

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

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

等價(jià)恒定雨強(qiáng)公式為

(18)

濕潤鋒運(yùn)動(dòng)公式為

(19)

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

地表下滲強(qiáng)度公式為

Iu=(Ru-Rg)cosα

(20)

壤中流強(qiáng)度公式為

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

(21)

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

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

蒸發(fā)強(qiáng)度公式為

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

(22)

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

等價(jià)恒定雨強(qiáng)公式為

Ru=K0nexp(-fNst)

(23)

濕潤鋒運(yùn)動(dòng)公式為

(24)

土壤層間水位深度公式為

(25)

地表下滲強(qiáng)度公式為

(26)

壤中流強(qiáng)度公式為

(27)

飽和地面徑流強(qiáng)度公式為

(28)

蒸發(fā)強(qiáng)度計(jì)算同公式(22)。

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

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

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

(29)

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

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

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

(30)

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

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

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

(31)

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

(32)

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

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

(33)

2.3 匯流模型

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

3 模型應(yīng)用

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

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

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

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

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

3.2 結(jié)果及分析

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

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

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

模型精度評定結(jié)果如圖6所示。從圖6(a)可以看出，不管是率定期還是驗(yàn)證期，洪峰流量和洪量的相對誤差均在±20%范圍內(nèi)，其中，洪峰流量相對誤差的絕對值平均為5.5%(率定期為5.2%、驗(yàn)證期為6.6%)，洪量相對誤差的絕對值平均為7.4%(率定期為7.3%、驗(yàn)證期為7.7%)。在洪水過程的評定方面，從圖6(b)可以看出，率定期和驗(yàn)證期ENS均值分別為0.83和0.85，EKG均值分別為0.83和0.82，ENS和EKG2個(gè)指標(biāo)的評價(jià)結(jié)果較為一致；另外，除20140513場次ENS為0.56外，其余場次ENS均大于0.7，率定期和驗(yàn)證期的ENS最大值可達(dá)0.95和0.91。上述結(jié)果表明，VRGL對洪峰流量、洪量以及洪水過程均有較好的模擬精度。作為示例，圖7提供了2場洪水的模擬與實(shí)測對比結(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) 本文提出了變動(dòng)產(chǎn)流層概念，將對壤中流的模擬從單一的土壤層間拓展到土壤-風(fēng)化基巖界面，充分反映了山坡水文試驗(yàn)的最新認(rèn)識(shí)。通過對壤中流形成及發(fā)展過程的精細(xì)刻畫，實(shí)現(xiàn)了蓄滿-超滲產(chǎn)流及其轉(zhuǎn)換機(jī)制的統(tǒng)一描述。

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

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