基于蓄滿超滲時空動態組合的網格新安江模型

杜若愚，姚 成，劉玉環，李致家，張 珂，朱躍龍

(1.河海大學水文水資源學院，江蘇 南京 210098； 2.河海大學計算機與信息學院，江蘇 南京 210098)

產流模式主要分為蓄滿產流和超滲產流[1]，我國南方濕潤地區以蓄滿產流為主，北方干旱地區以超滲產流為主，半干旱半濕潤地區的產流模式兩者兼而有之[2-3]。半干旱半濕潤地區降雨時空分布不均，下墊面條件空間異質性大，產流模式混合多變，即蓄滿產流和超滲產流隨時空變化明顯，導致多數水文模型洪水預報結果欠佳[4-5]。

半干旱半濕潤地區，是濕潤地區與干旱地區的過渡地帶，如果降雨落在具有很高下滲能力的土壤，就會像濕潤地區那樣發生壤中流和地下徑流，土壤蓄滿后發生飽和地面徑流，即蓄滿產流；而當降雨落在下滲能力小，再加上土壤保水能力差，就會發生超滲產流，這就在空間尺度上形成了蓄滿產流與超滲產流同時發生的現象。而半干旱半濕潤地區洪水預報的難點在于確定產流模式，李致家等[6-8]認為產流模式是隨時空動態變化的，并在新安江模型的基礎上添加了超滲產流結構，完善了產流理論，使其能夠適應半干旱半濕潤地區復雜條件的混合產流洪水預報；李致家等[9]對基于子流域的蓄滿超滲空間組合進行研究，分別采用新安江模型和河北雨洪模型計算不同產流區的產流[10]，但這種產流計算為空間靜態組合，因此，Liu等[11-12]研究了基于網格的蓄滿超滲時空組合，考慮了產流模式隨時空變化的特點，基于TOKASIDE模型構建TOKASIDE-D模型，并成功應用于我國華北地區洪水預報。

TOKASIDE-D模型是基于物理基礎的分布式水文模型，其考慮了地形與運動波的蓄滿超滲產流、地下水運動、水庫入流、調蓄計算和滲漏等模塊，核心思路是產流模式隨時空的動態變化。由于網格新安江(Grid-XAJ)模型[13-18]自創建以來在全國多數流域得到了廣泛應用，因此本文借鑒TOKASIDE-D模型中產流模式時空動態組合的思路，構建基于蓄滿超滲時空動態組合的網格新安江(Grid-XAJ-SIDE)模型，并將Grid-XAJ、Grid-XAJ-SIDE和Grid-GA[19](網格格林-安普特)模型應用于半干旱地區綏德流域進行驗證和分析，探討蓄滿、超滲產流模式時空動態組合方法對流域洪水模擬結果的影響。

1 模 型 構 建

1.1 Grid-XAJ模型

Grid-XAJ模型以DEM網格為計算單元，主要分為蒸散發、蓄滿產流、分水源、匯流4個計算模塊[13-18]。模型假設每個網格的降雨和下墊面條件空間分布均勻，不存在張力水蓄水容量分布曲線和自由水蓄水容量分布曲線，只需考慮各要素在不同網格間的變異性。Grid-XAJ模型是蓄滿產流型分布式水文模型，廣泛應用于我國濕潤地區洪水預報。

1.2 Grid-GA模型

Grid-GA模型[19]主要分為蒸散發模塊、產流模塊和匯流模塊3個計算模塊，各個模塊均在以DEM網格為基礎的正交網格內進行。模型認為每個網格的降雨和下墊面條件空間分布均勻，不存在下滲能力分布曲線，Grid-GA模型是超滲產流型分布式水文模型，主要用于干旱半干旱地區洪水預報。

1.3 Grid-XAJ-SIDE模型構建

Grid-XAJ-SIDE模型是在Grid-XAJ模型的基礎上，使用徑流曲線數-地形指數法劃分流域初始蓄超產流區分布，并添加超滲產流模塊，同時在產流計算過程中動態識別蓄滿網格和超滲網格，分別采用Grid-XAJ模型的產流與分水源和Green-Ampt下滲公式計算產流，最后根據網格間的匯流演算次序，依次將不同徑流成分演算至流域出口。

Grid-XAJ-SIDE模型主要分為蒸散發、蓄超網格動態識別、產流、分水源、匯流5個計算模塊，模型計算流程如圖1所示(圖中Pt為降雨量，Ipt為降雨強度，Wt為網格土壤含水量，WM為網格張力水蓄水容量，ft為下滲能力，t為計算時段，n為總時段數)。

圖1 Grid-XAJ-SIDE模型計算流程Fig.1 Calculation flow-process of Grid-XAJ-SIDE model

Grid-XAJ-SIDE模型基于網格新安江模型假設，并添加以下假設條件：①網格內降雨、土壤性質、土地利用等特征均勻分布，不存在下滲能力分布曲線；②流域內網格在計算過程中，依據土壤含水量是否達到田間持水量，降雨是否超過下滲能力，動態識別蓄滿網格或超滲網格；③每個網格的土壤垂向分布均勻。

1.3.1 蒸散發計算

Grid-XAJ-SIDE模型采用3層蒸散發模型計算流域實際蒸散發量。3層蒸散發模型將每個網格的土壤分為上層、下層和深層，對應的土壤蓄水容量分別為WUM、WLM和WDM，具體計算公式詳見文獻[1]。

1.3.2 初始蓄超網格判定

Grid-XAJ-SIDE模型基于徑流曲線數和地形指數，對流域的蓄滿網格和超滲網格的初始空間分布進行確定，徑流曲線數源自SCS模型，表示某種土壤水分條件下的曲線數，主要用來描述地面超滲產流發生的難易程度，徑流曲線數越大，發生超滲徑流的可能性越大。初始蓄超網格判定采用的數據主要包含DEM、土壤類型和土地利用空間分布。劃分方法如下：計算每個網格內的徑流曲線數和地形指數，再根據徑流曲線數進行初始分類[20]，當徑流曲線數小于60時，將網格劃分為蓄滿網格，否則劃分為超滲網格。在此基礎上，再加入地形指數進一步修正，當蓄滿網格的地形指數小于7，將該網格修改為超滲網格；當網格的地形指數大于25，將該網格修改為蓄滿網格，最終得到整個流域的初始蓄超網格空間分布。

1.3.3 蓄超網格動態識別

與初始蓄超網格判定不同，蓄超網格動態識別是Grid-XAJ-SIDE模型在產流計算過程中，根據當前時段網格單元的土壤含水量、降雨強度與土壤下滲能力之間的關系，動態識別蓄滿網格和超滲網格，是一種動態過程；而初始蓄超網格判定是在模型開始計算之前，基于徑流曲線數-地形指數法，根據下墊面信息進行判定，是一種靜態過程。

蓄超網格動態識別的主要原則：當降雨強度大于下滲能力時，該網格為超滲網格；當網格的土壤含水量達到田間持水量后，該網格則為蓄滿網格。其中，流域初始土壤含水量的值由日模型計算得到。

1.3.4 產流及分水源計算

1.3.4.1 蓄滿產流及分水源

Grid-XAJ-SIDE模型對于蓄滿網格采用蓄滿產流模式，即在降雨過程中，直到土壤含水量達到田間持水量才能產流，而在此之前，所有來水均用于補充土壤含水量。將計算時段內網格單元的實測降雨先扣除相應時段的蒸散發量，即可得到實際用于產流計算的時段雨量Pe，則

(1)

式中：R為時段產流量，mm；W0為上一時段土壤含水量，mm。

在Grid-XAJ-SIDE模型中，每個網格單元的R均被劃分為3種徑流成分即地表徑流Rs、壤中流Ri以及地下徑流Rg。與產流計算一樣，在進行分水源計算時，認為每個網格內自由水蓄水容量分布均勻。分水源計算公式為

Ri=KiS

(2)

Rg=KgS

(3)

(4)

式中：Ki為自由水含量對壤中流的出流系數；Kg為自由水含量對地下水的出流系數；S為自由水含量，mm；SM為自由水蓄水容量，mm。

1.3.4.2 超滲產流

對于超滲網格，采用Green-Ampt下滲公式計算網格下滲能力：

(5)

式中：f(t)為土壤下滲能力，mm/h；Ks為飽和水力傳導度，mm/h；Ψ為濕潤鋒處的土壤吸力，mm；Δθ為飽和含水率與初始含水率之差；F(t)為累計下滲量，mm。計算時無需考慮下滲能力分布曲線，根據超滲產流原理，當降雨強度大于土壤下滲能力時，以下滲能力下滲，超出下滲能力的部分降雨形成地表徑流，計算步驟如下：① 根據式(5)計算各網格下滲能力，其中，初始土壤含水率由日模型計算得到，初始土壤F(t=0)賦極小值0.000 01，即認為初始時刻的下滲能力無限大；②計算每個網格的R：

(6)

③計算該網格當前時段的下滲量：i=Pe-R。該網格τ時段的累積下滲量為：F(t=τ)=F(t=τ-1)+i(t=τ)，其中i為當前時段下滲量，mm。將F(t=τ)代入式(5)，計算下一時段的土壤下滲能力。每個網格存在一個最大累積下滲量Fmax，當土壤累積下滲量大于最大累積下滲量時，土壤含水量達到田間持水量，多出的水量形成地下徑流。④重復步驟①，依次計算流域內每個網格在各時段下的徑流量。

1.3.4.3 土壤含水量更新

Grid-XAJ-SIDE模型中，蓄滿網格土壤含水量按照蓄滿產流算法進行更新。超滲網格土壤含水量更新方式：根據超滲產流原理，當降雨強度小于下滲能力時，降雨全部下滲到土壤中，地面不產生積水；當降雨強度大于下滲能力時，降雨分為兩部分，一部分按照下滲能力下滲到土壤中，補充土壤含水量；另一部分則留在了地面，成為超滲地表徑流。設網格當前時段下滲量為i，上一時段上層、下層和深層的土壤含水量分別為WU0、WL0和WD0，三層土壤的蓄水容量分別為WUM、WLM和WDM，按照先上層后下層的次序更新土壤含水量。入滲量優先補充上層土壤含水量，若上層土壤達到蓄滿，多余的水量繼續下滲補充下層土壤含水量，如果下層土壤同樣已經蓄滿，多余的水量繼續補充深層土壤含水量，直至整個土層土壤含水量達到田間持水量。

1.3.5 匯流計算

Grid-XAJ-SIDE模型匯流分為逐網格坡面匯流和河道匯流2個階段計算。坡面匯流采用一維擴散波方程計算，使用基于兩步法MacCormack算法[21]的二階顯式有限差分格式進行擴散波方程組的求解；河道匯流采用基于網格的馬斯京根河道匯流演算法[22]進行計算。

1.3.6 模型參數

Grid-XAJ-SIDE模型以Grid-XAJ模型為基礎，并添加了超滲產流計算模塊，實現了蓄超產流時空動態組合的定量模擬。由于下墊面的空間分布具有異質性，研究流域的網格參數與下墊面的土壤、土地利用、高程等分布具有緊密相關性，本文涉及的3個模型均為分布式水文模型，其網格的上層張力蓄水容量WUM、下層張力蓄水容量WLM、張力水蓄水容量WM、自由水蓄水容量SM、地下徑流出流系數KG、壤中流出流系數KI、土壤總孔隙度θe、土壤有效孔隙度θs、濕潤鋒處土壤吸力ψ和飽和水力傳導度Ks等參數的取值也為分布式。其中WUM、WLM、WM、SM、KG和KI均根據Grid-XAJ參數及其空間分布估計方法進行計算[23-24]；每個網格的θe、θs、ψ和Ks等超滲產流參數，均與土壤類型相關，根據Rawls等[25-27]的研究確定每個網格的參數值；蒸發折算系數K、深層蒸散發系數C、最大累積下滲量IO、地下徑流消退系數CG、壤中流消退系數CI、河網蓄水消退系數CS、滯后時間Lag和河道徑流逐網格馬斯京根法演算參數Kech、Xech等參數則基于流域實測水文資料，采用人工優選法在經驗取值范圍內進行率定。

2 模 型 驗 證

2.1 研究流域

選擇半干旱地區的綏德流域作為研究流域，流域高程、水系及站點分布如圖2所示。綏德水文站位于陜西省綏德縣，東經110°14′，北緯37°30′，建站于1959年6月，流域控制面積3 893 km2。該站多年平均降雨量為486 mm，多年平均氣溫為9.7 ℃，實測最大流量為2 350 m3/s，相應水位為818.59 m。綏德流域地形東西跨度大，主要以山地為主，總體呈現西高東低的趨勢；土壤以壤質砂土為主。

圖2 綏德流域地形、水系及站點分布Fig.2 Topography, river and station distribution of Suide Watershed

2.2 數據資料

Grid-XAJ-SIDE模型所需下墊面數據主要包括DEM、土壤類型和土地利用類型。DEM數據選擇由美國太空總署(NASA)與國防部國家測繪局(NIMA)聯合測量的SRTM數據，分辨率為30 m，通過地理空間數據云(http://www.gscloud.cn)下載。土壤類型數據來源于聯合國糧農組織(FAO)和維也納國際應用系統研究所(IIASA)構建的世界土壤數據庫HWSD (harmonized world soil database version)1∶100萬土壤柵格數據。土地利用類型數據通過美國地質勘探局(United States Geological Survey，USGS)下載，數據是分辨率為1 km 的矢量文件。

以1 h為計算步長建立模型，選取2010—2018 年站點資料較為齊全且洪峰較大的15 場洪水資料進行計算，將收集到的實測降雨和流量資料插值為1h步長數據，其中前10場洪水資料用于模型參數率定，后5場洪水用于模型驗證。

2.3 驗證結果與分析

將Grid-XAJ模型、Grid-GA模型和Grid-XAJ-SIDE模型應用于綏德流域進行應用檢驗和對比分析時，網格大小均采用1 km×1 km，既可以更好地與土壤、植被等數據的網格大小保持一致，也可以在保證應用精度的前提下，提高模型運算效率。Grid-GA模型和Grid-XAJ-SIDE模型對于超滲產流的計算，計算時段會影響模擬結果，但綏德流域原始資料大多為2 h步長，目前受限于流域內觀測資料的精度，3個模型計算步長均選擇1 h。

依據GB/T 22482—2008《水文情報預報規范》[28]評價模型模擬結果，選擇徑流深誤差、洪峰相對誤差和確定性系數等洪水模擬指標進行綜合評價。3個模型人工優選參數見表1，率定和驗證的精度見表2。

表1 綏德流域各模型參數值

表2 綏德流域各模型模擬結果對比

對于徑流深的模擬，率定期Grid-XAJ-SIDE模型、Grid-XAJ模型和Grid-GA 模型模擬合格率分別為90.0%、80.0%和100.0%，驗證期合格率分別為80.0%、80.0%和100.0%，模擬結果均較好，主要由于綏德流域實測徑流深普遍偏小，15場洪水中，只有20170723號洪水實測徑流深大于15 mm，且3個模型對于徑流深的模擬值與實測值誤差多數在3 mm以內，依據GB/T 22482—2008《水文情報預報規范》[28]，3個模型的徑流深合格率都比較高。綏德流域地處半干旱區，洪水歷時短、漲落快、總徑流深小但洪峰較大，因此準確模擬洪峰更為重要。在洪峰合格率方面，Grid-XAJ-SIDE模型率定期和驗證期模擬結果均最好，分別為70.0%和60.0%，相比于Grid-XAJ模型(率定期60.0%，驗證期40.0%)和Grid-GA模型(率定期30.0%，驗證期60.0%)有明顯提高。3個模型對于峰現時間的模擬結果相同。3個模型的確定性系數普遍較差，說明在綏德這樣的半干旱流域，確定性系數不能用來描述陡漲陡落的洪水過程。

20130726號洪水事件中(圖3(a))，Grid-XAJ-SIDE模型和Grid-XAJ模型模擬結果接近實測過程，此次洪水的前期土壤蓄水程度較高，缺水量小，是以蓄滿產流為主導的洪水過程，初始土壤飽和度如圖4(a)所示。而Grid-GA模型模擬結果與實測過程相差較大，這是因為Grid-GA模型認為土壤是半無限土柱，即不會達到“蓄滿”狀態，入滲到土壤中的降雨不會再產生徑流。在20180710號洪水事件中(圖3(b))，Grid-XAJ-SIDE模型和Grid-GA模型模擬結果明顯優于Grid-XAJ模型，這是因為此次洪水的前期土壤蓄水程度低，缺水量大，不易蓄滿，是以超滲產流為主導的洪水過程，初始土壤飽和度如圖4(b)所示。此次洪水過程中Grid-XAJ-SIDE模型和Grid-GA模型在洪水起漲階段模擬結果相近，且與實測過程接近，但在退水階段，Grid-GA模型退水過程過快，與實測過程偏差較大。Grid-XAJ-SIDE模型認為，即使是超滲產流主導的網格，當下滲量補充土壤含水量達到“蓄滿”狀態后，下滲的水量會產生壤中流和地下徑流，因而比Grid-GA模型更接近實測退水過程。

圖3 綏德流域部分洪水實測過程線與模擬過程線對比Fig.3 Comparison of observed and simulated hydrographs of partial flood events in Suide Watershed

圖4 綏德流域部分洪水場次初始土壤飽和度分布Fig.4 Initial soil saturation distribution of partial flood events in Suide Watershed

Grid-XAJ-SIDE模型通過對逐時刻的土壤含水量和降雨等因子進行實時“監控”，動態調整網格的產流方式。以20130726號洪水和20180710號洪水過程為例，展示在模型計算過程中蓄超網格的動態變化，如圖5所示。20130726號洪水過程的整體蓄水程度較高(圖4(a))，容易蓄滿，因此本場洪水從7月26日16:00開始降雨之后，蓄滿網格數量便迅速增多(圖5(a))，到7月27日2:00達到最大，此時模擬洪水達到頂峰，之后進入退水階段。從圖4(b)和圖5(b)可以看出，20180710號洪水過程整體蓄水程度不高，且下滲較小，由于降雨強度大，歷時短，隨著降雨的進行，從開始降雨(20180710T10:00)到洪峰時刻(20180711T14:00)蓄滿網格減少，整個洪水過程蓄滿網格數量少于超滲網格，是以超滲網格為主導的洪水過程，進而模擬的洪水過程也呈尖瘦形態。

圖5 綏德流域部分洪水蓄滿和超滲網格空間動態分布Fig.5 Spatial dynamic distribution of saturation-excess and infiltration-excess grids in Suide Watershed

Grid-XAJ-SIDE模型以網格為計算單元進行產匯流計算，并基于下墊面特性的分布信息推求模型參數的空間分布，因此模型在輸出流域出口斷面流量的同時，可以輸出任意網格單元的流量過程。以綏德流域20180710號洪水過程為例，由圖6可以看出，此次洪水事件前期主要發生在流域上游，隨著洪水過程的進行逐漸轉移至中下游，以洪水過程流量空間動態分布為依據，不僅可以反映流域各點的流量過程，還可以推測暴雨中心的轉移過程，也從側面反映出半干旱地區降雨空間分布的不均勻程度。

圖6 綏德流域20180710號洪水過程流量空間動態分布Fig.6 Spatial dynamic distribution of discharge for flood event 20180710 in Suide Watershed

通過3個模型的對比驗證與分析，Grid-XAJ-SIDE模型結構完善，結合了蓄滿型產流和超滲產流模型的優點，在下墊面條件比較復雜的半干旱地區提升了洪水模擬精度，對流域內的洪水發生過程實現了精細描述。

3 結 語

在網格新安江模型(Grid-XAJ)的基礎上，基于徑流曲線數和地形指數劃分流域初始蓄超分布，構建基于蓄滿超滲時空動態組合的網格新安江(Grid-XAJ-SIDE)模型，并以半干旱地區綏德流域實測數據對模型進行了驗證。研究結果表明，與Grid-XAJ模型和Grid-GA模型相比，Grid-XAJ-SIDE模型根據土壤含水量和降雨等動態產流影響因子的實時變化，更加精確地描述了網格內的產流模式，更好地結合了蓄滿產流模型和超滲產流模型的優點，使洪水模擬精度有了較為明顯的提升，對半干旱半濕潤地區洪水預報具有積極意義。