張 珂，周佳奇，張企諾，陳新宇，晁麗君，姚 成，李致家

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

全球巖溶面積約為220萬km2，巖溶區(qū)居住人口約10億人。中國西南部是世界上巖溶發(fā)育最典型的地區(qū)之一，具有集中連片的巖溶地貌，分布面積約54萬km2[1-2]。由于石漠化、植被稀少等原因，巖溶流域的土層厚度、土地覆蓋類型等下墊面條件空間異質(zhì)性大，流域內(nèi)巖溶裂隙發(fā)育程度不同，導(dǎo)致了流域產(chǎn)匯流機制復(fù)雜，洪水模擬難度大。巖溶流域中的表層巖溶帶以及地下深層溶隙、溶洞等都是地表水入滲、地下水賦存和運移的良好介質(zhì)，在產(chǎn)匯流過程中起著重要作用，探究地表水-巖溶水-土壤水的轉(zhuǎn)化和產(chǎn)匯流規(guī)律十分必要。許波劉等[3]將巖溶裂隙概化為5種，推導(dǎo)了5種裂隙的瞬時單位線公式，改進(jìn)了新安江模型的匯流；Zhou等[4]將集總式新安江模型和兩種基于蓄水庫的巖溶模型結(jié)合起來，建立了巖溶新安江模型(K-XAJ)，在漓江流域的日徑流模擬中取得了較好的效果；Jeannin[5]在Holloch洞穴根據(jù)觀測到的水位、流量、巖溶管道規(guī)模以及管道長度建立了水動力模型，該模型可以很好地模擬以管道流為主的巖溶系統(tǒng)；張志才等[6]根據(jù)喀斯特流域多孔介質(zhì)與裂隙水流特征，改進(jìn)了分布式水文-植被-土壤模型(DHSVM)，并在貴州普定喀斯特水文試驗站取得了較好的模擬效果。隨著遙感技術(shù)在水文領(lǐng)域的不斷應(yīng)用[7-8]，一些學(xué)者通過遙感影像分析了巖溶流域石漠化情況[9]。大部分學(xué)者將不同發(fā)育程度的巖溶裂隙概化為地下溶隙線性水庫來構(gòu)建集總式巖溶水文模型，而構(gòu)建的分布式概念性巖溶水文模型主要用于日徑流和月徑流的模擬，構(gòu)建的巖溶水力學(xué)模型需要較為精細(xì)的資料，一般用于實驗小流域的水文模擬。此外，還有一些研究將機器學(xué)習(xí)和智能模型用于水文模擬[10-11]。

巖溶流域下墊面空間異質(zhì)性較大，其洪水過程的精確模擬難度大。集總式的巖溶水文模型難以刻畫流域下墊面條件的復(fù)雜性[3-4]；水動力巖溶模型在實驗小流域模擬效果好，但是需要精細(xì)的數(shù)據(jù)資料，且對于面積更大的中小河流流域，模型的洪水模擬計算量龐大[5-6,12-13]。分布式和半分布式概念性巖溶模型可以在考慮下墊面異質(zhì)性的同時簡化計算，但是目前分布式巖溶模型大部分被用于日徑流和月徑流的模擬，用于場次洪水模擬的分布式巖溶模型較少[14-17]。針對以上問題，本文在柵格新安江模型(grid-Xin’anjiang model, GXAJ模型)[18]基礎(chǔ)上，加入了表層巖溶帶自由水蓄水庫和深層地下溶隙線性水庫，構(gòu)建了柵格巖溶分布式水文模型(grid-based karst distributed hydrological model, GKDHM模型)，并檢驗了該模型在典型巖溶流域(普廳河流域)的應(yīng)用效果。

1 研究方法

1.1 巖溶地貌影響下的降水徑流模擬

本文研究的巖溶地貌對流域水文過程的影響主要包括表層巖溶帶的調(diào)蓄過程以及深層地下溶隙對地下徑流的調(diào)蓄作用。通過構(gòu)建表層巖溶帶自由水蓄水庫和深層地下溶隙線性水庫刻畫巖溶地區(qū)的產(chǎn)匯流過程。

1.1.1表層巖溶帶與深層地下溶隙調(diào)蓄

圖1為降水經(jīng)過土壤和表層巖溶帶調(diào)蓄形成水流的示意圖。表土層水分下滲進(jìn)入表層巖溶帶，表層巖溶帶水分向下滲透，部分進(jìn)入深層地下水系統(tǒng)，部分由于表層巖溶帶下部垂向滲透性減小而形成側(cè)向水流[19]。

圖1 土壤與表層巖溶帶水流示意圖Fig.1 Schematic diagram of water flow insoil and epikarst

基于流域巖溶發(fā)育情況和石漠化情況，采用圖2所示的表層巖溶帶自由水蓄水庫模擬表層巖溶帶蓄水作用以及出流過程。通過流域巖溶發(fā)育情況判斷柵格內(nèi)是否有巖溶發(fā)育，無巖溶發(fā)育的柵格采用圖2(a)的自由水蓄水庫來刻畫，對于巖溶發(fā)育柵格，存在石漠化(圖2(b))和未石漠化(圖2(c))兩種情況，壤中流(QI)為上層土向下層土下滲過程中由于滲透性變?nèi)醵a(chǎn)生的側(cè)向徑流[20]，石漠化巖溶區(qū)土壤層較薄，故認(rèn)為石漠化巖溶柵格無壤中流，未石漠化巖溶柵格存在壤中流[4]。根據(jù)流域石漠化情況，石漠化柵格采用圖2(b)的自由水蓄水庫來刻畫，未石漠化的柵格采用圖2(c)的結(jié)構(gòu)。

圖1中表土層的調(diào)蓄過程通過圖2中蓄水容量為SM的自由水蓄水庫模擬，圖1中表層巖溶帶的調(diào)蓄過程通過圖2中蓄水容量為SMEK的自由水蓄水庫模擬，表土層位于表層巖溶帶上部，蓄水容量為SM和SMEK的兩個自由水蓄水庫為串聯(lián)關(guān)系，水流經(jīng)表層土自由水蓄水庫滲入表層巖溶帶自由水蓄水庫。圖1深層飽和水流帶中發(fā)育的深層地下溶隙，許多研究將其概化為地下溶隙線性水庫[3,21]，本文采用同樣的方法構(gòu)建圖2(b)和圖2(c)中消退系數(shù)為CKG的地下溶隙線性水庫，模擬深層地下溶隙對地下徑流的調(diào)蓄作用。

(a) 非巖溶區(qū) (b) 石漠化巖溶區(qū) (c) 未石漠化巖溶區(qū)圖2 表層巖溶帶自由水蓄水庫示意圖Fig.2 Schematic diagram of freewater reservoirs of epikarst

a.降水經(jīng)過植被截留，滿足張力水蓄水容量后，剩余的雨量進(jìn)入土壤自由水蓄水庫。在巖溶區(qū)域內(nèi)土壤自由水通過滲透過程補給表層巖溶帶自由水蓄水庫，滲透過程用格林-安普特公式描述[22-23]：

(1)

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

b.表層巖溶帶自由水蓄水庫對下滲的土壤自由水進(jìn)行調(diào)蓄。定義表層巖溶帶側(cè)向徑流量為REKI，表層巖溶帶自由水蓄水庫對地下水的日出流系數(shù)為KEKG，REKI和巖溶區(qū)域地下徑流RG的計算公式如下：

(2)

式中：KEKI為表層巖溶帶自由水蓄水庫對表層巖溶帶側(cè)向出流的日出流系數(shù)；SEK為表層巖溶帶自由水蓄水量，mm。

c.柵格內(nèi)表層巖溶帶側(cè)向徑流采用線性水庫進(jìn)行匯流計算。對于存在巖溶發(fā)育的柵格，引入巖溶地下徑流消退系數(shù)，反映深層地下溶隙與非巖溶區(qū)對地下徑流調(diào)蓄的差異，計算公式為

QEKI,t=CEKIQEKI,t-1+(1-CEKI)REKIU

(3)

QG,t=CKGQG,t-1+(1-CKG)RGU

(4)

式中：QEKI，t、QEKI，t-1分別為t和t-1時刻表層巖溶帶側(cè)向出流，m3/s；QG，t、QG，t-1分別為t和t-1時刻地下徑流，m3/s；CEKI為表層巖溶帶側(cè)向出流消退系數(shù)；CKG為巖溶區(qū)地下徑流消退系數(shù)；U為單位換算系數(shù)。

1.1.2表層巖溶帶蓄水容量空間分布

20世紀(jì)70年代Mangin首先用表層巖溶帶代表巖溶包氣帶上部相對含水比較豐富的部分[24]。隨著對表層巖溶帶的深入研究，表層巖溶帶成為巖溶學(xué)的重要概念。覃小群等[25]認(rèn)為表層巖溶帶是地表和地下各種巖溶形態(tài)構(gòu)成的不規(guī)則帶狀強巖溶化層。

表層巖溶帶在巖溶發(fā)育區(qū)域普遍存在。對于表層巖溶帶的厚度，眾多研究表明，在地形平緩的地帶，如洼地、巖溶谷地和岸坡等地帶，表層巖溶帶發(fā)育程度較好；在山頂山脊等陡坡地帶，表層巖溶帶發(fā)育規(guī)模較小[2,25-26]。

水的侵蝕性是影響表層巖溶帶空間差異的重要因素。地形影響區(qū)域水循環(huán)過程，進(jìn)而影響巖溶動力過程，對表層巖溶帶發(fā)育產(chǎn)生重要影響。山坡土壤水有向洼地匯集的趨勢，在地勢低凹地區(qū)土壤含水量相對較高，增加了水流溶蝕時間，促進(jìn)了表層巖溶帶向深部發(fā)育[2]。在巖溶峽谷兩側(cè)地帶和緩坡部位由于其邊坡的風(fēng)化裂隙較其他地貌部位發(fā)育強烈，特別是河流兩岸岸壁臨空面，大氣降雨可迅速進(jìn)入巖體內(nèi)，大量水流快速進(jìn)入巖體具有強大的侵蝕力，使得該區(qū)域表層巖溶帶發(fā)育較深。在山頂山脊及陡坡地帶，地形陡峭，降雨入滲條件差，雨水多形成坡面流流失，雨水侵蝕的概率較小、時間較短，這些地段表層巖溶帶發(fā)育較淺[26]。

地形指數(shù)反映的是徑流在流域中的累計趨勢以及重力使徑流順坡移動的趨勢[27]，平緩地區(qū)的地形指數(shù)大，在山頂山脊等陡坡地帶地形指數(shù)較小[28]。在巖溶發(fā)育區(qū)域內(nèi)建立表層巖溶帶與地形指數(shù)的線性關(guān)系，假設(shè)巖溶發(fā)育區(qū)域內(nèi)地形指數(shù)最大的柵格對應(yīng)的表層巖溶帶蓄水容量最大，地形指數(shù)最小的柵格對應(yīng)的表層巖溶帶蓄水容量最小，即

(5)

式中：TImax、TImin分別為巖溶發(fā)育區(qū)域內(nèi)柵格地形指數(shù)中的最大和最小值；SMEKmax、SMEKmin分別為表層巖溶帶自由水蓄水容量最大和最小值，mm；ζa、ζb為線性關(guān)系的兩個系數(shù)。

基于數(shù)字高程計算得到的每個柵格的地形指數(shù)TI，通過率定SMEKmax和SMEKmin并結(jié)合式(5)計算得到ζa與ζb后，基于下式可以計算每個柵格的表層巖溶帶自由水蓄水容量SMEK：

SMEK=ζaTI+ζb

(6)

圖3 GKDHM模型結(jié)構(gòu)Fig.3 GKDHM structure

1.2 GKDHM模型

許多研究表明巖溶流域產(chǎn)流機制為蓄滿產(chǎn)流[29-30]，部分適用于巖溶流域降水徑流模擬的水文模型是在新安江模型基礎(chǔ)上改進(jìn)而來的[3-4,21]。本文在GXAJ模型的基礎(chǔ)上，通過引入表層巖溶帶自由水蓄水庫構(gòu)建了GKDHM模型，其結(jié)構(gòu)見圖3。石漠化巖溶區(qū)土壤層薄，采用二層蒸散發(fā)計算方法模擬蒸散發(fā)[20,31]，未石漠化巖溶區(qū)和非巖溶區(qū)采用三層蒸散發(fā)計算方法模擬蒸散發(fā)。巖溶發(fā)育區(qū)域土壤下滲補給表層巖溶帶自由水的過程采用格林-安普特公式模擬。由于石漠化區(qū)域土壤發(fā)育差，故認(rèn)為石漠化巖溶區(qū)無壤中流，巖溶區(qū)域深層地下溶隙對匯流過程的影響通過引入考慮巖溶區(qū)域地下水消退系數(shù)的地下溶隙線性水庫進(jìn)行模擬。

圖3中地面徑流1、壤中流1和地下徑流1分別為非巖溶區(qū)域的地面徑流、壤中流以及地下徑流，對應(yīng)圖2(a)中的地表徑流QS、壤中流QI和地下徑流QG。圖3中的未石漠化土壤壤中流2是巖溶區(qū)中未石漠化地區(qū)的壤中流，對應(yīng)圖2(c)中的壤中流QI。圖3中地面徑流2、表層巖溶帶側(cè)向出流和地下徑流2分別是巖溶區(qū)的地面徑流、表層巖溶帶側(cè)向出流和地下徑流，對應(yīng)圖2(b)和(c)中的地表徑流QS、表層巖溶帶側(cè)向出流QEKI和地下徑流QG。

GKDHM模型在GXAJ模型基礎(chǔ)上新增了以下參數(shù)：用于格林-安普特公式的飽和水力傳導(dǎo)度fc、土壤孔隙度n、土壤有效孔隙度θe、濕潤鋒處土壤吸力Ψ(可以根據(jù)土壤性質(zhì)得到)、表層巖溶帶自由水蓄水容量SMEK(取值范圍10～50 mm)、表層巖溶帶自由水蓄水庫對表層巖溶帶側(cè)向出流的日出流系數(shù)KEKI、表層巖溶帶自由水蓄水庫對地下水的日出流系數(shù)KEKG、表層巖溶帶側(cè)向出流消退系數(shù)CEKI(取值范圍0.930～0.998)、巖溶區(qū)地下徑流消退系數(shù)CGK(取值范圍0.920～0.998)、表層巖溶帶側(cè)向出流馬斯京根法匯流參數(shù)KEEK、和表層巖溶帶側(cè)向出流馬斯京根法匯流參數(shù)XEEK(取值范圍0～0.5)。參數(shù)取值參考了新安江模型參數(shù)、GXAJ模型參數(shù)以及相關(guān)研究成果[4,20-21,31-32]。國內(nèi)外許多研究表明，表層巖溶帶的透水性隨著表層巖溶帶的深度增大而降低[2,33]，對于反映表層巖溶帶透水性的參數(shù)KEKI和KEKG，假設(shè)KEKI與KEKG之和同表層巖溶帶厚度呈線性關(guān)系，通過率定KEKI與KEKG之和的最大值與最小值，即可確定KEKI與KEKG之和的空間分布，假設(shè)KEKI與KEKG的比例為一個常數(shù)，通過率定該常數(shù)即可確定流域KEKI和KEKG的空間分布。馬斯京根法匯流參數(shù)參考了姚成[28]對GXAJ模型馬斯京根法匯流參數(shù)的研究。未介紹的參數(shù)可參考GXAJ模型[27，32，34]。

2 研究區(qū)概況與數(shù)據(jù)資料

2.1 研究區(qū)概況

以云南省文山壯族苗族自治州富寧縣范圍內(nèi)的普廳河為研究區(qū)，普廳河屬珠江流域西江水系右江二級支流，發(fā)源于富寧縣花甲鄉(xiāng)戈博村附近，于羅村口匯入右江。普廳河流域(23°30′N～23°50′N，105°20′E～105°40′E)是典型的滇東南巖溶高原與廣西盆地間的過渡地帶。富寧站多年平均徑流量1.66億m3(1984—2014年)，多年平均降水量 1 088 mm(1984—2014年)。普廳河流域為南亞熱帶季風(fēng)氣候，水汽主要來自印度洋孟加拉灣的西南暖濕氣流和太平洋南?；虮辈繛车臇|南暖濕氣流，具有明顯季節(jié)性，年平均氣溫19.3 ℃，最高月平均氣溫25.2 ℃，最低月平均氣溫10.8 ℃。流域內(nèi)建有清華洞水庫，該水庫為中型水庫，于2006年建成驗收。富寧水文站以上控制面積為373 km2，流域內(nèi)設(shè)置4處雨量站(圖4)。

圖4 普廳河流域概況Fig.4 General situation of Puting River Watershed

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

a.水文氣象資料。普廳河流域流量、降水量、蒸發(fā)量等水文氣象資料由云南省水文局提供，研究所用資料包括日尺度和小時尺度的站點降水、日尺度的站點蒸發(fā)、日尺度和小時尺度的富寧水文站流量資料，資料范圍為1984—2014年。普廳河流域于2006年修建了清華洞水庫，由于沒有水庫調(diào)度資料，并且流域面積較小，流域出口流量受水庫修建以及水庫蓄泄洪水的影響大，故選取了1984—2003年的數(shù)據(jù)用于模型模擬。

b.下墊面資料。表1列出了普廳河流域各類土壤所對應(yīng)的格林-安普特公式下滲參數(shù)值。GKDHM模型所需的新增流域下墊面信息為流域巖溶發(fā)育情況與石漠化情況。根據(jù)地質(zhì)云(http://geocloud.cgs.gov.cn/)平臺的《西南巖溶區(qū)石漠化分布圖》和《中國南部及東南亞地區(qū)巖溶發(fā)育程度圖》，對地質(zhì)圖進(jìn)行配準(zhǔn)提取和矢量化得到圖5。其余所需的下墊面信息與GXAJ模型相同，所有下墊面信息柵格尺度為 1 km×1 km。

表1 普廳河流域格林-安普特公式下滲參數(shù)值Table 1 Infiltration parameters for Green-Ampt infiltration equation of Puting River Watershed

(a) 流域巖溶分布

(b) 流域石漠化情況圖5 GKDHM模型新增下墊面信息Fig.5 Newly added underlying surface informationrequired by GKDHM

3 結(jié)果與討論

3.1 模擬結(jié)果比較

模型模擬結(jié)果采用徑流深合格率、洪峰合格率、確定性系數(shù)3個指標(biāo)進(jìn)行評估(表2)，合格標(biāo)準(zhǔn)參考GB/T 22482—2008《水文情報預(yù)報規(guī)范》。選取1984—2003年的20場洪水，其中前12場洪水用于率定模型參數(shù)，后8場洪水用于驗證模擬結(jié)果。

在徑流深模擬方面，GKDHM模型有較好的模擬效果，率定期合格率75.0%，驗證期合格率75.0%，相較于GXAJ模型率定期58.3%和驗證期50.0%的合格率有顯著提升。GKDHM模型率定期和驗證期平均相對誤差分別為15.5%和17.5%，對比GXAJ模型的18.6%和25.1%，GKDHM模型20場洪水的平均徑流深誤差降低了4.92%。

表2 普廳河流域模擬結(jié)果對比Table 2 Comparison of simulation results of Puting River Watershed

普廳河流域為山區(qū)中小流域，較大的洪峰易導(dǎo)致泥石流、滑坡等災(zāi)害的發(fā)生，能夠準(zhǔn)確模擬洪峰十分重要。對于洪峰模擬，GKDHM模型率定期合格率50.0%，平均相對誤差25.5%，驗證期合格率62.5%，平均相對誤差20.5%;GXAJ模型率定期合格率50.0%，平均相對誤差26.7%，驗證期合格率50%，平均相對誤差27.5%。GKDHM模型有效地提升了洪峰模擬合格率，20場洪水平均洪峰誤差降低了3.59%。

總體來看，GKDHM模型的模擬效果優(yōu)于GXAJ模型。率定期和驗證期的平均確定性系數(shù)，GKDHM模型的結(jié)果都高于GXAJ模型的結(jié)果。20場洪水中GKDHM模型模擬結(jié)果確定性系數(shù)大于0.5的場次為13場，GXAJ模型為11場。確定性系數(shù)的提升表明GKDHM模型能更好地模擬洪水過程。GKDHM模型分布式的結(jié)構(gòu)能更準(zhǔn)確地反映巖溶流域空間異質(zhì)性較大的下墊面情況對產(chǎn)匯流的影響。

表3 普廳河流域模型參數(shù)率定結(jié)果Table 3 Calibrated parameters of model for Puting River Watershed

3.2 模型參數(shù)分析

圖6 表層巖溶帶自由水蓄水容量空間分布Fig.6 Spatial distribution of free waterstorage capacity of epikarst

(a) KEKG (b) KEKI圖7 表層巖溶帶自由水蓄水庫參數(shù)空間分布Fig.7 Spatial distributions of parameters for freewater reservoirs of epikarst

(a) 1985090709號洪水(GKDHM模型) (b) 1985090709號洪水(GXAJ模型)

(c) 2002061209號洪水(GKDHM模型) (d) 2002061209號洪水(GXAJ模型)圖8 洪水流量過程與徑流成分Fig.8 Observed and simulated hydrographs and associated runoff components

根據(jù)圖5(a)，流域內(nèi)除一般的巖溶發(fā)育區(qū)外，還存在巖溶洞穴區(qū)域。基于式(4)，通過CKG1和CKG2兩個巖溶區(qū)地下水消退系數(shù)分別反映巖溶發(fā)育區(qū)和巖溶洞穴區(qū)域?qū)Φ叵滤恼{(diào)蓄作用。巖溶洞穴退水能力較強，地下水平均匯集時間短，地下水消退系數(shù)較小，最終率定的巖溶地下水消退系數(shù)CKG2小于CKG1。非巖溶區(qū)域的地下水消退系數(shù)為CG。

3.3 模擬流量徑流成分分析

圖8為GKDHM模型和GXAJ模型對兩場典型洪水的模擬結(jié)果,由圖8(b)(d)可以看出，GXAJ模型模擬流量的徑流成分以地表徑流和壤中流為主，地下徑流在整個洪水過程中所占比例很小，且退水緩慢。本流域有大范圍巖溶發(fā)育并且石漠化范圍較大，流域的徑流成分中地下徑流應(yīng)該占據(jù)一定的比重，且退水較快，因此GXAJ模型徑流成分模擬結(jié)果并不理想。由于引入格林-安普特下滲過程和表層巖溶帶自由水蓄水庫調(diào)蓄過程，圖8(a)(c)顯示，GKDHM模型模擬徑流成分以地表徑流和地下徑流為主，壤中流較少，因此GKDHM模型可以更好地刻畫巖溶流域土壤發(fā)育差、壤中流較少、地下徑流比例大和地下徑流退水快的特點。

4 結(jié) 語

a.模擬精度方面，GKDHM模型降低了徑流深誤差和洪峰誤差，在徑流深、流量過程方面的模擬精度較GXAJ模型有較大提高，洪峰流量模擬精度有一定的提升。GKDHM模型可以更好地反映巖溶流域表層巖溶帶以及深層地下溶隙的發(fā)育對于產(chǎn)匯流的影響，有效地提高了洪水模擬的精度。

b.模型結(jié)構(gòu)和參數(shù)方面，GKDHM模型通過構(gòu)建表層巖溶帶自由水蓄水庫和深層地下溶隙線性水庫，使模型參數(shù)的空間分布更能反映流域巖溶發(fā)育特性，模擬徑流成分更加符合巖溶流域水文特征。

c.實際應(yīng)用方面，本文提出的巖溶流域降水徑流模擬方法有效提高了普廳河流域洪水模擬精度，可為受表層巖溶帶以及深層地下溶隙影響的洪水預(yù)報提供方法與技術(shù)支撐，在我國西南中小河流巖溶流域有較大的推廣應(yīng)用價值。