梁桂星，覃小群，崔亞莉，陳 爽，黃奇波

(1.中國地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083；2.北京市地?zé)嵫芯吭?，北?102218；3.中國地質(zhì)科學(xué)院巖溶地質(zhì)研究所，廣西 桂林 541004)

在我國西南地區(qū)，巖溶出露面積約為114×104km2，約占中國總巖溶面積的1/3。該地區(qū)雖然雨量豐富，但土層較薄，石漠化嚴(yán)重，且?guī)r溶構(gòu)造發(fā)育，含水層垂向上水力聯(lián)系密切，水文過程變化迅速，儲(chǔ)水能力低，造成該地區(qū)水資源環(huán)境脆弱[1-2]。巖溶水作為西南地區(qū)主要分布的水資源，其有效利用關(guān)系著國計(jì)民生，因此亟需查明該地區(qū)水文循環(huán)，科學(xué)合理地對(duì)該類型地區(qū)水資源進(jìn)行評(píng)價(jià)。

國內(nèi)外有眾多學(xué)者在巖溶地區(qū)探索應(yīng)用水文模型。YU L L等[3]在美國佛羅里達(dá)巖溶地區(qū)運(yùn)用MODFLOW-DCM模型模擬地下水流，并分析了含水層介質(zhì)和補(bǔ)給條件對(duì)巖溶水動(dòng)力特性的影響；楊鄭秋[4]在湖南香花嶺地區(qū)應(yīng)用雙系統(tǒng)耦合模型MODFLOW-CFP模擬水文過程，并對(duì)該地區(qū)地下水資源進(jìn)行評(píng)價(jià)；劉建國等[5]以貴州六枝特區(qū)某項(xiàng)目為例進(jìn)行地下水?dāng)?shù)值模擬研究，探索了在FEFLOW離散特性功能模塊下巖溶管道流模型的刻畫方式。這類水文模型雖在巖溶地區(qū)得到較好的應(yīng)用，但多數(shù)是基于地下水系統(tǒng)對(duì)含水層空間結(jié)構(gòu)進(jìn)行概化，且未能細(xì)致刻畫地表水與地下水的相互轉(zhuǎn)化關(guān)系。為了綜合考慮流域內(nèi)土壤質(zhì)地、地形地貌、土地利用等地理要素對(duì)水文效應(yīng)的影響，且合理刻畫水文循環(huán)過程，本文選用分布式水文模型SWAT(Soil and Water Assessment Tool)[6]。

SWAT模型由美國農(nóng)業(yè)部研制開發(fā)，已經(jīng)被廣泛運(yùn)用在不同的土地利用、土壤類型等復(fù)雜地理特征要素條件下的流域內(nèi)地表水量、泥沙及水質(zhì)的評(píng)價(jià)，且取得較好的效果[7-8]。國內(nèi)外在巖溶地區(qū)應(yīng)用SWAT模型的研究較多，但多數(shù)未能充分考慮巖溶特性，有些研究雖然考慮了巖溶地區(qū)特殊的含水層空間結(jié)構(gòu)，但很少細(xì)化巖溶區(qū)表層巖溶帶對(duì)水分運(yùn)移的影響。表層巖溶帶是由于強(qiáng)烈的巖溶化過程，表層碳酸鹽巖形成各種犬牙交錯(cuò)的巖溶個(gè)體形態(tài)和微形態(tài)組合構(gòu)成不規(guī)則帶狀的強(qiáng)巖溶化層[9]，其主要特點(diǎn)是具有賦存和調(diào)蓄地下水的功能，表現(xiàn)在兩個(gè)方面：一是巖溶水量上的調(diào)蓄，增加入滲補(bǔ)給量；二是巖溶水流過程的調(diào)蓄，即通過復(fù)雜的表層巖溶系統(tǒng)延緩降雨后雨水在巖溶水系統(tǒng)停留的時(shí)間[10]。為了體現(xiàn)表層巖溶帶的特點(diǎn)，國內(nèi)學(xué)者對(duì)SWAT進(jìn)行改進(jìn)以近似代替表層巖溶帶對(duì)水文循環(huán)的影響。代俊峰等[2]對(duì)SWAT模型進(jìn)行改進(jìn)，結(jié)合南方丘陵特點(diǎn)建立了湖北省漳河灌區(qū)水文模型；WANG Y等[11]在SWAT中引入非線性地下水算法來削減巖溶地區(qū)由于寬大的巖溶裂隙對(duì)地下水流動(dòng)造成的影響；任啟偉[12]考慮巖溶特點(diǎn)改進(jìn)SWAT模型，對(duì)廣西刁江流域地表徑流進(jìn)行模擬并對(duì)流域水量進(jìn)行評(píng)價(jià)；陳爽[13]通過改進(jìn)SWAT模型地下水模塊，對(duì)左江崇左地區(qū)巖溶區(qū)和碎屑巖區(qū)地下水分別進(jìn)行評(píng)價(jià)。

目前，部分研究通過改進(jìn)SWAT模型來體現(xiàn)表層巖溶帶的調(diào)蓄特點(diǎn)，但是根據(jù)巖溶特性刻畫水分從土壤向地下水運(yùn)移過程的研究較少。本文將分別考慮上層滯水與表層巖溶水的相似性及研究區(qū)不同巖性分布，通過改進(jìn)SWAT模型土壤水模塊及地下水模塊，建立適用于漓江流域的分布式水文模型，為漓江流域水資源合理開發(fā)利用提供技術(shù)支撐。

1 研究方法

表層巖溶帶的概念是20世紀(jì)70年代由法國學(xué)者M(jìn)ANGIN A提出的，是指巖層表層具有大量溶蝕擴(kuò)大裂隙的強(qiáng)巖溶化帶。西南巖溶地區(qū)巖溶在垂直空間上發(fā)育存在差異性，勞文科等[14]在研究湖南洛塔地區(qū)時(shí)發(fā)現(xiàn)表層巖溶帶與下層飽水帶之間，存在一定厚度的巖溶弱發(fā)育帶，使得整個(gè)巖溶剖面具有雙層水位的現(xiàn)象。這種雙層水位的現(xiàn)象與上層滯水有很大相似性，上層滯水經(jīng)常分布于砂層中的黏土夾層之上和石灰?guī)r中溶洞底部有黏性土充填的部位，由于局部隔水層的存在而出現(xiàn)雙層水位現(xiàn)象。上層滯水的排泄主要是蒸發(fā)及沿隔水層邊緣下滲；由于下層巖溶發(fā)育性較差且?guī)r體連貫性較好，表層巖溶帶主要是沿著下層弱透水層以泉的形式排泄及通過弱透水層向下層包氣帶滲漏。

適合SWAT模型的土壤水量平衡方程為：

SWi=SWi-1+Rday,i-Qlat,i-Ea,i-Wrchrg,sh

(1)

式中：SWi——第i天土壤含水量/mm；

SWi-1——第i-1天土壤含水量/mm；

Rday,i——第i天滲入土壤中水量/mm；

Qlat,i——第i天土壤剖面補(bǔ)給河流的量/mm；

Ea,i——第i天土壤水蒸發(fā)量/mm；

Wrchrg,sh——第i天土壤底部滲透到淺層含水層的量/mm。

淺層含水層的水量均衡方程為：

aqsh,i=aqsh,i-1+Wrchrg,sh-Qgw-Wrevap-Wpump,sh

(2)

式中：aqsh,i——第i天淺層含水層的儲(chǔ)水量/mm；

aqsh,i-1——第i-1天淺層含水層的儲(chǔ)水量/mm；

Wrchrg,sh——第i天淺層含水層的補(bǔ)給量/mm；

Qgw——第i天匯入主河道的地下水量/mm；

Wrevap——第i天因土壤水分不足而進(jìn)入土壤帶的量/mm；

Wpump,sh——第i天淺層含水層的抽水量/mm。

巖溶地區(qū)產(chǎn)匯流過程有兩種形式，分別為降雨通過較為寬大的巖溶裂隙及管道形成快速徑流和通過細(xì)小裂隙形成的慢速徑流。SWAT模型對(duì)土壤中上層滯水有較為詳細(xì)的刻畫，但是不能體現(xiàn)由于巖溶發(fā)育差異而存在的不同狀態(tài)的裂隙流。為了體現(xiàn)巖溶地區(qū)表層巖溶帶的調(diào)蓄能力和水分在不同裂隙狀態(tài)下的運(yùn)動(dòng)，需要對(duì)SWAT模型土壤水模塊及地下水模塊進(jìn)行改進(jìn)。

首先刻畫水分在表層巖溶帶中的運(yùn)動(dòng)：

Wst,i=Wst,i-1+Rday-Wloss

(3)

式中：Wst,i——第i天表層巖溶帶水量/mm；

Wst,i-1——第i-1天表層巖溶帶水量/mm；

Rday——進(jìn)入表層巖溶帶的水量/mm；

Wloss——表層巖溶帶排泄的水量/mm。

其中表層巖溶帶排泄的水量可表示為：

Wloss=Wst-FC

(4)

式中：Wst——表層巖溶帶含水量/mm；

FC——表層巖溶帶最大可持水量/mm。

然后對(duì)水分從表層巖溶帶進(jìn)入淺層含水層過程進(jìn)行刻畫，可用式(5)表示：

Wrchrg,sh=Qfast+QSlow

(5)

式中：Qfast——以快速流進(jìn)入淺層含水層的量/mm；

QSlow——以慢速流進(jìn)入淺層含水層的量/mm。

Qfast=Vcrk

(6)

式中：Vcrk——寬大裂隙體積/mm。

表層巖溶帶的慢速裂隙流過程與變水頭的水柱過程相似，從表層巖溶帶滲入淺層含水層的水分存在滯后效應(yīng)。

(7)

式中：ΔT——模擬步長/24 h；

TTperc——衰減系數(shù)。

各層的衰減系數(shù)不同，計(jì)算公式為：

(8)

式中：SAT——表層巖溶帶的飽和含水量/mm；

K——表層巖溶帶飽和水力傳導(dǎo)率/(mm·h-1)。

2 模型應(yīng)用

2.1 區(qū)域概況

研究區(qū)位于廣西桂林市境內(nèi)漓江流域(圖1)。漓江發(fā)源于桂林市興安縣西北部，流域面積為10 889.13 km2，其中碳酸鹽巖區(qū)面積為4 379.27 km2，碎屑巖區(qū)面積為6 509.86 km2。漓江流域年均降水量為1 543.2 mm，年均氣溫為18.7 ℃，年均相對(duì)濕度75.6%。

圖1 研究區(qū)概況Fig.1 General situation of the study area

區(qū)內(nèi)含水巖組主要由融縣組灰?guī)r組成，根據(jù)含水組巖性及地下水動(dòng)力特征，研究區(qū)劃分出巖溶水、基巖裂隙水及孔隙水三大類。碎屑巖區(qū)以碎屑巖裂隙水為主，碳酸鹽巖區(qū)主要以碳酸鹽巖裂隙溶洞水及碳酸鹽巖溶洞裂隙水為主，孔隙水主要分在桂林市漓江兩岸。漓江河水天然狀況下一般不存在倒灌補(bǔ)給作用，無論枯水期或豐水期，地下水位均高于河水位。研究區(qū)內(nèi)雖然巖溶發(fā)育，但是大規(guī)模的地下河流較少，多以細(xì)小地下河分支存在于各個(gè)局部區(qū)域上。由于地下河與地表河流的相似性，可以統(tǒng)一劃為地表河流。

2.2 基礎(chǔ)數(shù)據(jù)

建立漓江流域SWAT水文模型需要的基礎(chǔ)地理數(shù)據(jù)主要包括DEM、土壤類型、土地利用類型、氣象、徑流、水庫、水系及巖性分布等數(shù)據(jù)。該研究區(qū)涉及5個(gè)水庫，除了收集到青獅潭水庫和思安江水庫月出流數(shù)據(jù)外，其他3個(gè)水庫均未收集到出流數(shù)據(jù)，僅有水庫屬性參數(shù)。為了盡可能反映水文地質(zhì)條件，本文在無水庫天然條件下建立水文模型，根據(jù)模型在天然條件下模擬的流域出流量作為水庫的出流量，使漓江流域水文模型更加完善。水庫及水文站位置分布見圖2(a)。在模型應(yīng)用中，利用GIS通過從DEM和土地利用類型中獲取地理特征參數(shù)，且通過GIS鏈接土壤數(shù)據(jù)庫，從而獲取土壤質(zhì)地的空間分布及結(jié)構(gòu)參數(shù)，其中主要包括土壤寬裂隙體積比、飽和水力傳導(dǎo)率、持水率、顆粒含量等，土壤類型分布見圖2(c)。SWAT模型雖然涉及的參數(shù)數(shù)據(jù)很多，但不屬于參數(shù)化模型，需要率定的參數(shù)很少。

2.3 模型建立

模型輸入的數(shù)據(jù)坐標(biāo)系統(tǒng)格式均為Beijing_1954_GK_Zone_19N投影坐標(biāo)系統(tǒng)。通過考慮巖性分布情況及水文站分布位置，對(duì)研究區(qū)子流域及水文響應(yīng)單元進(jìn)行劃分，為不同巖性流域水文分析做鋪墊。其中，通過導(dǎo)入矢量水系及流域內(nèi)水文觀測站，選擇流域出口，劃分了138個(gè)子流域。該流域碳酸鹽巖與碎屑巖呈交叉分布，巖性不同含水層空間結(jié)構(gòu)就有較大差異性；高度巖溶化的碳酸鹽巖地區(qū)巖溶構(gòu)造分布密集，水力聯(lián)系密切，多以捷徑式補(bǔ)給含水層；碎屑巖區(qū)垂向上土壤剖面水力聯(lián)系性差，多以活塞式推進(jìn)方式補(bǔ)給含水層。巖性不同的流域水文參數(shù)及水分運(yùn)移過程不同，因此需要針對(duì)不同巖性進(jìn)行分析驗(yàn)證。根據(jù)研究區(qū)巖性分布確定了碎屑巖區(qū)子流域56個(gè)，巖溶區(qū)子流域82個(gè)。研究區(qū)根據(jù)土地利用類型、土壤種類及坡度等地理要素總共劃分了875個(gè)水文響應(yīng)單元(HRU)，其中碎屑巖區(qū)265個(gè)HRU，巖溶區(qū)610個(gè)HRU。圖2(d)展示了模型通過輸入DEM數(shù)據(jù)、土地利用數(shù)據(jù)、土壤類型數(shù)據(jù)、坡度等級(jí)數(shù)據(jù)、水文水系分布數(shù)據(jù)劃分的水文響應(yīng)單元的空間分布。

3 結(jié)果與討論

3.1 敏感性分析及參數(shù)率定

巖性不同其水文地質(zhì)參數(shù)也有較大差異性，通過巖性劃分的子流域參數(shù)也有所不同，因此需要分別對(duì)不同巖性的子流域進(jìn)行參數(shù)校準(zhǔn)與率定后，才能建立與實(shí)際水文、水文地質(zhì)條件相符的數(shù)值模型。

SWAT模型涉及水文循環(huán)的流域參數(shù)很多，本次研究區(qū)選擇15個(gè)常規(guī)的流域參數(shù)進(jìn)行敏感性分析，利用SWAT-CUP中SUFI-2模塊對(duì)參數(shù)進(jìn)行1 000次迭代運(yùn)算，確定了巖溶區(qū)7個(gè)敏感性參數(shù)和碎屑巖區(qū)10個(gè)敏感性參數(shù)，各區(qū)域敏感性參數(shù)見表2和表3。GW_DELAY(地下水滯后系數(shù))：通過滲透或側(cè)向流穿過表層巖溶帶及土壤剖面最底層進(jìn)入淺層含水層的滯后時(shí)間，取決于潛水面的埋深及包氣帶中地層巖性特征；CN2(SCS徑流曲線數(shù))：SCS曲線數(shù)控制著地表徑流，是土壤滲透性、土地利用和前期土壤水分條件的函數(shù)；GW_REVAP(地下水再蒸發(fā)系數(shù))和RCHRG_DP(含水層滲漏系數(shù))分別主要取決于潛水面埋深及深層含水層的地層巖性，分別控制著淺層含水層的水量及深層含水層的水量。

利用建立的分布式水文模型對(duì)漓江流域2008—2016年逐月徑流過程進(jìn)行初步模擬。其中2008—2013年作為模型的參數(shù)率定期，2014—2016年作為模型參數(shù)的驗(yàn)證期。引用NS(納什系數(shù))和R2(決定系數(shù))來評(píng)價(jià)模型徑流與實(shí)測徑流的擬合效果，NS及R2大于0.5既滿足模型模擬要求。選取碎屑巖區(qū)靈渠(三)站及巖溶區(qū)桂林(三)站、潮田站、陽朔站等4個(gè)水文站對(duì)模型進(jìn)行識(shí)別驗(yàn)證，模型擬合結(jié)果見圖3，月徑流擬合評(píng)價(jià)結(jié)果見表4。

研究區(qū)土地利用類型以林地和耕地為主，分別占總面積的63.97%和22.31%；居民區(qū)主要分布在桂林市區(qū)周圍，所占比例較小；下墊面變化幅度不大對(duì)研究區(qū)產(chǎn)匯流影響較小。從圖3可以看出，整體上SWAT模型對(duì)漓江流域的徑流模擬效果較好，各水文站月模擬流量與實(shí)測流量趨勢擬合基本一致，符合流域的實(shí)際情況，但對(duì)部分較大洪峰模擬效果不太理想。究其原因主要有兩點(diǎn)：一是利用天然徑流量代替無水庫出流數(shù)據(jù)的水庫出流量，存在一定誤差；二是通過較為稀疏的國家氣象站插值的CMADS降水?dāng)?shù)據(jù)精度略低于局地降水?dāng)?shù)據(jù)精度。從具體水文站徑流模擬可以看出，桂林(三)、潮田站及陽朔站的R2均大于0.85，NS均大于0.75，整體徑流模擬效果較好，但是在枯水期模擬徑流略低于實(shí)測徑流。靈渠(三)站模擬效果較差，在枯水期模擬值大于觀測值，R2和NS均小于其他三個(gè)水文站的值。原因有二：一方面靈渠縣建有人工水庫，人為開采強(qiáng)烈，豐水期水量豐富，人為開采對(duì)徑流擬合的影響遠(yuǎn)小于枯水期；另一方面由于靈渠(三)站處于支流上，徑流較小，豐水期容易造成流量的峰值缺失。

表2 漓江流域模型參數(shù)敏感性分析Table 2 Sensitivity analysis of parameters of the Lijiang River Basin model

注：SUFI-2算法通過t檢驗(yàn)和p因子方法來判斷參數(shù)的敏感性，其中t的絕對(duì)值越大參數(shù)越敏感；p值表示t值的顯著性，p因子越小參數(shù)敏感性越顯著。

圖3 漓江流域各水文站模型改進(jìn)前后月徑流模擬比較Fig.3 Comparison of monthly runoff simulation before and after model improvement in the Lijiang River Basin

表3 模型參數(shù)率定結(jié)果Table 3 Calibration results of the model parameters

3.2 改進(jìn)前后對(duì)比

利用修改后的SWAT模型對(duì)漓江流域2008—2016年進(jìn)行模擬，并與改進(jìn)前和實(shí)測月徑流進(jìn)行對(duì)比，結(jié)果見圖3。改進(jìn)后的模型調(diào)整了枯水期的徑流，相對(duì)改進(jìn)前徑流模擬值有明顯的增加，與實(shí)測值吻合度較為理想。改進(jìn)后各水文站徑流擬合評(píng)價(jià)見表4，改進(jìn)后決定系數(shù)R2和納什系數(shù)NS有不同程度的提高，模擬的精度得到較大提高。由于流量的繼承性，碎屑巖區(qū)產(chǎn)生的徑流量對(duì)水文站徑流量的貢獻(xiàn)率(流量貢獻(xiàn)率是指該站上游碎屑巖地區(qū)產(chǎn)生的徑流量占該水文站總徑流量的比例)在50%以上，且?guī)r溶區(qū)地表徑流也占河道總徑流量一定比例，使得基流產(chǎn)生的滯后性在總徑流上反映得不明顯。

為反映表層巖溶的存在造成地下水補(bǔ)給的滯后性，繪制了各水文站所在子流域基流總量過程線(圖4)。由圖可見，位于碎屑巖區(qū)靈渠(三)站所在子流域改進(jìn)前后基流沒有發(fā)生明顯變化，而處于巖溶子區(qū)水文站子流域的基流量變化曲線均反映了不同程度的滯后性和消峰作用，表明改進(jìn)后的模型較好地刻畫了表層巖溶帶調(diào)蓄作用，提高了模型的精度。

表4 模型改進(jìn)前后月徑流過程模擬參數(shù)評(píng)價(jià)Table 4 Evaluation of monthly runoff process simulation parameters before and after model improvement

圖4 模型改進(jìn)前及改進(jìn)后基流對(duì)比Fig.4 Basic flow comparison before and after model improvement

4 結(jié)論

(1)在巖溶地區(qū)使用SWAT模型需要對(duì)其進(jìn)行改進(jìn)。本文添加了表層巖溶帶模塊，分別刻畫了水分以快速流與慢速流進(jìn)入淺層地下水的過程，實(shí)現(xiàn)了巖溶區(qū)復(fù)雜的水文地質(zhì)結(jié)構(gòu)水循環(huán)模擬，體現(xiàn)了表層巖溶帶的調(diào)蓄作用，符合巖溶流域水循環(huán)的特點(diǎn)。

(2)利用改進(jìn)后的模型對(duì)漓江巖溶流域進(jìn)行了模擬，結(jié)果表明：改進(jìn)后的模型與改進(jìn)前的相比較，徑流模擬精度顯著提升，尤其是枯水期徑流模擬得到明顯改善。

(3)針對(duì)不同巖性及空間結(jié)構(gòu)差異性，通過分析基流的滯后性，刻畫出巖溶區(qū)表層巖溶帶的存在對(duì)地下水補(bǔ)給的影響，可服務(wù)于巖溶區(qū)地下水資源評(píng)價(jià)。