考慮基巖裂隙作用的流域水文過程模擬

石 清，陳 喜，魏玲娜，張志才

(1.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇南京 210098;2.中國水電顧問集團北京勘測設(shè)計研究院，北京 100024;3.南京信息工程大學水文氣象學院，江蘇 南京 210044)

美國華盛頓大學開發(fā)的分布式水文-土壤-植被模型(DHSVM)包含了地表與土壤水流的描述，可以較好地模擬洪水過程［1］。然而，已有研究發(fā)現(xiàn)在進行長期水文過程模擬時，現(xiàn)有模型對枯期徑流擬合普遍較差［2］。國內(nèi)新安江模型用于流量過程模擬時，退水段的流量過程主要受地下水消退系數(shù)控制，該系數(shù)對退水段影響較大［3-4］。在參數(shù)調(diào)試過程中，退水段的擬合精度有時與整個流量過程的擬合精度相矛盾，故在模擬過程中需綜合考慮兩者之間的關(guān)系，盡量做到兼顧兩者的擬合精度［5］。

在我國南方廣大山丘區(qū)，枯期徑流不僅來源于土壤層，還來源于下覆基巖裂隙水?；鶐r裂隙水是我國分布最廣泛的地下水類型之一，淺部裂隙地下水都與上覆土壤水有不同程度的循環(huán)交替［6］，進而影響產(chǎn)匯流過程。筆者以東江流域內(nèi)星豐流域為研究區(qū)，基于DHSVM 2.0版本，根據(jù)流域地形、土壤、植被覆蓋率動態(tài)變化特征以及基巖裂隙水對徑流的調(diào)節(jié)作用，建立考慮植被-土壤-基巖裂隙作用的分布式水文模型，模擬紅壤丘陵區(qū)流域水文過程。

1 流域概況及觀測資料

星豐流域地處東江上游，為高丘陵山區(qū)，面積42.6 km2，海拔150～508 m，平均比降為3.14%。流域內(nèi)土壤主要為麻赤紅壤、頁赤紅壤和潴育水稻土，其中麻赤紅壤占72%，土層較厚。按國際制土壤分類標準，流域土壤概化為壤土、砂壤土和黏土3種類型。經(jīng)室內(nèi)實驗分析，土壤密度為1.32～1.45 g/cm3，孔隙率為41% ～45%，土壤含水率為12.3% ～32.7%。

星豐流域位于地下水二級功能區(qū)東江河源和平地下水水源涵養(yǎng)區(qū)，基巖裂隙水是枯期徑流的主要來源［7］。研究區(qū)內(nèi)基巖多屬變質(zhì)巖類裂隙含水巖組，石英砂巖占77.4%，粉砂巖占22.6%?；◢弾r風化后形成的風化殼比較深厚，一般在10 m以上，深厚者可達30～40 m。風化殼含有很多石英砂，其黏結(jié)力差、透水性強、結(jié)構(gòu)松散［8］。

星豐流域林地占流域面積的91%，耕地占8%，其余為建筑用地和水面。野外實地考察發(fā)現(xiàn)，星豐流域現(xiàn)狀覆被以松、杉混雜的針葉林為主?；?001年1月8日至2008年12月18日共184期250 m×250m空間分辨率的MODIS植被指數(shù)EVI數(shù)據(jù)(MOD13Q1)，利用Mu等［9］2007年提出的算法，反演不同類型的植被覆蓋率F(圖1(a))?；贛ODIS數(shù)據(jù)(MOD13Q1，MOD15A2，MOD09A1)，求得不同類型的植被覆蓋率F(圖1(a))、典型覆被葉面積指數(shù)(L)的月變化曲線(圖1(b))以及各種植被類型的短波冠層反射率A(圖1(c))。

圖1 不同植被的F，L和A月變化過程Fig.1 Monthly variations of fractional coverage F，leaf area index L，and canopy reflection A for different types of vegetation

采用流域內(nèi)4個雨量站和1個水文站2004—2008年實測逐日降雨和流量觀測資料，并移用研究區(qū)附近河源氣象站氣溫、風速、相對濕度、長短波輻射等氣象資料進行模型計算。

2 DHSVM及改進

2.1 DHSVM結(jié)構(gòu)

DHSVM以地面高程模型(DEM)柵格節(jié)點為中心，將流域劃分為若干計算柵格單元。根據(jù)柵格單元地形特征，描述短波輻射、降雨、氣溫、坡面流及柵格單元土壤和植被特性時空變化，在每一個計算時段內(nèi)，根據(jù)各柵格的能量和質(zhì)量平衡方程聯(lián)立解、柵格單元間坡面流和壤中流匯流演算方法，建立柵格單元水文聯(lián)系。根據(jù)河道分級演算，計算河網(wǎng)匯流過程。DHSVM實現(xiàn)了柵格尺度上流域水文過程及其與下墊面植被、土壤等影響因素的系統(tǒng)描述。

DHSVM中降雨分別被上、下冠層植被截留，直至最大截留能力Icj:

式中:rj——冠層截留系數(shù);Fj——冠層覆蓋率。下標j若為o則表示上冠層，若為u則表示下冠層。

冠層截留量變化為

式中:St——t時刻冠層截流量;P——降雨量;EIj——植被潮濕部分蒸發(fā)量(即截留雨量蒸發(fā))。

采用雙層Penman-Monteith方法模擬上、下冠層散發(fā)量，結(jié)合土壤水分條件計算土壤實際蒸發(fā)量［10］，三部分之和作為實際總蒸散發(fā)量。

a.植被潮濕部分蒸發(fā)量:

其中

式中:Epj——潛在蒸散發(fā)量;Awj——潮濕部分面積，干燥部分則為(1-Awj);Δtw——以Epj蒸發(fā)截留雨量所需時間。

b.植被干燥部分散發(fā)量:

式中:Ej——植被表面蒸騰量;ζ——飽和水汽壓～氣溫關(guān)系曲線的斜率;γ——空氣濕度常數(shù);rcj——冠層阻力;raj——空氣動力學阻力;Δt——計算時段長。

c.土壤實際蒸發(fā)量:

其中

式中Fe為土壤解吸率，由土壤對表層的輸水速率確定，是土壤類型和上層土壤含水量的函數(shù)。d.土壤層水量平衡采用3層根帶模式模擬水分在上、中、下非飽和土壤層中的運動:

式中:dk——第k層土壤的厚度;n——土壤層數(shù);θk——第k層土壤的含水率;If——計算時段內(nèi)滲入土壤層中的水量;Qv——流入下層的水量;fo，k，fu，k——第k層土壤中對上冠層、下冠層的根系比;Vexk——地下水位上升時地下水對土壤水的補給;Qsin，t，QS，t——計算時段初側(cè)向壤中流的入流量、出流量。

假定土壤飽和滲透系數(shù)隨土壤深度呈指數(shù)遞減，推導出土壤導水系數(shù)Tx，y的計算公式如下:

式中:Kx，y——柵格(x，y)表層土壤側(cè)向飽和滲透系數(shù);fx，y——指數(shù)衰減系數(shù);zx，y——地下水位埋深;Dx，y——土壤層厚度。

采用運動波方法計算側(cè)向飽和壤中流，用Wigmosta［11］發(fā)展的準三維路徑模式計算坡面流。坡面地表流、壤中流逐柵格匯流并流向流域低洼處匯集進入河道，利用線性槽蓄法進行河道匯流演算。

2.2 模型改進

分析星豐流域2003—2008年實測日流量過程，依據(jù)飽和水運動Depuit-Boussinesq方程對退水過程較長的典型年份(2003年、2004年、2006年)退水過程進行擬合(圖2)。由圖2可知，退水特征明顯表現(xiàn)為2個階段，原因在于南方紅壤丘陵區(qū)多為二元水文地質(zhì)結(jié)構(gòu)［12］，退水段徑流主要來源于土壤層及基巖裂隙水。上層土壤層滲透性按式(9)進行衰減，退水較快;下層基巖裂隙退水過程較穩(wěn)定。因此，筆者在飽和土壤層之下添加裂隙含水層模塊，綜合考慮土壤水和基巖裂隙水之間的水力聯(lián)系，對DHSVM進行改進。

采用Wigmosta等［13］提出的地形驅(qū)動影響下飽和壤中流的計算方法，利用地下水位(而不是地面高程)計算柵格之間的水力梯度。由于水頭的變化，需要對每一個計算時段內(nèi)的地下水位和相應(yīng)的坡度、坡向、流向、水力梯度進行計算。裂隙含水層的水分運動分為垂向和側(cè)向兩部分進行計算，地下水流量的計算公式

圖2 星豐流域枯期徑流退水過程規(guī)律Fig.2 Low-flow curve fitted by Depuit-Boussinesq model for Xingfeng Basin

如下:

式中:fi——柵格(x，y)在i方向的流線寬度;GGW——水力梯度;Kl——裂隙含水層側(cè)向水力傳導率;S——柵格單元內(nèi)含水層含水量;Kd——含水層下基巖的垂向傳導率。

裂隙含水層的水量平衡方程為

式中:St，St+Δt——t，t+ Δt時刻裂隙含水層含水量;(QGin，t-QG，t)——裂隙含水層中基流的凈增水量(單元柵格中側(cè)向流入量QGin，t減去流出量QG，t);Qv，GW——上覆土壤層向裂隙含水層滲漏的水量;Vex，GW——裂隙含水層水位上升對其上覆土壤層的補給量。

緊鄰裂隙層的上覆土壤層水量平衡方程式改寫為

當裂隙含水層沒有到達飽和時，上覆飽和壤中流以下滲的形式對裂隙含水層進行補給，直至裂隙含水層飽和。當柵格單元內(nèi)裂隙含水層的含水量大于其最大容量時，多余的水量將“返還”給上覆土壤層。除了和上覆土壤進行水分交換，裂隙含水層水分還會向下滲透，逸出含水層，形成“深層損失”，“深層損失率”由基巖的水力傳導率決定。當柵格單元內(nèi)的地下含水層水位高于河床，裂隙含水層水量將會被河道截留。

3 模型應(yīng)用

3.1 模型率定與驗證

根據(jù)DEM將流域劃分為50 m×50 m柵格，共162行、156列。模型計算參數(shù)率定期為2003—2006年，驗證期為2007—2008年，計算步長為2 h。根據(jù)模型計算與實測流域出口斷面流量過程，采用人工試錯法對模型參數(shù)進行率定。通過確定性系數(shù)和徑流深相對誤差來評測模型日徑流模擬的精度。表1、表2列出了模型主要土壤、基巖裂隙以及植被特征參數(shù)的取值。

表1 土壤及基巖裂隙特征參數(shù)Table 1 Soil property parameters

表3為利用DHSVM及本文改進的DHSVM進行參數(shù)率定與模型驗證的結(jié)果。如不考慮基巖裂隙水，率定期枯水年份(2003年、2004年)計算徑流量R顯著偏小，豐水年份(2005年、2006年)則偏大。這是因為DHSVM本身在枯水處理方面有一定缺陷，存在退水衰減過快的問題，枯季徑流模擬效果不佳。在模型率定中，由于枯水年徑流較小，模擬徑流深偏小，為了在連續(xù)模擬中保持水量平衡，使豐水年徑流模擬結(jié)果相應(yīng)地偏高。進一步分析發(fā)現(xiàn)，裂隙含水層對年內(nèi)徑流分配具有一定的調(diào)節(jié)作用，改進后DHSVM模擬年徑流確定性系數(shù)有顯著提高(表3)，流量過程線的吻合度較高(圖3)。

表2 植被特征參數(shù)Table 2 Vegetation property parameters

表3 DHSVM改進前后蒸散發(fā)及徑流模擬結(jié)果Table 3 Comparison of evapotranspiration and streamflow before and after DHSVM revision

圖3 實測與模擬流量過程線Fig.3 Observed and simulated streamflow

利用實測流量和模型改進前、后的模擬流量序列，采用期望值公式，將流量按從大到小進行排頻(圖4)，顯然，考慮裂隙作用對枯期流量模擬結(jié)果改進顯著。進一步選取頻率為75%的流量值Q75%作為枯季徑流的特征值，統(tǒng)計出不同歷時下實測和模型改進前、后模擬的枯季流量低于Q75%的發(fā)生次數(shù)(圖5)。通過比較可以發(fā)現(xiàn)，改進后枯期徑流模擬精度比改進前有較大提高。

3.2 參數(shù)靈敏度分析

選用擾動分析法，計算參數(shù)在小范圍內(nèi)變化所導致的模型輸出變化率，通過敏感度因子Sen加以衡量，計算公式如下:

式中:y，y'——參數(shù)變化前、后輸出值;x——變化前參數(shù)值;Δ——參數(shù)變化值。

選取各參數(shù)敏感度因子的算數(shù)平均值作為該參數(shù)的靈敏度值，模型計算輸出值選用上覆土壤層含水量、日平均蒸散發(fā)量、日平均流量以及汛期日平均流量、枯季日平均流量。

圖4 實測與模擬流量歷時曲線Fig.4 Observed and simulated flow duration curves

分析表明:對上覆土壤層含水量最敏感的參數(shù)是孔隙度Φ和田間持水量θfc，其他參數(shù)敏感度都較小;對日平均蒸散發(fā)量較敏感的參數(shù)有植被參考高度H，F(xiàn)，θfc、最小氣孔阻力Rsmin;對日平均徑流量最敏感的參數(shù)排序為H，F(xiàn)，θfc，Rsmin，其中汛期日徑流量敏感性系數(shù)較大的參數(shù)主要有H，F(xiàn)，Φ和裂隙含水層孔隙度GP，而枯期日徑流量敏感性系數(shù)較大的參數(shù)主要是Φ，θfc，GP。由此可見，上覆土壤含水量及蒸散發(fā)主要受土壤、植被參數(shù)影響，基巖裂隙含水層參數(shù)對汛期和枯期徑流起調(diào)節(jié)作用，減少汛期流量、增加枯期流量，對枯期流量影響顯著，但對全年流量影響不大。

圖5 各時段實測和模擬流量低于Q75%的數(shù)目Fig.5 Occurrences of observed and simulated flow below Q75%

4 結(jié) 語

根據(jù)東江星豐流域植被、土壤、基巖特征，在考慮裂隙含水層對徑流的調(diào)節(jié)功能情況下對DHSVM進行改進。流域水文過程模擬及參數(shù)靈敏性分析表明，基巖裂隙水對徑流過程模擬效果影響顯著，模型改進后對徑流過程的模擬精度有較大提升，尤其對枯期徑流的模擬更接近實際觀測結(jié)果?；鶐r裂隙含水層參數(shù)對汛期和枯期徑流起調(diào)節(jié)作用，減少汛期流量、增加枯期流量，對枯期流量影響尤為顯著。

［1］WIGMOSTA M S，VAIL L W，LETTENMAIER D P.A distributed hydrology-vegetation model for complex terrain［J］.Water Resources Research，1994，30(6):1665-1679.

［2］BATHURST JC.Physically-based distributed modelling of an upland catchment using the systeme hydrologique europeen［J］.Journal of Hydrology，1986，87(1/2):79-102.

［3］ZHAO Renjun.The Xin＇anjiang model applied in China［J］.Journal of Hydrology，1992，135(1/2/3/4):371-381.

［4］趙人俊.流域水文模擬:新安江模型與陜北模型［M］.北京:水利電力出版社，1984.

［5］黃國如，陳永勤.基于新安江水文模型的東江流域枯水徑流模擬［J］.華南理工大學學報:自然科學版，2006，34(11):93-98.(HUANG Guoru，CHEN Yongqin.Low flow simulation of the Dongjiang Basin based on Xin ＇anjiang watershed hydrological model［J］.Journal of South China University of Technology:Natural Science Edition，2006，34(11):93-98.(in Chinese))

［6］陳宇，溫忠輝，束龍倉.基巖裂隙水研究現(xiàn)狀與展望［J］.水電能源科學，2010，28(4):62-65.(CHEN Yu，WEN Zhonghui，SHU Longcang.Status and prospect of research on bedrocKfissure water［J］.Water Resources and Power，2010，28(4):62-65)(in Chinese))

［7］廣東省水利廳.廣東省地下水功能區(qū)劃［R］.廣州:廣東省水利廳，2009.

［8］鄧南榮，李定強，王繼增，等.廣東省東江流域土壤侵蝕空間分布特征研究［J］.中國水土保持，1999(5):21-23.(DENGNanrong，LIDingqiang，WANGJizeng，et al.The characteristic of spatial distribution of soil erosion in Dongjiang Basin of Guangdong Province［J］.Soil and Water Conservation in China，1999(5):21-23.(in Chinese))

［9］MU Q，HEINSCH F A，ZHAO M，et al.Development of a global evapotranspiration algorithm based on MODIS and global meteorology data［J］.Remote Sensing of Environment，2007，11(4):519-536.

［10］EAGLESON P.Climate，soil，and vegetation:3.A simplified model of soil moisture movement in the liquid phase［J］.Water Resources Research，1978，14(5):722-730.

［11］WIGMOSTA M S，PERKINS W A.Simulating the effects of forest roads on watershed hydrology［J］.Water Science and Application，2001，2:127-143.

［12］劉昆，何捷，焦樹林，等.粵東花崗巖流域覆蓋型風化殼風化特征比較［J］.生態(tài)環(huán)境，2008，17(5):1937-1941.(LIU Kun，HE Jie，JIAO Shulin，et al.Comparison on the weathering characteristics of the alluvium-covered weathering crust for the granite catchments in eastern Guangdong province［J］.Ecology and Environment，2008，17(5):1937-1941.(in Chinese))

［13］WIGMOSTA M S，LETTENMAIER D P.A comparison of simplified methods for routing topographically driven subsurface flow［J］.Water Resources Research，1999，35(1):255-264.