劉玉環(huán)，劉志雨,2，李致家，張 珂，劉開磊

(1. 河海大學水文水資源學院，江蘇 南京 210098；2. 水利部信息中心(水利部水文水資源監(jiān)測預報中心)，北京 100053；3. 淮河水利委員會水文局(信息中心)，安徽 蚌埠 233000)

洪水是中國最常見的自然災害之一，特別是在中國華北地區(qū)(半濕潤)，陡漲陡落的洪水極易造成嚴重的生命和財產(chǎn)損失[1-3]。在實際生產(chǎn)中，洪水防治重要的非工程措施之一就是利用水文模型進行洪水預報。研究發(fā)現(xiàn)，大部分水文模型在濕潤地區(qū)應用效果較好，但是在半濕潤半干旱地區(qū)模擬效果不佳[4-6]。這是由于這類地區(qū)的下墊面和降水時空分布不均勻，超滲產(chǎn)流和蓄滿產(chǎn)流隨時空變化的現(xiàn)象尤為明顯[7]。另外，華北地區(qū)秋冬季降雨較少，用水量較多，地下水開采量大于補給量，地下水位下降。汛期剛開始時，地下水埋深較深，徑流在匯流過程中不斷滲漏，用于抬升地下水位[8]。因此，華北地區(qū)的洪水預報比濕潤地區(qū)更具挑戰(zhàn)性[9]。

目前，經(jīng)過國內(nèi)外學者多年的研究，在通過大量的模型改進后，混合產(chǎn)流問題取得了階段性成果，從單一產(chǎn)流模型如CASC2D模型[10]、SCS模型[11]、TANK模型[12]等，到混合產(chǎn)流模型包括垂向混合徑流模型[13-14]、增加超滲的新安江模型[15-16]等，專門針對海河流域的新安江-海河模型以及河北雨洪模型[17]等，再到靈活框架模型如FUSE[18]、FARM[19]及SCCM[20]等。但是，這些改進后模型的產(chǎn)流模式在計算時都是固定不變的，沒有系統(tǒng)地動態(tài)識別洪水過程中產(chǎn)流模式隨時間變化的情況。

李致家等[21]認為混合產(chǎn)流模擬的主要困難是識別隨時空變化的產(chǎn)流機制。對于一個特定流域的某場洪水來講，流域地形地貌、土壤屬性、植被等下墊面特征是基本不變的，但降雨和土壤含水量等因子會隨時間而動態(tài)變化，并且這類因子的變化與產(chǎn)流模式的空間分布是相互關(guān)聯(lián)的。因此，針對這一現(xiàn)象，本文在TOKASIDE模型(topgraphic kinematic approximation and saturation-infiltration double excess grid-based distributed model)的基礎上，通過降雨和土壤含水量之間的判定關(guān)系，動態(tài)識別蓄滿產(chǎn)流和超滲產(chǎn)流網(wǎng)格，并添加滲漏模塊，構(gòu)建基于物理基礎的分布式蓄超空間動態(tài)組合模型(TOKASIDE-D)，并在我國華北地區(qū)的北辛店流域進行驗證和分析，探究流域產(chǎn)流區(qū)蓄滿、超滲產(chǎn)流模式空間動態(tài)組合方法的應用效果。

圖1 TOKASIDE-D模型結(jié)構(gòu)Fig.1 Structure diagram of TOKASIDE-D model

1 模 型 構(gòu) 建

1.1 TOKASIDE 模型

TOKASIDE模型是劉志雨等[22]提出的基于地形與運動波的蓄滿超滲產(chǎn)流機制的網(wǎng)格化分布式模型。該模型是在TOPKAPI模型的基礎上發(fā)展起來的，進一步考慮了地下水運動、水庫入流與調(diào)蓄計算等模塊以及超滲產(chǎn)流機制。TOKASIDE模型在空間上采用非線性運動波方程法，將每一個計算單元中的水文過程概化為3個“結(jié)構(gòu)上相似”的非線性水庫方程，用以描述水流在土壤層、飽和坡面以及河道中運動波方程的融合。3個方程分別代表上層土壤中的排水、飽和土壤(或不透水土壤層)中的地面徑流和河道徑流。

1.2 TOKASIDE-D模型

TOKASIDE模型在產(chǎn)流計算模塊中雖然考慮了流域超滲產(chǎn)流，但是這一考慮是針對流域內(nèi)所有網(wǎng)格，即扣除蒸發(fā)后的凈雨先完成超滲過程，后繼續(xù)下滲進行蓄滿過程，這與實際情況中蓄滿產(chǎn)流與超滲產(chǎn)流交替主導的現(xiàn)象不符。在華北地區(qū)開展徑流模擬探究時，需要對蓄滿產(chǎn)流和超滲動態(tài)交替以及地下水超采的情形進行考慮。針對這一問題，對TOKASIDE模型進行改進，彌補假設條件中缺失的情況，動態(tài)識別蓄滿產(chǎn)流和超滲產(chǎn)流網(wǎng)格，以提高模型在華北地區(qū)的模擬精度。

在TOKASIDE模型基礎上，通過降雨和土壤含水量之間的判定關(guān)系，動態(tài)識別蓄滿產(chǎn)流和超滲產(chǎn)流網(wǎng)格，添加滲漏模塊，構(gòu)建基于物理基礎的TOKASIDE-D模型，模型結(jié)構(gòu)見圖1(圖中p為降雨量，pe為降雨強度，θt為網(wǎng)格土壤含水率，θS為網(wǎng)格田間持水率，f為穩(wěn)定下滲率，Q為模擬流量)。將流域網(wǎng)格初步劃分為“蓄滿網(wǎng)格”和“超滲網(wǎng)格”兩類。蓄滿網(wǎng)格和超滲網(wǎng)格分別采用TOKASIDE的土壤水非線性水庫方程和Green-Ampt下滲公式[23]計算，另外針對地下水超采出現(xiàn)的滲漏問題，在模型中添加滲漏模塊。

1.2.1 基本假設

基于TOKASIDE模型的假設，根據(jù)華北地區(qū)下墊面情況，現(xiàn)添加以下假設條件：(a)網(wǎng)格內(nèi)的降雨、土壤性質(zhì)、土地利用等特征均勻分布，不存在下滲能力分布曲線。(b)土壤分成2層(飽和層、下滲層)計算。我國北方半干旱地區(qū)的特色是土層厚，年降雨量僅在400～800 mm，并且多集中于7—9月，所以其土壤基本不可能達到飽和狀態(tài)。但其表層土壤還是會達到飽和狀態(tài)，所以此處將土壤分成2層進行模擬計算[24]。蓄滿和超滲發(fā)生在飽和層(土壤水模塊)；下滲層則為無限下滲土層(滲漏模塊)，主要針對地下水超采形成的滲漏問題。(c)蓄超網(wǎng)格劃分：河道附近地勢平坦，土壤水含量相對較高，河道網(wǎng)格(三級以上)全部為蓄滿網(wǎng)格；非河道網(wǎng)格可為蓄滿網(wǎng)格或超滲網(wǎng)格。(d)蓄超網(wǎng)格轉(zhuǎn)化：非河道網(wǎng)格在計算過程中，通過判斷網(wǎng)格土壤含水量是否達到田間持水量、降雨是否超過下滲率，動態(tài)調(diào)整網(wǎng)格的蓄滿或超滲屬性。

1.2.2 土壤水模塊修改

在原始TOKASIDE模型中，網(wǎng)格計算分為非河道網(wǎng)格和河道網(wǎng)格。從圖1可以看出，在2種網(wǎng)格內(nèi)，均需要完成土壤水、地表水計算，河道網(wǎng)格另需進行河道匯流演算。因此，產(chǎn)流機制的轉(zhuǎn)換主要發(fā)生在非河道網(wǎng)格的土壤水模塊。

降雨未發(fā)生時，在計算過程中當非河道網(wǎng)格土壤含水率未達到田間持水率時，為超滲網(wǎng)格，達到田間持水率后轉(zhuǎn)為蓄滿網(wǎng)格。降雨發(fā)生時，降雨強度超過網(wǎng)格下滲率為超滲網(wǎng)格，若網(wǎng)格為蓄滿則修改為超滲；否則原網(wǎng)格屬性不變。蓄滿產(chǎn)流主要包含飽和地面徑流和壤中流，而超滲產(chǎn)流計算產(chǎn)生超滲地表徑流。

a. 當土壤含水率未達到田間持水率時(θ<θs)，計算網(wǎng)格為超滲網(wǎng)格：

(1)

式中：Vs——網(wǎng)格單元的土壤含水量，m3；X——網(wǎng)格單元尺寸，m；Quo、Qus——從上游區(qū)域進入當前網(wǎng)格單元i的地表、地下徑流量，m3；hout——時段末土壤水深度，m；Ks——飽和水力導水率，m/s；θ0——初始土壤含水率，m3/m3；Sf——濕潤鋒面處土吸力，m；I——累積下滲量，m。

如果降雨強度pe≤f，超滲地表徑流為hexf=0(hexf為超滲土壤水出流深度，m)，下滲量全部補充土壤含水量：

hout=hin+pet-hexf

(2)

如果pe>f，超滲地表徑流為hexf=pe-f，下滲量全部補充土壤含水量至田間持水量，產(chǎn)生壤中流hsoil：

hsoil=hout-hmax

(3)

式中：t——時間，s；hin——時段初土壤水深度，m；hmax——土壤水最大深度，m。

b. 當土壤含水率達到田間持水率時(θ≥θs)，計算網(wǎng)格轉(zhuǎn)換為蓄滿網(wǎng)格，其計算流程如下：

(4)

式中：Cs——當?shù)氐膫鲗剩琺/s；αs——土壤特性決定的參數(shù)。

時段初土壤水深度hin=hmax；壤中流hsoil=hout-hmax；飽和地表徑流hexf=hout-hsoil-hmax；更新時段末土壤水深度hout=hmax。

c. 當降雨強度pe>f時，蓄滿網(wǎng)格修改為超滲網(wǎng)格，則計算模塊由蓄滿模塊轉(zhuǎn)為超滲模塊，按照式(1)～(3)進行計算。

1.2.3 滲漏模塊

華北地區(qū)地下水超采嚴重，地下水位埋深大，建模時考慮垂向土壤深層的滲透。當土壤水的含水量超過蓄水能力時，穿過飽和層以上的厚土壤層進行垂直向下運移。假定上層土壤層的下滲率是土壤水含量的一個函數(shù)，則下滲水量為

(5)

式中：Pr——下滲水量，m/s；Ksv——土壤垂向飽和導水率，m/s；θsat——土壤飽和含水量,m3/m3；αp——取決于土壤類型的指數(shù)(沙石:αp?11；泥土:αp?25)。

1.2.4 地表水模塊

土壤表層飽和后的剩余降水構(gòu)成地表水模塊的輸入，地表坡面徑流計算與土壤部分的描述相似，公式為

(6)

式中：Vo——地表徑流量，m3；ro——由水量平衡方程求解得到的余下水量(或土壤水回流量)，m/s；Co——與地表水流曼寧公式相關(guān)的系數(shù)；αo——使用曼寧公式得出的指數(shù)，取3/5。

1.2.5 河道演算模塊

河道徑流的運動波近似方程采用運動學方法進行描述，在這種方法中動量方程通過曼寧公式近似寫成：

(7)

式中：Vc——單元網(wǎng)格的河道水量，m3；rc——旁側(cè)入流，包括地表徑流量和土壤水流量，m2/s；Quc——來自上游單元河道網(wǎng)格的流量，m3/s；S0——河底坡降；nc——河道曼寧糙率，s/m2/3；Cx——濕周，m。

1.2.6 模型參數(shù)

TOKASIDE-D模型的主要參數(shù)包括計算單元土壤類型對應的土壤厚度L(m)、土壤橫向飽和水力導水率Ksh(m/s)、垂向飽和水力導水率Ksv(m/s)，土地利用類型對應的地表曼寧糙率ns，河道分級對應的河道曼寧糙率nc等。Green-Ampt下滲公式的參數(shù)：飽和水力導水率Ks(m/s)，濕潤鋒面處的土壤吸力Sf(m)。在構(gòu)建TOKASIDE-D模型時，將TOKASIDE模型與Green-Ampt下滲公式中物理含義相同的參數(shù)合并，比如：TOKASIDE模型的垂向飽和水力導水率Ksv和Green-Ampt下滲公式的飽和水力導水率Ks，從而減少模型參數(shù)的數(shù)量。TOKASIDE-D是基于物理基礎的水文模型，在進行參數(shù)優(yōu)化率定時，大部分參數(shù)可根據(jù)土壤和土地利用性質(zhì)得來，而其余參數(shù)使用人工優(yōu)選法和SCE-UA自動優(yōu)選法[25]結(jié)合率定得到。

2 研究流域與數(shù)據(jù)資料

2.1 研究流域

北辛店水文站為海河流域大清河水系清水河的控制站，控制流域面積1 650 km2(圖2)。該站實測最高洪水位為19.70 m，最大洪峰流量為710 m3/s，多年平均年徑流量為0.29億m3，多年平均年降水量為457.1 mm。該流域?qū)偬猩缴角皟A斜沖積平原區(qū)，地勢西高東低，緩緩傾斜，土壤以沙壤土為主。流域西部以山地丘陵為主，幾乎全為次生林或次生溫帶灌草叢，東部為平原區(qū)，地勢平坦開闊，土質(zhì)疏松，土層深厚，已開墾為農(nóng)田。

圖2 北辛店流域地形及站點分布Fig.2 Topography and station distribution of Beixindian Watershed

近年來，受人類活動影響，海河流域下墊面條件發(fā)生了較大變化，流域產(chǎn)匯流規(guī)律也發(fā)生一定程度的改變。受地下水過量開采的影響，海河流域地下水埋深較深，不易蓄滿，產(chǎn)流機制呈現(xiàn)出“先超后蓄”的特點，超滲地面徑流產(chǎn)流機制發(fā)揮了重要作用。

2.2 數(shù)據(jù)資料

本文所需下墊面數(shù)據(jù)為數(shù)字高程地圖(DEM)、土壤類型和土地利用類型數(shù)據(jù)。DEM從地理空間數(shù)據(jù)云網(wǎng)站下載，選擇GDEMV2 DEM數(shù)字高程數(shù)據(jù)產(chǎn)品，分辨率為30 m。土壤數(shù)據(jù)來源于聯(lián)合國糧農(nóng)組織(FAO)，1∶100萬土壤柵格數(shù)據(jù)資料。土地利用數(shù)據(jù)可從美國地質(zhì)勘探局(United States Geological Survey，USGS)下載，數(shù)據(jù)為分辨率為1 km的矢量文件。

研究流域水文數(shù)據(jù)包括降雨及流量數(shù)據(jù)(時段為1 h)，由河北省水文水資源勘測局提供。氣溫數(shù)據(jù)來自中國國家氣象數(shù)據(jù)中心。流量數(shù)據(jù)來自北辛店水文站，降雨資料來自黃龍寺、順平等6個雨量站(圖2)。選用北辛店流域1980—2002年的12場次洪水數(shù)據(jù)，其中1980—1996年的9場為率定期，1997—2002年的3場為驗證期。

3 TOKASIDE-D模型的應用檢驗

將TOKASIDE模型、Green-Ampt模型、TOKASIDE-D模型用于北辛店流域進行應用檢驗。依據(jù)傳統(tǒng)的水文模擬與預報精度評定準則，并參考《水文情報預報規(guī)范》[26]的規(guī)定，結(jié)合半干旱地區(qū)洪水特征，選擇以下評價指標：徑流深誤差合格率，該誤差以實測徑流深的20%作為許可誤差判定預報徑流深是否合格；洪峰誤差合格率，以實測洪峰的20%作為許可誤差；峰現(xiàn)時差合格率，以峰現(xiàn)時間小于3 h為許可誤差；平均確定性系數(shù)，評價洪水實測過程與預報過程之間的擬合程度。北辛店流域3種模型參數(shù)率定的結(jié)果見表1，模擬結(jié)果見表2。

表1 北辛店流域各模型的優(yōu)化參數(shù)

表2 北辛店流域3個模型的模擬結(jié)果統(tǒng)計指標

由表2可得：徑流深誤差合格率方面，率定期TOKASIDE-D模型最高(56%)，而Green-Ampt模型最低(33%)，驗證期3個模型均只有1場合格，合格率33%；洪峰誤差合格率方面，率定期TOKASIDE-D模型也是最高(78%)，TOKASIDE模型最低(44%)，驗證期TOKASIDE與Green-Ampt模型合格率均為33%，而TOKASIDE-D模型高于兩者。峰現(xiàn)時差方面，TOKASIDE模型和TOKASIDE-D模型的率定期和驗證期合格率一致；但均低于Green-Ampt模型的合格率，說明TOKASIDE-D模型在模擬洪峰時間方面與原TOKASIDE模型差異不大，土壤入滲后補充土壤水的過程延長了飽和地面徑流的出現(xiàn)時間，進而使得模擬洪峰時間出現(xiàn)滯后。TOKASIDE-D確定性系數(shù)比原模型有提高，說明改進模型對洪水形狀的模擬與觀測值更接近。

此外，選取2個典型的洪水事件1996080420和1980081700進行詳細分析，其洪水實測過程線與模擬過程線如圖3所示。首先，從2個洪水事件的實測值可以看出，2個洪水的洪峰值差不多，前者洪水起漲快，退水慢；后者的洪水過程陡漲陡落，呈尖瘦形態(tài)。從實測過程線可以大致推斷出1996080420洪水的降水量大于1980081700，前者降水更多的是入滲到土壤中，而后者洪水明顯由超滲產(chǎn)流主導，降水入滲量少，基本上發(fā)生了超滲產(chǎn)流。

圖3 北辛店流域部分洪水事件實測過程線與模擬過程線對比Fig.3 Comparison of observed and simulated hydrographs of partial flood events in Beixindian Watershed

對比3個模型的模擬過程線，在1996080420洪水事件中(圖3(a))，TOKASIDE-D模型的模擬洪水過程線在起漲階段以及洪峰值效果最好，模擬洪水的整體形態(tài)與實測洪水最接近。TOKASIDE模型和Green-Ampt模型的模擬洪水起漲速度都很快，但模擬洪峰值都比實測洪峰大得多，且退水速度快，與實測退水過程偏差較大，說明這2個模型在模擬時將大部分凈雨攔在了地表，轉(zhuǎn)為超滲地表徑流，入滲補充土壤含水量的部分偏少；而TOKASIDE-D模型在漲水階段的土壤蓄水量在退水時被“釋放”出來。但相對實測的退水流量還是偏小，這主要是下滲的水量在滲漏模塊的作用下繼續(xù)下滲補充地下水，從而出現(xiàn)部分水量不足的情況。在1980081700洪水事件中(圖3(b))，3個模型的模擬洪水過程線都呈現(xiàn)陡漲陡落的形態(tài)，與實際洪水過程線擬合較好。這是因為1980081700的洪水過程陡漲陡落，呈尖瘦形態(tài)，是比較典型的超滲產(chǎn)流洪水過程，觀察降雨過程，暴雨中心在流域中下游，因此得以快速匯流至流域出口，并沒有太多的下滲損失。另外，洪水發(fā)生在1980年，人類活動和水利工程相對較少，比較接近天然的降雨徑流過程，因此3種模型在該場次洪水模擬的效果相對較好。

與原TOKASIDE模型對所有網(wǎng)格設置超滲“門檻”不同的是，TOKASIDE-D模型通過對每時刻的土壤含水量和降雨等動態(tài)因子進行“監(jiān)控”，動態(tài)調(diào)整網(wǎng)格的產(chǎn)流機制。以1996080420和1980081700這2個典型的洪水事件為例，展示洪水計算過程中產(chǎn)流機制的動態(tài)變化過程，如圖4和圖5所示。從圖4可以看出，1996080420洪水從8月4日20:00開始降水，一直到5日的2:00，降雨8 h，流域基本還是超滲狀態(tài)，在3:00流域出口附近網(wǎng)格基本蓄滿，在隨后的3 h內(nèi)，網(wǎng)格蓄滿數(shù)量大量增加，直到6日5:00達到最大，此時模擬洪水也達到了峰值，隨后，網(wǎng)格土壤含水量開始下降，洪水過程處于退水階段，蓄滿網(wǎng)格數(shù)量逐漸減少。在這里也可以解釋，TOKASIDE-D模型在1996080420洪水的模擬中，網(wǎng)格蓄水程度高，退水時網(wǎng)格中有大量的壤中流和地下徑流產(chǎn)生，從而在退水階段模擬優(yōu)于其他2個模型。從1996080420洪水的實測過程可知其由超滲主導，這一點從圖5也可以看出，整個洪水過程蓄滿網(wǎng)格數(shù)量確實遠遠少于1996080420洪水，超滲網(wǎng)格占主導地位，進而模擬的洪水過程線呈現(xiàn)尖瘦形態(tài)。

圖4 1996080420洪水的蓄滿網(wǎng)格與超滲網(wǎng)格動態(tài)分布Fig.4 Dynamic distribution of saturation-excess and infiltration-excess grids for flood event 1996080420

圖5 1980081700洪水的蓄滿網(wǎng)格與超滲網(wǎng)格動態(tài)分布Fig.5 Dynamic distribution of saturation-excess and infiltration-excess grids for flood event 1980081700

綜上可知，Green-Ampt模型在北辛店模擬效果相比其他2個模型較差，主要因為產(chǎn)流機制不符合當?shù)氐奶卣?。北辛店流域降雨初期土壤水含量低，容易發(fā)生超滲產(chǎn)流，而后土壤入滲能力強，土壤蓄水量大，轉(zhuǎn)為蓄滿產(chǎn)流主導。由于Green-Ampt模型為超滲產(chǎn)流模型，凈雨大多發(fā)生產(chǎn)流，轉(zhuǎn)為超滲地表徑流，導致模擬洪水過程遠大于實測過程。TOKASIDE模型雖然表現(xiàn)尚可，但是其模擬洪峰與實測洪峰相差較大；而TOKASIDE-D模型結(jié)合了蓄滿產(chǎn)流和超滲產(chǎn)流2種產(chǎn)流機制，并且根據(jù)土壤含水量以及降雨的時空分布，動態(tài)調(diào)整產(chǎn)流機制，能更加準確地把握北辛店流域的產(chǎn)流機制的過程，適應性更強。該地區(qū)的滲漏嚴重，雖然增加了滲漏模塊，但部分場次洪水模擬結(jié)果依舊較差，在后續(xù)研究中，進一步探究該流域的存在的問題以及獨有特征，尋找合適的解決方法。

4 結(jié) 語

在基于物理基礎的分布式模型TOKASIDE模型的基礎上，通過降雨和土壤含水量之間的判定關(guān)系，動態(tài)識別蓄滿產(chǎn)流和超滲產(chǎn)流網(wǎng)格，并添加滲漏模塊，構(gòu)建了基于物理基礎的分布式蓄超空間動態(tài)組合TOKASIDE-D模型，并在北辛店流域進行試驗應用。研究結(jié)果表明，與原始模型(TOKASIDE模型、Green-Ampt模型)相比，TOKASIDE-D蓄超空間動態(tài)組合模型在北辛店流域為代表的華北地區(qū)洪水模擬精度有較為明顯的提升。

TOKASIDE-D模型能夠很好地結(jié)合蓄滿模型和超滲模型的優(yōu)點，并且根據(jù)土壤含水量與降雨等因子的變化，動態(tài)調(diào)整產(chǎn)流機制，并針對地下水超采考慮地下水的滲漏。隨著模型結(jié)構(gòu)與原理的改進完善，TOKASIDE-D模型能夠準確捕捉降雨、土壤含水量等動態(tài)因子對產(chǎn)流模式的控制效應，在我國北方地區(qū)的中小流域，有著很大潛力的研究價值以及較高的應用價值。