李彬權(quán)，梁忠民，付宇鵬，王 軍，胡義明

(1. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；2. 河海大學(xué)水科學(xué)研究院，江蘇 南京 211106；3. 水利部珠江水利委員會水文水資源局，廣東 廣州 510611)

在全球氣候變化與高強度人類活動影響的大環(huán)境下，原本依據(jù)天然降雨徑流關(guān)系的水文模型或預(yù)報模式已經(jīng)無法滿足當(dāng)前復(fù)雜下墊面環(huán)境以及對預(yù)報精準(zhǔn)度日益增長的需求[1- 3]。當(dāng)前變化環(huán)境下洪水預(yù)報面臨新的問題和挑戰(zhàn)，主要表現(xiàn)在[4- 5]：① 氣候變化導(dǎo)致降水特征發(fā)生變化，水循環(huán)加快，極端事件增加，沿海防洪能力下降；② 下墊面變化改變了流域產(chǎn)匯流規(guī)律，破壞資料系列的一致性，模型參數(shù)的代表性不復(fù)存在，特別是資料缺乏地區(qū)的洪水預(yù)報問題仍是重大難題。氣候變化和人類活動破壞了水文系列的代表性，資料不能滿足模型參數(shù)率定的要求[6]。因此，亟需創(chuàng)新模型的理論方法，以反映下墊面變化對產(chǎn)匯流過程的擾動；發(fā)展分布式水文模型的物理確定方法，以解決資料缺乏地區(qū)模型參數(shù)的確定難題。

1995年，Garrote和Bras[7]在Cabral下滲理論[8]基礎(chǔ)上，通過研究濕潤鋒演進(jìn)及相對不透水層(變動飽和帶)產(chǎn)生與發(fā)展過程，建立土壤水分剖面變化的數(shù)學(xué)表達(dá)方程，提出了一種描述徑流發(fā)生位置及多源徑流成分的新型產(chǎn)流模式，并逐步發(fā)展成為完整的分布式水文模型RIBS(Real- Time Interactive Basin Simulator)[9]。這種新型產(chǎn)流模式有別于傳統(tǒng)的蓄滿- 超滲2種產(chǎn)流機制水平向或垂向組合方式，避免人為設(shè)定某種產(chǎn)流機制的主觀性，更接近真實的產(chǎn)流過程，但仍難以突破當(dāng)前分布式水文模型的最大應(yīng)用瓶頸，即參數(shù)高度依賴率定，在資料缺乏地區(qū)應(yīng)用精度受限[10- 11]。

本文綜合RIBS模型中的變動飽和帶產(chǎn)流模式和網(wǎng)格水滴匯流方法，構(gòu)建基于單元網(wǎng)格剖分結(jié)構(gòu)的分布式產(chǎn)匯流模型，分析模型參數(shù)敏感性，建立敏感性產(chǎn)流參數(shù)與地形參數(shù)、土壤類型數(shù)據(jù)之間的統(tǒng)計關(guān)系，結(jié)合野外坡面流觀測試驗確定匯流參數(shù)，在多個實際流域進(jìn)行應(yīng)用檢驗。

1 基于變動飽和帶的產(chǎn)匯流模型構(gòu)建

構(gòu)建的模型中產(chǎn)流模塊為RIBS模型中變動飽和帶產(chǎn)流模式[9]，匯流模塊借鑒了文獻(xiàn)[12]提出的“網(wǎng)格水滴”匯流計算思路，由此構(gòu)建得到的產(chǎn)匯流模型結(jié)構(gòu)如圖1(a、b、c分別代表③ 、⑤ 、⑦ 單元內(nèi)的河段)所示。本節(jié)簡要介紹模型產(chǎn)匯流計算的主要過程，詳細(xì)內(nèi)容可參閱文獻(xiàn)[7- 9,12]。

圖1 基于變動飽和帶的產(chǎn)匯流模型結(jié)構(gòu)Fig.1 Framework of the rainfall- runoff model based on variable saturation zone concept

1.1 基于變動飽和帶概念的產(chǎn)流模型介紹

1.1.1 產(chǎn)流模型基本假設(shè)

下滲模型包括非飽和下滲和飽和下滲2個階段。由于飽和水力傳導(dǎo)度隨深度減小，隨著下滲量的累積，土壤將有可能形成變動飽和帶。忽略毛管力作用[13]，則重力是非飽和下滲的主要作用力；同時在場次洪水模擬預(yù)報過程中假定地下水位恒定不變。以下墊面網(wǎng)格單元地表處為原點建立坐標(biāo)系，與地面平行的方向為x方向(坡度下降方向為正)，與地面垂直的方向為z方向(向下為正)。假定土壤飽和水力傳導(dǎo)度(Ks)隨深度呈指數(shù)衰減，即

Ksz=K0zexp(-fd)Ksx=K0xexp(-fd)

(1)

式中：d為土壤深度，mm；Ksz為z方向的土壤飽和水力傳導(dǎo)度，mm/h；K0z為z方向的地表飽和水力傳導(dǎo)度，mm/h；f為飽和水力傳導(dǎo)度隨深度衰減系數(shù)，取值為0～1 mm-1；Ksx和K0x分別為x方向的土壤飽和水力傳導(dǎo)度和地表飽和水力傳導(dǎo)度，mm/h。引入各向異性比率αr=K0x/K0z>1，描述z與x兩方向上飽和水力傳導(dǎo)度的關(guān)系[8]。

采用Brooks和Corey方法[14]計算z、x方向的非飽和土壤水力傳導(dǎo)度分別為：

(2)

(3)

式中：Kz、Kx分別為z、x方向的非飽和土壤水力傳導(dǎo)度；ε為土壤孔隙分布指數(shù)；λ為土壤孔隙分布指數(shù)；Se為關(guān)于土壤含水率(θ)的量綱一變量；θs為飽和土壤含水率，一般取值為0.3～0.6；θr為殘余土壤含水率，一般取值為0.01～0.15。

1.1.2 恒定雨強的下滲模型

在恒定降雨強度條件下，土壤水分剖面通常會經(jīng)歷4個連續(xù)的狀態(tài)過程(圖2)，濕潤鋒的位置深度記為NW，變動飽和帶的上邊界深度記為NU，地下水深度記為NG，4個狀態(tài)下各深度之間關(guān)系為：① 非飽和狀態(tài)NW≤NU

圖2 土柱單元不同水分狀態(tài)示意Fig.2 Schematic diagram of different moisture states of soil column unit

在變動飽和帶尚未產(chǎn)生條件下，飽和水力傳導(dǎo)度隨深度遞減，因此當(dāng)恒定降雨強度I≤K0z時，一定存在一個深度Nb，使得該深度z方向上的飽和水力傳導(dǎo)度與雨強相等，即Kz(θs,Nb)=I， 代入式(2) 得到[8]：

(4)

在此深度之下，飽和水力傳導(dǎo)度小于恒定雨強，濕潤鋒下降速率小于恒定雨強，因而水分在此深度以下開始逐漸積累，產(chǎn)生變動飽和帶。在濕潤鋒到達(dá)Nb之前，濕潤鋒與土柱表面之間處于非飽和狀態(tài)。定義一個比較小的恒定雨量為I0(取I0=0.1 mm/h)，在持續(xù)時間較長的恒定雨強I>I0條件下，土壤水分剖面分成受I影響和濕潤鋒NW之下仍只受I0影響的2個不連續(xù)區(qū)域。

在I>K0z和NW>Nb2種條件下，濕潤鋒以上的土壤會處于飽和狀態(tài)。定義Keq為飽和區(qū)域內(nèi)飽和水力傳導(dǎo)度的調(diào)和平均值：

(5)

NW和NU的計算公式為：

(6)

(7)

式中：α為網(wǎng)格單元地表與水平面的夾角；θI0,NW和θI,NU分別為受I0影響的NW處土壤含水率和受I影響的NU處土壤含水率。

在I

(8)

1.1.3 變化雨強的下滲模型

為簡化變化雨強條件下土壤水分剖面計算的難度，模型假定[7]：在變化雨強下只有第1個濕潤鋒向土壤深處運動，此后下滲水量均補充第1個濕潤鋒，直到土壤達(dá)到完全飽和狀態(tài)。為此，引入等效降雨強度(Ie)，該值將使得濕潤鋒以上非飽和區(qū)域的土壤含水量與變化降雨強度條件下實際土壤含水量相等。

應(yīng)用迭代法求解土壤水分剖面方程可得到Ie的數(shù)值[7]。產(chǎn)生變動飽和帶后，濕潤鋒控制方程2個變量NW、NU可由下式確定：

(9)

(10)

式中：θIe,NU為受Ie影響的NU處土壤含水率。

1.1.4 產(chǎn)流計算公式

在時間步長Δt內(nèi)下墊面網(wǎng)格的地表產(chǎn)流量為

(11)

式中：Rf為地表產(chǎn)流量，mm。

相應(yīng)的壤中流產(chǎn)流量為

(12)

式中：Rint為壤中流產(chǎn)流量，mm。

因此，在假定地下水位不變條件下產(chǎn)流總量為地表產(chǎn)流量與壤中流產(chǎn)流量之和。

1.2 網(wǎng)格水滴匯流模型

所謂“網(wǎng)格水滴”，即將流域上降雨后有產(chǎn)流量的網(wǎng)格視作水滴，在重力驅(qū)動和阻力影響下沿地面坡度方向運動，每個水滴將在何時到達(dá)流域出口斷面取決于路徑長度和在路徑軌跡上運動速度[12]。由于坡面與河網(wǎng)匯流的水力條件有所差別，坡面匯流和河網(wǎng)匯流采用不同的流速公式。

1.2.1 坡面匯流

根據(jù)流域DEM得到地形地貌信息，利用D8法確定流向和匯流路徑長度。坡面流速主要受地形坡度、土地利用、地表糙率等因素影響，采用美國水土保持局提出的坡面流速公式為

(13)

式中：S為坡度；μ為坡度冪指數(shù)；Kv為坡面流速系數(shù)，m/s；ns為地表糙率，取值為0～1；m為糙率冪指數(shù)，一般取值為0.5～1.0。

1.2.2 河網(wǎng)匯流

對于“網(wǎng)格水滴”匯流路徑上的河道網(wǎng)格而言，采用曼寧公式計算水滴的流速，公式為

(14)

式中：Cv為河網(wǎng)流速系數(shù)，m/s；nc為河道糙率；Rc為河道水力半徑。

根據(jù)式(13)和式(14)確定的水滴流速，結(jié)合匯流路徑長度即可得到水滴匯流時間為

(15)

式中：ti為第i個網(wǎng)格水滴的匯流時間，s；di為第i個網(wǎng)格水滴的匯流路徑長度，m；Vi(x)為第i個網(wǎng)格水滴的流速。

2 研究流域

選取20個流域的歷史場次洪水資料用于模型參數(shù)確定和應(yīng)用檢驗，各流域的出口斷面水文站地理位置和流域控制面積見表1。其中，前4個流域?qū)儆谔|苕溪支流余英溪姜灣水文站以上流域(姜灣流域)的子流域，是以和睦橋?qū)嶒炚綶15]為基礎(chǔ)發(fā)展而成的野外山坡實驗流域。姜灣流域4個子流域的DEM空間分辨率為100 m，模型模擬時間步長為30 min；其他流域的DEM空間分辨率為1 km，模型模擬時間步長為1 h。

表1 選用的20個流域地理位置和控制面積

3 模型參數(shù)的確定方法

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

采用普適似然不確定性估計方法(GLUE)[16]分析6個產(chǎn)流參數(shù)的敏感性，先驗分布為均勻分布，似然函數(shù)為確定性系數(shù)(Ens)，其閾值取0.5。各流域的敏感性分析結(jié)果大致相同，參數(shù)K0z和f的較小改變均會導(dǎo)致Ens值的較大改變，θs、θr、ε和αr4個參數(shù)取值變化對模型精度影響不大。因此，參數(shù)K0z和f為敏感性參數(shù)，其他產(chǎn)流參數(shù)不敏感。

3.2 產(chǎn)流參數(shù)的確定

3.2.1 地表飽和水力傳導(dǎo)度

選擇姜灣流域開展野外雙環(huán)入滲觀測試驗，建立地表飽和水力傳導(dǎo)度實測值與地形參數(shù)的統(tǒng)計關(guān)系。在實驗流域河谷內(nèi)選擇40個試驗點，每個試驗點開展5次以上平行的雙環(huán)入滲觀測試驗，分析得到平均穩(wěn)定下滲率，代表試驗點所在下墊面網(wǎng)格(100 m×100 m)地表飽和水力傳導(dǎo)度的平均值。

利用GIS軟件提取得到姜灣流域下墊面網(wǎng)格的3類10個地形參數(shù)：① 一階地形因子：坡度(Tslo)、坡向(Tasp)；② 二階地形因子：剖面曲率(Tpro)、平面曲率(Tpla)、坡度變率(TSOS)、坡向變率(TSOA)；③ 復(fù)合地形因子：地表切割深度(Tcde)、粗糙度(Trou)、地形起伏度(Ttop)、高程變異系數(shù)(Tvar)。利用入滲試驗的前30組試驗資料建立地表飽和水力傳導(dǎo)度與10個地形因子之間的統(tǒng)計關(guān)系，并利用后10組試驗資料進(jìn)行驗證。建立的統(tǒng)計關(guān)系式見式(16)，試驗實測值與公式計算值對比結(jié)果如圖3所示，式(16)的決定系數(shù)為R2=0.76，表明其具有較高的參數(shù)估計精度。

K0z=11 692.152+0.054Tasp-272.04Tpla-640.902Tpro+84.252Tslo-19.338Tcde-
11 529.006Trou-3.354Ttop+135.792Tvar-34.326TSOS+0.894TSOA

(16)

3.2.2 飽和水力傳導(dǎo)度隨深度衰減系數(shù)

參數(shù)f與土壤粒徑隨深度的變化有關(guān)，因此可通過建立f與不同深度土壤粒徑的關(guān)系來確定該參數(shù)值。由于數(shù)據(jù)分辨率限制，難以采用試驗測定方式獲取下墊面網(wǎng)格參數(shù)f的真實值。本文通過模型率定獲取流域空間平均的參數(shù)f值，建立率定的參數(shù)f值與不同深度土壤類型的統(tǒng)計關(guān)系式。土壤資料來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心，包括表層(0～30 cm)和深層(30～100 cm)的砂土、粉砂土、黏土3種土壤類型數(shù)據(jù)。選用表1中前14個流域場次洪水資料率定參數(shù)f值，建立與不同深度土壤類型的統(tǒng)計關(guān)系式見式(17)，決定系數(shù)為R2=0.995；選用表1中后6個流域資料進(jìn)行驗證，對比結(jié)果如圖4所示。

f=-0.008Tsand-0.006 8Tsilt+0.000 01Tclay+0.006 3Ssand+0.004 1Ssilt+0.009Sclay

(17)

式中：Tsand為表層砂土含量， %；Tsilt為表層粉砂土含量， %；Tclay為表層黏土含量， %；Ssand為深層砂土含量， %；Ssilt為深層粉砂土含量， %；Sclay為深層黏土含量， %。

驗證流域的式(17)計算值與參數(shù)f率定值非常接近，相對誤差為-1.6%～6.2%，平均相對誤差絕對值為2.8%，表明該關(guān)系式較合理，結(jié)果具有較高精度。

3.2.3 非敏感性產(chǎn)流參數(shù)

非敏感性產(chǎn)流參數(shù)包括θs、ε、θr和αr。其中，θs和ε來源于根據(jù)遙感資料反演的全球土壤水力參數(shù)數(shù)據(jù)集[17](空間分辨率為30″)，θr可基于表層土壤粒徑分布數(shù)據(jù)利用經(jīng)驗公式[17]估算得到。根據(jù)已有研究經(jīng)驗[7- 9]，參數(shù)αr可取10。

3.3 匯流參數(shù)的確定

3.3.1 坡面匯流參數(shù)

坡面匯流參數(shù)包括Kv、ns、nc、m和μ。其中，ns、nc根據(jù)土地利用類型和文獻(xiàn)[18]中給出的糙率經(jīng)驗值來確定，其他3個參數(shù)根據(jù)姜灣野外實驗流域的坡面流觀測試驗分析確定。選用姜灣流域和睦橋?qū)嶒炚?場降雨(20180724、20180726、20180803、20180813、20180816)7個山坡共35組坡面流速觀測數(shù)據(jù)，根據(jù)式(13)建立方程組，求解得到Kv=0.009 9 m/s，m=1.089，μ=0.65。

3.3.2 河網(wǎng)匯流參數(shù)

河網(wǎng)匯流模塊的唯一參數(shù)為Cv，在假定河道糙率固定條件下取決于Rc，而Rc與水面寬(河寬)有關(guān)。借鑒已有文獻(xiàn)研究成果[19- 20]，在資料缺失地區(qū)忽略河道糙率，僅考慮水力半徑的影響，Cv的估算公式為

(18)

式中：B為流域出口處河寬，m；c1、c2、c3為需要確定的參數(shù)。

采用表1中1—14流域場次洪水資料率定式(18)，后6個流域資料進(jìn)行驗證，對比結(jié)果如圖5所示。流域出口斷面河寬結(jié)合水文站基本情況、野外實測和衛(wèi)星影像綜合估算，率定得到參數(shù)值為c1=0.074 5、c2=466.930 5、c3=0.362 2，關(guān)系式的決定系數(shù)為R2=0.925。驗證流域的式(18)計算值與Cv率定值非常接近，相對誤差均小于12%。

圖3 K0z估算公式在姜灣實驗流域率定與驗證結(jié)果Fig.3 Calibration and validation of the estimate formula of K0z in the Jiangwan experimental watershed

圖4 飽和水力傳導(dǎo)度隨f估算公式在20個流域的率定與驗證結(jié)果Fig.4 Calibration and validation of the estimate formula of coefficient of attenuation of f in 20 basins

圖5 Cv估算公式在20個流域的率定與驗證結(jié)果Fig.5 Calibration and validation of the estimate formula of Cv in 20 basins

4 模型的應(yīng)用驗證

從場次洪水過程模擬角度檢驗?zāi)Ｐ蛻?yīng)用效果，包括2個敏感性產(chǎn)流參數(shù)(K0z、f)確定方法與模型率定參數(shù)結(jié)果的對比分析。① 地表飽和水力傳導(dǎo)度確定方法驗證：選擇姜灣流域梅家塘、范塢里、古竹灣、佛堂村4個子流域2018年7月至2020年8月共16場較大洪水進(jìn)行模擬(11場用于率定，5場用于驗證)。② 飽和水力傳導(dǎo)度隨深度衰減系數(shù)確定方法驗證：選擇表1中后6個流域的歷史場次洪水進(jìn)行模擬。模型模擬效果采用Ens、洪峰相對誤差(Epeak)、洪峰滯時(ΔTpeak)和洪量相對誤差(Evol)4個指標(biāo)進(jìn)行評定。

4.1 參數(shù)K0z確定方法的驗證

對比分析姜灣流域2套場次洪水模擬結(jié)果，見表2：① 參數(shù)K0z根據(jù)式(16)確定；② 參數(shù)K0z來源于全球土壤水力參數(shù)數(shù)據(jù)集[17]；其他參數(shù)采用統(tǒng)一的率定值。對比各斷面場次洪水模擬結(jié)果發(fā)現(xiàn)，除洪峰滯時指標(biāo)外，式(16)的模型結(jié)果比遙感數(shù)據(jù)集中K0z值的結(jié)果精度更高，表明本文提出的參數(shù)K0z確定方法較為合理。圖6給出了佛堂村子流域的驗證場次20190901號洪水過程模擬結(jié)果，模擬時段步長Δt=30 min。

表2 參數(shù)K0z確定方法在姜灣流域場次洪水模擬中驗證結(jié)果

圖6 佛堂村子流域20190901號洪水過程模擬結(jié)果Fig.6 Simulated hydrographs for the flood #20190901 of the Fotangcun sub- basin

4.2 參數(shù)f確定方法的驗證

對比七鄰等6個流域2套場次洪水模擬結(jié)果，見表3：① 參數(shù)f根據(jù)式(17)確定，其他參數(shù)通過率定獲得；② 所有模型參數(shù)均通過率定獲得。結(jié)果表明，6個流域的率定參數(shù)模擬結(jié)果均略優(yōu)于式(17)的模型結(jié)果；總體上看，本文提出的參數(shù)f確定方法對應(yīng)的場次洪水模擬的平均確定性系數(shù)為0.83、洪峰和洪量誤差的平均絕對值分別為10.07%和6.86%、洪峰滯時誤差的平均絕對值為2.61 h，滿足精度要求。圖7給出了七鄰、東灣流域2場驗證場次洪水的模擬結(jié)果，模擬時段步長Δt=1 h。

表3 參數(shù)f確定方法在6個流域場次洪水模擬中驗證結(jié)果

圖7 七鄰流域和東灣流域洪水過程模擬結(jié)果Fig.7 Simulated the flood hydrographs the Qilin Basin and the Dongwan Basin

5 結(jié) 論

利用姜灣野外實驗流域和16個實際流域的資料，建立地表飽和水力傳導(dǎo)度和飽和水力傳導(dǎo)度隨深度衰減系數(shù)2個敏感性產(chǎn)流參數(shù)與下墊面特征數(shù)據(jù)的定量統(tǒng)計關(guān)系，并建立匯流參數(shù)的經(jīng)驗估計關(guān)系式。主要結(jié)論為：

(1) 姜灣流域4個干流斷面的實例應(yīng)用表明，利用下墊面地形特征參數(shù)估算地表飽和水力傳導(dǎo)度的方法在洪水模擬中精度較高，16場洪水的確定性系數(shù)、洪峰(洪量)誤差和洪峰滯時誤差均滿足規(guī)范的精度要求；與使用遙感資料推求地表飽和水力傳導(dǎo)度值得到的洪水模擬結(jié)果相比，除洪峰滯時誤差略差外，其他3個精度指標(biāo)均更優(yōu)。

(2) 根據(jù)20個實際流域的資料對飽和水力傳導(dǎo)度隨深度衰減系數(shù)的估計關(guān)系式進(jìn)行率定和驗證，結(jié)果表明利用土壤質(zhì)地數(shù)據(jù)估算深度衰減系數(shù)值的方法在場次洪水模擬中的精度指標(biāo)均滿足規(guī)范要求，但比單純進(jìn)行模型參數(shù)率定的結(jié)果略差。

(3) 提出的利用下墊面特征數(shù)據(jù)估算模型參數(shù)的方法在資料缺乏地區(qū)場次洪水模擬預(yù)報中具有適用性。