張 珂，張企諾，陳新宇，李致家，黃鵬年，姚 成，王 晟，周佳奇

(1.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098；2.河海大學(xué)長江保護與綠色發(fā)展研究院，江蘇 南京 210098; 3.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；4.中國氣象局-河海大學(xué)水文氣象研究聯(lián)合實驗室，江蘇 南京 210098;5.南京信息工程大學(xué)水文與水資源工程學(xué)院，江蘇 南京 210044)

自20世紀50年代以來，我國興建了大量水利工程[1]。水庫、塘壩等蓄水工程阻攔了部分降雨和地表徑流，對場次洪水尤其是中小河流洪水具有明顯的蓄洪削峰作用，從而改變了下游徑流的時空分布[2]。地下水是水資源的重要組成部分，由于地表水資源匱乏，地下水在我國北方尤其華北地區(qū)是主要的供水水源[1]。1980年以后，隨著經(jīng)濟社會高速發(fā)展，人類活動對下墊面的影響日益顯著[3]。地下水的持續(xù)超采造成了地下水位的持續(xù)下降，包氣帶隨之不斷增厚，同時大范圍地下水漏斗的形成增大了降雨徑流過程的入滲損失量，使流域徑流量大幅減少。水利工程和地下水超采等人類活動極大地改變了流域產(chǎn)匯流規(guī)律，因此提出相應(yīng)的模擬方法來科學(xué)刻畫這種特殊降雨徑流關(guān)系變得尤為重要，相關(guān)研究已經(jīng)取得了一些進展。王中根等[4]構(gòu)建了海河流域地表水SWAT與地下水MODFLOW松散耦合模型，模擬結(jié)果符合流域管理的應(yīng)用需求；張建中等[5]建立了模擬地下水埋深變化和灌溉耕地產(chǎn)流特點的蓄滿-超滲產(chǎn)流水文模型，在牤牛河等7個華北平原區(qū)典型流域的應(yīng)用取得了較高的模擬精度；李致家等[6]在新安江模型基礎(chǔ)上增加了對地表和地下攔蓄的模擬，構(gòu)建了半分布式新安江-海河模型，有效地提高了海河典型流域的洪水模擬精度。以上研究為構(gòu)建適用于地表水利工程和地下水超采影響流域的水文模型打下了基礎(chǔ)，但是對上述地表和地下兩種人為干擾的刻畫還不夠精細。

與濕潤流域相比，半濕潤流域降雨時空分布不均，下墊面空間異質(zhì)性對徑流的影響更大，洪水破壞性也更大[7-13]。目前，北方流域的水文模型多為經(jīng)驗?zāi)Ｐ突虬敕植际剿哪Ｐ蚚14-15]，難以精細刻畫流域降雨空間分布不均的特點和土壤植被、地形地貌等下墊面因子的空間異質(zhì)性[16-18]以及地表水利工程控制區(qū)域和地下水埋深的空間分布。針對以上問題，本文在柵格新安江模型[19-20]基礎(chǔ)上，引入地表地下雙人工調(diào)蓄模塊，提出了考慮水利工程和地下水超采影響下的降雨徑流模擬方法，在此基礎(chǔ)上構(gòu)建柵格新安江-地表地下雙人工調(diào)蓄分布式水文模型(Gridded Xin’anjiang-dual anthropogenic aboveground and underground regulation distributed hydrological model，GXAJ-DAR模型)，并檢驗了該模型在受人類活動影響嚴重的華北典型流域的應(yīng)用效果。

1 模型構(gòu)建

1.1 人類活動影響下的降雨徑流模擬

本文研究的人類活動影響主要包括地表水利工程運行與地下水超采，將提出相應(yīng)的方法模擬庫塘閘壩對降雨和地表徑流的直接攔蓄作用以及地下水漏斗對產(chǎn)匯流的間接攔蓄效應(yīng)。

1.1.1地表人工調(diào)蓄模擬

基于水利普查資料，獲取流域庫塘閘壩蓄水能力及其空間分布信息，采用“地表調(diào)蓄水庫”方法刻畫庫塘閘壩對地表徑流的阻攔作用[6,21]。具體計算方法如下：

a.利用ArcGIS提取流域內(nèi)各庫塘閘壩的控制區(qū)域，控制區(qū)域內(nèi)柵格產(chǎn)生的部分地表徑流將被地表調(diào)蓄水庫所阻攔，本研究暫不考慮水利工程的調(diào)節(jié)性能?？刂茀^(qū)內(nèi)的每一個柵格點都會受到水庫的影響，但具體影響機制較為復(fù)雜，因此對其做概化處理，即以水利工程總庫容的一半作為其蓄水容量，將該蓄水容量平均到其控制區(qū)內(nèi)的所有柵格，進而得到控制區(qū)域每一柵格的“水庫蓄水容量”：

(1)

式中：Rc為單個水利工程控制區(qū)域各柵格的水庫蓄水容量，mm；V為水利工程防洪庫容，萬m3；N為控制區(qū)域的柵格數(shù)量；A為柵格面積，km2。

b.以模擬時段初的土壤含水量為指標(biāo)計算地表調(diào)蓄水庫的初始蓄水量：

(2)

式中：Rc0為柵格初始水庫蓄水量，mm；W為模擬時段初柵格的包氣帶土壤含水量，mm；Wm為柵格包氣帶張力水蓄水容量，mm。

c.網(wǎng)格模擬單元在一次洪水過程中的地表徑流攔蓄量估算公式為[1]

Rv=min(Rc-Rc0,Rs)

(3)

式中：Rv為柵格地表徑流攔蓄量，mm；Rs為洪水過程中柵格地表徑流產(chǎn)流量，mm。

非水利工程控制區(qū)域的柵格產(chǎn)流量計算不受地表調(diào)蓄水庫控制；控制區(qū)域內(nèi)進入地表調(diào)蓄水庫的地表徑流將消耗于蒸散發(fā)或用于灌溉等。當(dāng)?shù)乇碚{(diào)蓄水庫蓄滿后，產(chǎn)生的地表徑流將不再受到攔蓄影響。

1.1.2地下人工調(diào)蓄模擬

基于流域的數(shù)字高程地圖、地表水位和地下水位資料，采用“地下調(diào)蓄水庫”方法刻畫產(chǎn)流對地下水超采引起包氣帶變化的非線性響應(yīng)以及地下水漏斗對產(chǎn)匯流的間接調(diào)蓄效應(yīng)[6,21]。具體計算過程分為兩個步驟：

a.當(dāng)下滲雨量滿足包氣帶張力水蓄水容量時，多余的雨量將進入自由水蓄水庫。由于地下水超采的影響，一部分自由水蓄水庫中的水將向下滲漏以補充地下水，抬升地下水位。該滲漏過程采用格林-安普特下滲公式[22-23]進行刻畫，即:

(4)

式中：f(t)為柵格自由水蓄水庫滲漏能力，mm/h；K為柵格飽和水力傳導(dǎo)度，mm/h；ψ為柵格濕潤鋒處土壤吸力，mm；Δθ為柵格土壤飽和含水率與初始時刻含水率之差；F(t)為柵格累計滲漏量，mm；t為時間。

b.自由水蓄水庫的滲漏水量將進入地下調(diào)蓄水庫，累計下滲水量超過地下調(diào)蓄水庫閾值S0的部分以地下水出流形式(出流系數(shù)Kg)形成地下徑流。其中，地下調(diào)蓄水庫閾值S0的計算方法如下：

(5)

式中：Zr為流域出口河道水位，m；Zi為流域出口所在柵格的相對不透水層高程，m；Zg為柵格多年平均地下水位高程，m；ΔH為柵格與流域出口位置的相對高程差，m；ΔH′為模擬時段初地下水位測站幾何中心位置的地下水位高程與其多年平均地下水位高程之差，m；n為柵格土壤孔隙度。

1.2 GXAJ-DAR模型

基于上述產(chǎn)流模擬方法，在柵格新安江模型(Gridded Xin’anjiang model，GXAJ模型)[24]基礎(chǔ)上構(gòu)建適用于受人類活動影響流域降雨徑流模擬的GXAJ-DAR模型。不透水面上的直接徑流與表層自由水蓄水庫蓄滿后產(chǎn)生的飽和地面徑流匯合，水利工程控制區(qū)內(nèi)柵格產(chǎn)生的地面徑流經(jīng)地表調(diào)蓄水庫攔蓄，超出蓄水容量部分以地面徑流形式匯流至流域出口；表層自由水蓄水庫中的部分水量將補充地下調(diào)蓄水庫，達到閾值后，地下調(diào)蓄水庫將以某一出流系數(shù)出流形成地下徑流；模型其余部分與GXAJ模型一致，如圖1所示。與基于泰森多邊形的半分布式新安江-海河模型(XAJ-HH模型)相比，GXAJ-DAR模型考慮了流域降雨空間分布不均的特點和土壤植被、地形地貌等下墊面因子的空間異質(zhì)性以及地表水利工程控制區(qū)域和地下水埋深的空間分布，可以更精細地刻畫人類活動影響下的降雨徑流過程。

圖1 GXAJ-DAR模型結(jié)構(gòu)示意圖

GXAJ-DAR模型在GXAJ模型基礎(chǔ)上新增以下參數(shù)：穩(wěn)定下滲率fc、土壤孔隙度n、濕潤鋒處土壤吸力ψ、水利工程防洪庫容Vf、流域出口位置相對不透水層高程Zi。fc、n、ψ可根據(jù)土壤性質(zhì)得到，Vf可根據(jù)實際水利工程信息獲得，Zi采用人工優(yōu)選法率定得到。姚成[25]研究發(fā)現(xiàn)，在分辨率為1 km的柵格內(nèi)認為蓄水容量分布均勻可以滿足洪水模擬精度要求，因此本研究中的GXAJ-DAR模型將不考慮參數(shù)和變量在柵格內(nèi)部的空間變異性問題。

以下將GXAJ-DAR模型與XAJ-HH模型在研究流域進行應(yīng)用檢驗，比較分析兩個模型的模擬效果。

1.3 模型率定與評估方法

GXAJ-DAR模型的集總式參數(shù)和XAJ-HH模型的參數(shù)均采用人工優(yōu)選法在經(jīng)驗取值范圍內(nèi)進行率定。依據(jù)GB/T 22482—2008《水文情報預(yù)報規(guī)范》，選擇徑流深合格率、洪峰合格率、確定性系數(shù)3個指標(biāo)來綜合評判各模型模擬結(jié)果(其中確定性系數(shù)大于0.5時假定為合格)，并對本文提出的人類活動影響下的降雨徑流模擬方法進行適用性評價。

2 模型驗證

2.1 研究流域概況

海河流域為受地下水超采嚴重影響的典型流域，主要分布于我國渤海以西，黃河以北，位于 35°N～43°N、112°E～120°E之間[26]，流域面積達31.78萬km2，長度為450 km，其中60%為山區(qū)，40%為平原。多年觀測數(shù)據(jù)顯示，海河流域年均降水量僅548 mm，由于氣候、地形等多種因素的影響，夏季暴雨集中，冬春雨雪稀少，年降水量分布具有很強的地域性及季節(jié)差異性，同時還伴隨著連續(xù)枯水及連續(xù)豐水的周期性變化。

選擇河北省保定市清水河流域(38.6°N～39.2°N、114.8°E～115.5°E)作為研究流域，該流域位于海河流域大清河水系，屬于半濕潤流域。流域出口北辛店水文站多年平均徑流量為0.29億m3(1971—2010年)，多年平均降水量為457.1 mm(1971—2010年)，最大年降水量為 833.7 mm。實測最高洪水位為 19.70 m，最大洪峰流量為1966年的710 m3/s，其次是1989年的 686 m3/s(由于洪水漫灘，實際洪峰流量估計值為 1 500 m3/s)。

清水河流域面積1 650 km2，流域內(nèi)設(shè)有雨量站11個，其中有完整水雨情資料的水文站1個，雨量站6個；設(shè)有地下水位測站17個，其中有完整地下水埋深資料的3個(圖2)。北辛店水文站上游建有庫塘閘壩9座，其中中型水庫1座(龍?zhí)端畮?、小型水庫7座(西顯口、北固城水庫等)、塘壩1座(水頭水庫)(圖2)。

圖2 清水河流域概況

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

a.水文氣象資料。清水河流域的降水量、蒸發(fā)量、流量和地下水埋深等水文氣象資料由河北省水文局提供，研究所用資料信息包括日尺度和小時尺度的站點降水量和出口流量、小時尺度的地表水位、月尺度的地下水位以及日尺度的蒸發(fā)資料，資料年限為1980—2002年。依據(jù)2000年前洪峰流量大于25 m3/s、2000年后洪峰流量大于10 m3/s的原則，選取該流域1980—2002年間的12場洪水進行模擬研究(2000年后河道基本斷流)，其中8場用于率定，4場用于驗證。

b.下墊面資料。研究所需流域下墊面信息包括水利工程控制區(qū)域分布和多年平均地下水位相對高程分布(圖3)，后者是基于河北省水文局提供的多年地下水位等值線圖插值得到的相對于3個地下水位測站幾何中心的地下水位高程分布。其余所需下墊面信息與GXAJ模型相同，網(wǎng)格尺度為1 km×1 km。GXAJ-DAR模型的新增土壤參數(shù)如表1所示。

(a) 水利工程控制區(qū)域分布

表1 GXAJ-DAR模型新增土壤參數(shù)值

2.3 模擬精度比較

根據(jù)河北省水文局提供的信息，采用清水河流域1980年后10場洪水資料對現(xiàn)有基于降雨徑流的經(jīng)驗?zāi)Ｐ瓦M行參數(shù)率定和驗證，結(jié)果顯示洪峰合格率為40%，本研究的模擬統(tǒng)計結(jié)果如表2所示。根據(jù)GB/T 22482—2008《水文情報預(yù)報規(guī)范》，當(dāng)徑流深實測值的20%小于3 mm時，即徑流深實測值小于15 mm時，許可誤差取3 mm。對于清水河流域的12場洪水，除1989072101號洪水外，其余場次洪水的實測徑流深均小于15 mm，另外有7場洪水實測徑流深小于3 mm，因此兩個水文模型對于徑流深的模擬誤差可基本控制在3 mm以內(nèi)，徑流深模擬合格率均較高，在率定期與驗證期均達到100%(表2)；盡管如此，由圖4可知GXAJ-DAR模型的徑流深模擬誤差總體而言小于XAJ-HH模型，經(jīng)計算，兩者的平均相對誤差分別為45%和47%。由此可見，GXAJ-DAR模型的徑流深模擬效果稍優(yōu)于XAJ-HH模型。

表2 清水河流域模擬結(jié)果統(tǒng)計

(a) 徑流深模擬誤差

在洪峰模擬方面，XAJ-HH模型在率定期和驗證期的模擬合格率分別為37.5%和25.0%，與現(xiàn)有的經(jīng)驗?zāi)Ｐ拖啾扔兴蛔悖鳪XAJ-DAR模型相較于XAJ-HH模型和經(jīng)驗?zāi)Ｐ驮诤榉迥M效果方面有顯著提升，率定期和驗證期的合格率分別達到62.5%和50.0%(表2)。進一步分析各個場次的洪峰流量模擬誤差發(fā)現(xiàn)，GXAJ-DAR模型的洪峰流量模擬相對誤差明顯小于XAJ-HH模型(圖4)，經(jīng)計算，兩者的平均相對誤差分別為31%和50%。由此可見GXAJ-DAR模型的洪峰流量模擬效果優(yōu)于XAJ-HH模型和經(jīng)驗?zāi)Ｐ汀?/p>

在流量過程模擬方面，驗證期GXAJ-DAR模型確定性系數(shù)大于0.5的洪水模擬場次多于XAJ-HH模型，兩者的確定性系數(shù)合格率分別為50.0%和25.0%(表2)，經(jīng)計算，兩者的平均確定性系數(shù)分別為0.42和-0.08。圖4(c)為兩個模型模擬結(jié)果確定性系數(shù)的分布，其中XAJ-HH模型對第2002062412號洪水的模擬結(jié)果與實測過程相差巨大，故該確定性系數(shù)異常值點(-5.4)未在圖中?？傮w而言，GXAJ-DAR模型的確定性系數(shù)高于XAJ-HH模型，說明GXAJ-DAR模型對流量過程的模擬精度明顯高于XAJ-HH模型。

分別挑選率定期和驗證期洪峰流量最大和最小的兩場洪水做典型模擬流量過程展示，見圖5。對于兩場大洪水(圖5(a)(b))，兩個模型在洪峰、洪量和洪水過程三方面的模擬精度均較好，差異不顯著。對于兩場小洪水(圖5(c)(d))，在率定期，GXAJ-DAR模型的洪峰模擬精度高于XAJ-HH模型；在驗證期，GXAJ-DAR模型的洪峰和洪量模擬精度均明顯高于XAJ-HH模型。由此可見，GXAJ-DAR模型在模擬小洪水方面相比于XAJ-HH模型具有更顯著的優(yōu)勢。

(a) 1989072101號洪水(率定期)

綜上所述，GXAJ-DAR模型的綜合模擬效果優(yōu)于XAJ-HH模型和現(xiàn)有的經(jīng)驗?zāi)Ｐ停蛇_到丙級預(yù)報精度。GXAJ-DAR模型分布式的參數(shù)和結(jié)構(gòu)能更準(zhǔn)確地刻畫半濕潤地區(qū)產(chǎn)匯流過程對流域下墊面特征和降雨、蒸發(fā)、土壤濕度等水文氣象因子的時空響應(yīng)，有效提高洪水模擬精度。

2.4 模型參數(shù)對比

兩個模型的集總式參數(shù)率定結(jié)果如表3所示。清水河流域的土壤含水量較濕潤地區(qū)低，降雨過程中包氣帶往往無法蓄滿，且降雨空間分布不均，因此，產(chǎn)流過程對張力水蓄水容量Wm和自由水蓄水容量Sm的空間分布較為敏感[6]。XAJ-HH模型在各個流域分塊內(nèi)采用同樣的張力水蓄水容量分布曲線和自由水蓄水容量分布曲線來刻畫蓄水容量的空間分布，由于清水河流域地形起伏較大，各個分塊的下墊面特性差別明顯，因此這種方法的計算結(jié)果不能很好地反映流域?qū)嶋H情況；GXAJ-DAR模型則根據(jù)地形指數(shù)和土壤類型來計算每個柵格內(nèi)的平均Wm，在此基礎(chǔ)上根據(jù)植被類型計算柵格內(nèi)的平均Sm，能較好地反映流域下墊面實際物理特性。

表3 清水河流域模型參數(shù)率定值

2.5 地下攔蓄量空間分布

圖6和圖7分別為洪量最大的1989072101號洪水和洪量最小的1997073108號洪水期間地下攔蓄量的空間分布。1989072101號洪水期間，由于累計降水量較大且前期土壤含水量較高(圖6(a))，除流域東南角因初始包氣帶和自由水蓄水庫缺水量相對較大而導(dǎo)致地下攔蓄量相對較小外，流域絕大部分區(qū)域的地下攔蓄量均達到較高水平(圖6(c))，其中流域的西北部的沙壤土區(qū)和部分中部粉壤土區(qū)的下滲能力相對較高，使得這部分的地下攔蓄量又明顯大于其他區(qū)域。

(a) 包氣帶缺水量

(a) 包氣帶缺水量

1997073108號洪水期間的地下攔蓄量空間分布呈現(xiàn)出很大的差異。該場降雨雨量小，降雨初期流域大部分區(qū)域的張力水和自由水蓄水容量均較低(圖7(a)(b))，因此入滲雨量無法滿足下墊面包氣帶和自由水蓄水庫缺水量。由于以上原因，流域的地下攔蓄量呈現(xiàn)較低水平，許多區(qū)域甚至出現(xiàn)無地下攔蓄現(xiàn)象，地下攔蓄作用主要集中在流域的西北區(qū)域和水系附近區(qū)域(圖7(c))，究其原因，這些區(qū)域的初始土壤濕度較大，因此入滲水量足以使其達到田間持水量，多余部分便可進入自由水蓄水庫開始產(chǎn)流與滲漏過程。

3 結(jié) 論

a.在模擬精度方面，XAJ-HH模型在清水河流域的模擬效果與現(xiàn)有的經(jīng)驗?zāi)Ｐ拖啾扔兴蛔悖籊XAJ-DAR模型的徑流深模擬效果稍優(yōu)于XAJ-HH模型，在洪峰流量和流量過程方面的模擬精度較XAJ-HH模型和經(jīng)驗?zāi)Ｐ陀忻黠@提高，尤其在模擬小洪水方面具有更顯著的優(yōu)勢。這說明GXAJ-DAR模型能更精細地刻畫降雨徑流過程對流域下墊面特征和降雨、蒸發(fā)、土壤濕度等水文氣象因子的時空響應(yīng)，有效提高洪水模擬精度。

b.模型結(jié)構(gòu)與參數(shù)方面，GXAJ-DAR模型中柵格尺度的地表和地下雙調(diào)蓄結(jié)構(gòu)能基于水利工程控制區(qū)域和地下水埋深的實際空間分布對地面徑流和地下徑流進行攔蓄；GXAJ-DAR模型對下墊面參數(shù)(Wm、Sm)的計算方法能更好地反映流域下墊面特性的實際空間分布，Sm平均值的率定結(jié)果更為合理。

c.從應(yīng)用層面來看，本文提出的GXAJ-DAR模型能有效提高清水河流域人類活動影響下的洪水模擬精度，為地表地下人類活動雙重干擾下的洪水預(yù)報提供理論和技術(shù)支撐，尤其在我國北方流域具有較大的推廣應(yīng)用價值。