考慮山坡和洼地水力聯系的半分布式喀斯特水文模型
——以陳旗流域為例

薛冰賢，陳 喜,2，陳 曦，張志才，程勤波

(1. 河海大學 水文水資源與水利工程科學國家重點實驗室, 南京 210098；2. 天津大學 表層地球系統科學研究院, 天津 300072)

我國南方喀斯特地區有著獨特的地貌和水文特征，山坡土層薄、持水能力弱，植被以林地、灌叢為主；洼地土層較厚，多為農田，因此，山坡和洼地滯留水量是植被耗水、農業灌溉用水的重要水源。山坡、洼地土層以下溶溝、溶隙和管道發育，是水土流失的重要通道；暴雨時山坡水流向洼地匯集，極易形成陡漲陡落雨洪過程，引發洪澇災害。因此，喀斯特地區坡、洼系統水文過程的模擬與預測對生態保護、水資源利用和洪澇災害防治具有重要意義。

喀斯特地區不同大小裂隙和管道的導水介質導水性和水流特征顯著不同，因此地下水流模擬通常采用雙重介質模型，如在孔隙介質地下水滲流(基質流)模型基礎上增加管道流計算，常用的有MODFLOW-CFP[1]模型；另一類是把導水裂隙概化為等效多孔介質[2]，并考慮導水介質各向異性模擬喀斯特流域地下水動態或水文循環過程，這類分布式模型通常應用于具有詳細巖溶裂隙分布、結構、導水性等觀測資料的實驗流域。集總式、半分布式水文模型通常在匯流計算層次考慮不同介質對水流的調蓄作用[3]概化水流特征，如把蓄水體概化為兩個或多個水箱[4]，建立基于落水洞的巖溶半分布式水文模型[5]，模擬和預測降雨-徑流響應特征，但對不同介質產流入滲機制和不同地貌單元(如山坡與洼地)以及它們之間水力聯系的模擬功能還不足，且喀斯特地區入滲補給與非喀斯特不同，具有表層巖溶帶裂隙以及漏水洞、漏斗等快速入滲和水流運動通道。因此，發展兼具資料條件和喀斯特地貌特征的半分布式模型具有重要意義。

本文考慮峰叢喀斯特地貌特征，將流域劃分為山坡、洼地兩個單元，并考慮裂隙和管道的產流和入滲機制以及飽和帶快速流和慢速流的交互作用，構建半分布式水文模型。利用流域出口斷面流量以及洼地地下水位等觀測資料進行模型參數率定和驗證，分析坡、洼系統多重介質徑流量分布，為峰叢喀斯特地區水文過程模擬和水資源量評估提供支撐。

1 模型與方法

將峰叢喀斯特流域劃分為山坡(H)、洼地(D)兩個單元(面積分別為AH和AD)，每個單元內根據蓄水、導水能力不同，劃分為具有基質流特征的土壤與細小裂隙區(所占面積比例為α)以及具有管道流特征的大裂隙管道(漏斗)區(所占面積比例為(1-α))；垂向劃分為非飽和帶和飽和帶，進行單元內產流、入滲補給和蓄泄關系以及兩個單元之間水力聯系計算(圖1)。

1.1 非飽和帶產流和入滲補給計算

在具有基質流特征的土壤與細小裂隙區，與南方非喀斯特地區產流機制相似，采用新安江模型蓄水容量曲線[6]計算時段產流R。由于喀斯特裂隙發育深，R一部分向下入滲補給飽和帶，即IS=ks×R(ks為慢速水箱補給系數)；剩余部分(R-Is)側向匯至周邊大裂隙和管道(漏斗)區。

在具有管道流特征的大裂隙管道(漏斗)區，降落在該區的降雨量P扣除少部分損失后(aP，a為折扣系數)形成的自由水(產流)，加上上述土壤與細小裂隙區部分匯入的徑流(R-Is)進入管道；并進行管道蓄量VF,t調蓄計算，管道中一部分蓄量向下補給飽和含水層，即IF=kf×VF,t(kf為快速水箱補給系數)。當管道蓄量VF,t大于其最大蓄量VM時，管道溢出產生地表徑流Rs，即Rs=VF,t-VM；否則Rs=0。

1.2 飽和帶蓄泄關系演算

基質區和管道區入滲補給飽和含水層中慢速、快速水箱，進行蓄泄演算。在兩個單元之間，山坡(H)慢速、快速水箱出流分別流向洼地(D)慢速、快速水箱，再經過洼地調蓄流入地表或地下河出口。模型結構如圖1所示。

山坡(H)、洼地(D)單元慢速或快速水箱蓄泄演算如下。

對于山坡(H)單元，水箱j(j等于F或S，代表快速或慢速水箱)蓄量VHj計算公式為：

(1)

對于洼地(D)單元，水箱j蓄量VDj計算公式為：

(2)

任一山坡或洼地單元i，EXi計算公式[7]為：

(3)

式中：ke為水箱交換系數。

假定任一單元i(i等于H或D，代表山坡或洼地單元)中水箱j的出流Rij與蓄量Vij成線性關系，即：

Rij=ηij×Vij/Ai

(4)

式中：ηij為單元i中水箱j出流系數。

Eg采用阿維里揚諾夫公式[8]計算，并考慮地下水埋深空間分布的非均勻性，通常假定服從伽瑪分布[9]，Eg計算公式為：

(5)

式中：m為計算單元分區數；Dm為最大平均地下水埋深，m；D0為初始平均地下水埋深，m；Dmax為潛水蒸發極限埋深，m；Γ(γ)為伽瑪函數；γ、β分別為伽瑪函數的形狀、尺度參數；kc為蒸發折算系數；Ep為流域潛在蒸發量，mm。

洼地地下水平均埋深D計算公式為：

D=(VDS+VDF)/AD/sy

(6)

式中：sy為給水度，本文取0.05。

模型應用于西南典型喀斯特小流域，面積很小，因此忽略徑流的匯流過程，洼地地表徑流Rs以及慢速、快速水箱出流Rij進入河道從流域出口斷面流出。

2 流域概況及資料

陳旗流域位于黔中丘原盆地區的普定巖溶盆地(26°15′36″ ～26°15′56″、東經105°43′30″ ～105°44′42″)，屬于亞熱帶季風濕潤氣候區，流域面積0.9 km2。利用DEM數據將流域劃分為山坡、洼地兩個單元，面積分別為0.73 km2、0.17 km2。流域多年平均降水量1 300 mm，降雨主要集中在5-8月份，年均氣溫15.1 ℃。具有貴州典型的峰叢山體、峰叢洼地地貌及喀斯特水文地質特征[10]。

流域內在陽坡、火燒坡位置(圖2)分別設立了HOBO小型氣象站，獲取實時降雨、氣溫、氣壓、風速、濕度、太陽輻射等氣象資料，利用彭曼公式計算潛在蒸散發量。流域出口處修建了三角堰并分別放置了HOBO水位自計儀觀測水位，獲取5 min間隔的實時水位數據，通過堰流公式計算獲得流域的流量數據。利用流域洼地內4眼觀測井的水位數據，推求流域洼地地下水埋深。

圖2 陳旗流域圖Fig.2 Chenqi watershed

模型率定期為2016年10月8日至2017年10月30日，歷時9 299 h；驗證期為2017年10月30日至2018年6月12日，歷時5 406 h；計算步長1 h。

3 參數率定

3.1 目標函數

目標函數是為了確定模擬過程與實測過程的吻合程度，首先分別采用流域出口流量和洼地地下水埋深的納什效率系數NSE[11]為單目標函數，再結合流量和地下水埋深的納什效率系數構建多目標函數[12]，計算公式為：

(7)

式中：σo,Q、σo,D分別為實測流量Q、地下水埋深D的標準差；n為數據時序長度；NSEQ、NSED分別為Q、D的納什效率系數。

3.2 參數敏感性分析與優化

采用全局敏感性分析方法，分析各個參數以及它們之間相互作用對模型輸出的影響。采用Monte Carlo法隨機生成各個參數大量的樣本[13](本文生成50000組樣本)，利用蒙特卡洛分析工具箱(MACT)中的RSA方法對參數和目標函數進行歸一化處理，并計算和繪制累計頻率分布圖[14]，利用頻率分布圖計算參數敏感性指數(見表1)。敏感性指數大于0.1時，參數敏感；敏感性指數小于0.04時，參數不敏感[12]。最后采用SCE-UA算法[15]對敏感參數進行優選得到最優參數值，結果見表1。

表1 參數敏感性指數及其率定值Tab.1 Parameter sensitivity index and rating parameters

注：*表示敏感，□表示不敏感。

4 結果分析

由表1中參數的敏感性指數可知，以流量為目標，敏感參數為：土壤與細小裂隙(基質流)面積所占比例(α)、慢速、快速水箱入滲補給常數(ks、kf)、快速水箱出流系數(hF)，以及山坡單元蒸發折算系數(kc)、管道區降水折扣系數(a)。以地下水埋深為目標，敏感參數為：山坡單元的α、a、ks以及洼地單元的hS。采用多目標函數，敏感參數為：山坡和洼地單元α、ks、kf、hS以及山坡單元kc、a。由此可見，無論是以流量或地下水動態變化的單目標，還是反映兩者變化的多目標，影響產流的慢速、快速水箱入滲補給參數(ks、kf)以及基質流面積所占比例(α)為敏感性參數；基質區平均張力水容量(wm)、蓄水容量曲線指數(b)均不敏感。

本文對多目標函數敏感的參數進行優選，得到最優參數值(見表1)。山坡、洼地單元之間差異大的參數包括影響水量平衡的蒸發折算系數(kc)、平均張力水容量(wm)以及影響出流過程的慢速、快速水箱出流系數(ηS、ηF)。相比于山坡，洼地土層厚，巖石裸露率較低，基質區面積(α)大，張力水容量(wm)大；由于地形平緩，出流系數(ηF、ηS)小，慢速與快速水箱交換弱(ke小)。山坡管道區最大蓄量VM達70.75 mm，即次降雨量滿足VM時產生地表徑流，這與該流域坡面集水區實測地表徑流形成的累計降雨深基本一致[16]，說明參數率定結果反映坡地、洼地水文地質特征。

流量模擬過程如圖3，率定期NSEQ為0.92；實測、模擬徑流深分別為300.7 mm、334.0 mm，誤差為11.1%；場次洪水過程漲水、退水過程以及峰現時間基本吻合，其中，洪峰最大的洪水過程峰現時間推遲1h，兩場大洪水洪峰模擬值偏小，誤差分別為-5.6%、-7.7%。驗證期NSEQ為0.91；實測、模擬徑流深分別為175.6 mm、198.0 mm，誤差為12.7%；場次洪水基本吻合，峰現時間基本一致，其中洪峰誤差較大的一場洪水誤差為10.4%。各場次洪水峰現時間與洪峰誤差都在許可誤差[6]范圍以內。根據模擬結果(圖3)，暴雨-洪峰模擬誤差相對較大(2017/06/12、2017/07/09)。分析其原因，快速流系統入滲過程隨降雨強度變化而變化，尤其在暴雨期間，落水洞集中補給地下管道作用增強，而模型對快速流系統入滲過程進行了概化，未考慮入滲隨降雨特征的變化。此外，地下管道匯水過程也隨流量變化而變化，模型未考慮匯流過程也是造成誤差的重要原因之一。

利用洼地地下水埋深資料對模擬結果進行驗證(圖4)，率定期NSED為0.81，平均埋深誤差為0.01 m；驗證期NSED為0.83，平均埋深誤差為0.02 m。模擬與實測過程基本吻合。

圖3 模擬與實測流量過程與洪水過程Fig.3 Simulated and measured discharge process and flood process

圖4 地下水埋深模擬與實測過程線及兩者相關關系Fig.4 Simulated and measured groundwater depth process and correlation

山坡和洼地單元蒸散發量以及產流量模擬結果統計值如表2。山坡林地、灌叢等覆被好，實際蒸散發量大，徑流深小；但山坡面積大，計算期內產水量占總水量的76.0%。流域地表徑流、快速徑流、慢速徑流占總水量比例分別為3.8%、71.8%、24.4%，地表徑流和快速徑流量之和遠比慢速徑流量大，反映了喀斯特地區裂隙管道快速流對陡漲、陡落流量過程的控制作用。

表2 流域水文過程模擬結果Tab.2 Simulation results of Watershed hydrological process

5 結 語

本文針對峰叢喀斯特流域山坡、洼地地貌特征，構建了反映導水介質慢速、快速水流入滲補給和蓄泄特征以及山坡-洼地水力聯系的半分布式水文模型，通過典型流域實測數據對模型進行驗證。研究結果表明：

模型可以較好地模擬喀斯特地區降雨徑流過程、洼地地下水動態變化以及不同地貌單元對流域出口流量的貢獻，參數率定值能夠反映山坡與洼地地形和水文地質特征對產流過程的影響；山坡單元是峰叢—洼地喀斯特地貌區主要水源，約占總產流的3/4；快速徑流約占總產流的2/3，是喀斯特山區陡漲、陡落流量過程主要控制因素。本研究為我國南方喀斯特地區水資源利用和洪澇形成過程預測提供定量分析方法。

模型還有許多不足，模型應用于西南喀斯特小流域，忽略了流域的匯流過程；同時，沒有考慮洼地中農業生產活動對產流和地下水過程的影響，在今后的研究中將加以改進。

□