羅智鋒 王 文 陳 喜
(1.河海大學水文水資源與水利工程科學國家重點實驗室,南京 210098;2.河海大學水文水資源學院,南京 210098)
水文尺度問題自20世紀90年代初被正式提出后,就受到廣泛關(guān)注和重視.水文科學的理論研究與實踐證明,不同尺度的水文循環(huán)機理有很大差異,如何考慮流域水文過程的時空不均勻性是尺度問題的關(guān)鍵[1].從空間上來說,影響這種不均勻性的主要因素有流域地形、土壤、植被等數(shù)據(jù)質(zhì)量.不同來源、不同網(wǎng)格分辨率的數(shù)據(jù)對流域特征描述的不一致,直接影響著模型的模擬精度.以往的研究表明,植被土壤數(shù)據(jù)的分辨率選擇對水文模擬影響較?。?],多數(shù)研究集中于分析DEM 網(wǎng)格水平分辨率的影響.DEM 是描述地表形狀的連續(xù)變化的數(shù)據(jù)源,用于提取坡度、坡向、匯流路徑長、集水面積、流域邊界等流域特征.DEM 網(wǎng)格分辨率深深影響著以地形為基礎的分布式水文模型.一些研究利用不同水平分辨率的DEM 提取地形特征,發(fā)現(xiàn)低分辨率導致平均坡度減小,流域面積和地形指數(shù)增加[3-4].不少研究者分析了Topmodel的尺度效應,發(fā)現(xiàn)其對網(wǎng)格分辨率變化非常敏感[4-5],網(wǎng)格大小的增大,導致平均地形指數(shù)相應增大,地下水埋深相應增加,模擬洪峰流量增加.另外一些研究也分析了不同網(wǎng)格尺度下的參數(shù)率定問題,如Molnar[6]通過比較分布式水文模型CASC2D 模擬結(jié)果對不同DEM 分辨率的敏感性,發(fā)現(xiàn)粗網(wǎng)格分辨率也可以進行水文模擬,但是要減小網(wǎng)格大小的影響,需要調(diào)整坡地匯流與河道糙率相關(guān)的參數(shù).
另一方面,水文模型的結(jié)構(gòu)具有尺度依賴性,隨意選擇網(wǎng)格大小影響著因模型結(jié)構(gòu)原因帶來的計算誤差.例如,MIKESHE 和LISEM 模型網(wǎng)格尺度效應研究都發(fā)現(xiàn),模型運動波求解時會產(chǎn)生數(shù)值誤差,這種誤差使洪水波發(fā)生擴散作用,且隨網(wǎng)格增大而加強,從而導致洪峰模擬隨著網(wǎng)格增大而減小[7-8].
網(wǎng)格尺度從這兩方面影響著模擬效果.因此,了解網(wǎng)格尺度變化對水文模擬的影響,對于提高模擬精度具有重要意義.本文通過多尺度試驗,分析DHSVM 模型在不同尺度下的不同徑流過程的響應特征,尋求合適的模擬尺度,為模型應用提供參考依據(jù).
選擇美國中部的伊利諾伊流域為研究流域,該流域位于阿肯色斯與俄克拉馬州的邊界(見圖1),曾被選為分布式水文模型比較計劃(DMIP)流域之一.流域氣候為大陸性濕潤氣候,處于冷暖氣流交替出現(xiàn)的地區(qū),年平均最大、最小氣溫分別為22℃、9℃.多年平均降水為1 200mm,多年平均水面蒸發(fā)為1 066 mm,40%的降水發(fā)生于3~6月,接下來的夏季為干季,9~11月為濕季,30%的降水發(fā)生在此階段,冬季流域為干旱季節(jié).流域出口為Tahlequah站,實際控制面積為2 483km2,多年平均流量為29m3/s,月平均最高流量發(fā)生在3~6月,最低流量發(fā)生于7~8月,基流占年平均流量的29.7%~72.5%(1936~2007).
流域內(nèi)有6個氣象站(如圖1所示),提供逐小時氣溫、長短波輻射、水汽壓、風速、降水等氣象數(shù)據(jù).植被分類數(shù)據(jù)來自國際地圈生物圈計劃(IGBP)的1km分辨率數(shù)據(jù)庫,土壤分類數(shù)據(jù)來自美國農(nóng)業(yè)部(USDA)全國土壤測量中心的STATSGO 數(shù)據(jù)庫,流域主要的植被類型是落葉闊葉林(91.8%),主要的土壤類型為粉質(zhì)壤土(46.6%),其次是粉質(zhì)黏壤土(38.1%).其他數(shù)據(jù)包括流域DEM 數(shù)據(jù)(30m×30m)、植被覆蓋度、水面蒸發(fā)、流量、河道斷面等數(shù)據(jù),大部分由NOAA 水文辦公室DMIP網(wǎng)站提供.
圖1 伊利諾伊流域位置圖
DHSVM(Distributed Hydrology Soils and Vegetation Model)模型由Mark Wigmosta于1994年提出[9],用于模擬中小尺度流域(通常小于10 000km2)內(nèi)土壤、植被和地形對地表以及地表以下水流運動的影響,目前最新的版本為DHSVM 3.0,為本研究采用版本.模型由7個模塊組成:蒸散發(fā)、地面降雪和融雪、冠層截雪和積雪融化、不飽和土壤水運動、飽和壤中流、飽和坡面流和河道流量演算.模型采用基于Penman-Monteith公式的雙層的樹冠計算截留和蒸散發(fā),采用質(zhì)量、能量平衡模型計算積融雪,水分的垂向運動采用一維的達西流公式,計算經(jīng)過多層土壤的不飽和壤中流.
DHSVM 采用近似運動波的方法逐網(wǎng)格計算壤中流,網(wǎng)格內(nèi)水流可以向周圍相鄰8個方向網(wǎng)格(0~7)流動.網(wǎng)格(i,j)在k 方向上的飽和壤中流的輸移率為:
式中,ωi,j,k是k 方向上網(wǎng)格流線密度;βi,j,k是k 方向上地下水位線的坡度;Ti,j(z,D)是網(wǎng)格輸移率.式(1)中的輸移率的計算公式如下:
式中,Ki,j是網(wǎng)格表層土壤側(cè)向飽和水力傳導率;fi,j是垂向衰減系數(shù);Di,j是網(wǎng)格土壤層厚度;zi,j是地下水位埋深;網(wǎng)格單元飽和壤中流的總出流等于式(1)計算的各方向的出流量之和.
模型坡面匯流也采用類似飽和壤中流計算的逐網(wǎng)格方法,并同時考慮了霍頓超滲產(chǎn)流、蓄滿產(chǎn)流與回歸流3種地面徑流.與壤中流計算不同的是坡面匯流中的流速采用的是定值,等于網(wǎng)格大小除以時間步長,意味著一個時間步長內(nèi)流出此網(wǎng)格的水量等于初始時刻的網(wǎng)格的儲水量.
模型河道的流量演算采用相對簡單但是很穩(wěn)健的線性槽蓄法,改進的馬斯京根法也可用于河道演算.線性槽蓄法對于大小不一、地形多變的流域,都能得到比較滿意的模擬結(jié)果,在DHSVM 的應用中,大部分的河道匯流演算均采用此方法.模型更多結(jié)構(gòu)介紹可參考文獻[9].
DHSVM 模型構(gòu)建需要處理大量數(shù)據(jù),除了上述提到的地形、植被、土壤數(shù)據(jù),還需要河網(wǎng)和土壤深度數(shù)據(jù),可通過自帶的處理程序生成.將柵格數(shù)據(jù)處理成統(tǒng)一的空間步長,時間序列數(shù)據(jù)處理成統(tǒng)一的時間步長,以此為驅(qū)動數(shù)據(jù),構(gòu)建模型.參考Wigmosta應用于一個相似大小流域的模擬效果[9],并經(jīng)過初步模擬結(jié)果比較,確定空間步長為200m、時間步長為3h以構(gòu)建模型.選取1994~1997年作為模型的率定期,1998~1999年為驗證期,對模型進行手動率定與驗證.
參數(shù)率定參照Wigmosta應用DHSVM 模型的經(jīng)驗,選擇最小植被氣孔阻抗、飽和側(cè)向傳導系數(shù)和側(cè)向傳導率垂向遞減指數(shù)3個最為敏感參數(shù)進行調(diào)整.其他植被土壤參數(shù)分別采用Land Data Assimilation System(LDA)提供的植被分類標準基本參數(shù)表與NOAA 水文辦公室提供的全球5min土壤參數(shù)分類基準參數(shù)表設定.率定過程,首先調(diào)整對水量平衡影響較大的最小植被氣孔阻抗,使總徑流偏差最小,然后再調(diào)整對洪水峰型影響較大的側(cè)向傳導系數(shù)與垂向遞減指數(shù),使模擬徑流過程盡量接近實測徑流過程.模擬效果采用確定性系數(shù)DC、徑流深相對誤差BIAS指標進行評價,如下式:
對模型率定期和驗證期的逐日平均徑流量的實測值和模擬值進行對比驗證,結(jié)果顯示率定期的確定性系數(shù)為0.71,總徑流偏差為-0.56%;驗證期效率系數(shù)為0.79,總徑流偏差為7.12%,如表1所示.從率定期和驗證期的洪水過程線(圖2a、2b)可以看出模型對洪峰的模擬不太滿意,低估了洪峰值.從總體上來看,模型在該流域具有適用性.率定后的模型參數(shù)用于尺度效應研究.
表1 模型率定期與驗證期模擬結(jié)果
圖2 模擬與實測流量比較
本文采用兩種方法來分析模型的網(wǎng)格尺度效應,一種是靜態(tài)參數(shù)法,即保持率定好的參數(shù)不變,利用3種不同網(wǎng)格大?。?00m×100m、200m×200m 和500m×500m)的DEM 提取流域特征,分別建立模型,分析空間模擬尺度變化對水文模擬的影響.另一種方法是變化參數(shù)法,選取模型最小植被氣孔阻抗、側(cè)向傳導系數(shù)和側(cè)向傳導率垂向遞減指數(shù)3個敏感參數(shù),每次只改變其中一個參數(shù),其他保持不變,分析不同參數(shù)對網(wǎng)格大小的敏感性.
將原始30m×30m 分辨率DEM 通過ARCGIS最鄰近插值方法分別重采樣到100m、200m、500m,并提取流域特征,統(tǒng)計結(jié)果見表2.由表中可見隨著網(wǎng)格大小的增加,DEM 高程各項(最大值、最小值、平均值與標準差)基本無變化,而坡度則隨著網(wǎng)格大小增加而減小,說明大網(wǎng)格對流域描述會產(chǎn)生地形的坦化現(xiàn)象,地形特征的空間異質(zhì)性變小.這與部分研究得到的結(jié)果相類似[2-4,10],即不同網(wǎng)格大小的DEM 得出的流域面積、高程大體上一致,但對與流域坡度有關(guān)的參數(shù)的影響較大.
表2 不同DEM 分辨率提取的地形特征
靜態(tài)參數(shù)情況下,選取一場洪水(1997 年12 月31日06時~1998年1月16日21時),比較各網(wǎng)格下的模擬徑流,發(fā)現(xiàn)隨著網(wǎng)格的增大,模擬洪峰流量相應增大(如圖3所示).其中500m 分辨率對降水最為敏感,不僅洪峰大,峰現(xiàn)時間也相應提前.而從基流部分的模擬結(jié)果來看,隨著網(wǎng)格的增大,基流流量呈減小趨勢.
圖3 不同網(wǎng)格尺度徑流模擬結(jié)果
為進一步分析引起不同網(wǎng)格大小下模擬徑流的差異的原因,將模型1997年12月31日~1998年1月16日降水、蒸發(fā)、不同產(chǎn)匯流成分統(tǒng)計見表3.
表3 不同網(wǎng)格尺度的產(chǎn)匯流模擬差異
從表3可見,隨著網(wǎng)格大小的變化,降水、蒸發(fā)項影響均較小,徑流深變化總體不大,但是徑流成分有較大變化.網(wǎng)格分辨率越細,壤中流比例越大,其中100m 網(wǎng)格模擬的壤中流比例最大,達到了0.72.由此可知,網(wǎng)格大小影響著匯流過程中不同徑流成分.不同網(wǎng)格尺度下匯入到河道的地表徑流、壤中流模擬結(jié)果如圖4所示.
圖4 匯入到河道的地表徑流、壤中流模擬結(jié)果比較
如2.1所述,DHSVM 模型計算坡面匯流采用定流速的計算方法,與平均坡度無關(guān),而與網(wǎng)格大小和河網(wǎng)密度有關(guān).為了消除河網(wǎng)密度對模擬結(jié)果的影響,生成河網(wǎng)時,各分辨率已采用相同的閾值,以保證河網(wǎng)密度相同.因此坡面匯流只與網(wǎng)格大小有關(guān),即網(wǎng)格越小,地表徑流要經(jīng)過更多的網(wǎng)格,即更長的時間才能被河網(wǎng)截留.在一場洪水中,大網(wǎng)格由于匯流路徑短,而降水往往是發(fā)生在有限的幾個時間步長,因此有多部分的水量通過地表徑流直接匯到河網(wǎng)[10],引起地表徑流增加.而對于壤中流采用的演算方法,每個時間步長流出網(wǎng)格各個方向的水量與局地坡度有關(guān),另外還與側(cè)向傳導系數(shù)、遞減系數(shù)有關(guān).當設定同樣的模型參數(shù),壤中流匯流只與坡度有關(guān).大網(wǎng)格引起地形坦化,平均坡度減小,壤中流出流緩慢,網(wǎng)格內(nèi)更容易滯蓄水量.與小網(wǎng)格相比,壤中流出流較小,模擬的平均土壤含水量較大.
圖5比較了1998年01月10日12:00(前期有降水,流域濕潤)各分辨率下模擬的地下水位埋深空間分布情況.由圖可見,隨著網(wǎng)格的增大,流域內(nèi)地下水位埋深普遍變淺,驗證了網(wǎng)格大小對流域土壤水分空間分布影響較大.理論上,對于一場洪水的模擬,因為網(wǎng)格增大,導致地表徑流增加,壤中流減小,這兩種機制呈中和作用,但是從圖4可以看出,洪水過程中,大網(wǎng)格的地表徑流程陡漲陡落趨勢,且占總徑流的比例更大,從而容易產(chǎn)生峰現(xiàn)時間提前,峰值提高的洪水.
圖5 不同網(wǎng)格模擬地下水埋深空間分布
比較本文與以往Topmodel模型網(wǎng)格尺度影響的研究結(jié)果發(fā)現(xiàn),網(wǎng)格大小對兩個模型模擬結(jié)果有相似的影響.對于Topmodel模型,隨著分辨率的降低,流域平均坡度值偏小,計算的地形指數(shù)偏大,導致模擬洪峰流量增大,平均地下水埋深變小,地表徑流占總徑流比例增大[11].模擬的土壤含水量空間分布也有類似的結(jié)果,大網(wǎng)格水量滯蓄導致更大的飽和面積[2].從兩者模型匯流模塊的結(jié)構(gòu)分析入手,可以發(fā)現(xiàn)其深層原因.DHSVM 模型的壤中流匯流方法又被稱為顯式匯流方法(Explicit Routing)[12],與Topmodel采用的基于統(tǒng)計理論的隱式方法(Implicit Routing)在機理上有一定程度的相似性,它們都假設土壤側(cè)向飽和水力傳導率隨土層深度呈指數(shù)遞減關(guān)系,且基于類似運動波方法計算壤中流輸移率.所不同的是Topmodel采用統(tǒng)計-動力方法計算流域土壤水的分布,簡化了計算量,使模型運算效率比DHSVM 高,但又不失其物理基礎.而DHSVM 模型采用逐網(wǎng)格方法將水流匯至流域出口,充分考慮了單元間的水流累積對徑流模擬的影響,物理概念更明確.Tague和Band[12]在一個小流域上比較過兩種匯流方法,發(fā)現(xiàn)兩種方法都能達到較好的徑流模擬效果,但隱式方法對土壤含水量的空間分布的模擬不如顯式方法精確,而且隱式方法模擬徑流量對地形數(shù)據(jù)的坦化和土壤傳導系數(shù)的變化更為敏感.因此,兩種相似的匯流方法在地形均化影響下,具有一致的網(wǎng)格尺度效應.值得注意的是,兩者的地表徑流計算方法并不相同,DHSVM 模型地表徑流影響直接來源于網(wǎng)格尺度的選擇,而Topmodel則來源于地形均化的影響.
為了分析網(wǎng)格大小對參數(shù)敏感性的影響,選擇最小植被氣孔阻抗、側(cè)向傳導系數(shù)和垂向遞減指數(shù)3個參數(shù),每次只改變一個參數(shù)值,其他參數(shù)保持前面率定好的不變,并假設改變的參數(shù)在空間上均勻分布,運行模型,統(tǒng)計1994~1999年各網(wǎng)格的模擬效果.
由圖6可見,各網(wǎng)格大小下的參數(shù)敏感性不一致.從確定性系數(shù)來看,最小植被阻抗與垂向遞減指數(shù)的敏感性在100m、200m 網(wǎng)格尺度下比較接近,500m 網(wǎng)格尺度下各參數(shù)敏感性比較大,且與其他網(wǎng)格尺度變化規(guī)律不一致.側(cè)向傳導系數(shù)的敏感性隨著網(wǎng)格大小的增大而增大,說明網(wǎng)格大小對側(cè)向傳導率的敏感性影響較大.由此可見,側(cè)向傳導系數(shù)是尺度效應較大參數(shù).而從水量誤差來看,各網(wǎng)格下的參數(shù)敏感性比較接近,網(wǎng)格越大,模擬水量負偏越大.
圖6 不同網(wǎng)格下的參數(shù)敏感性分析
Topmodel參數(shù)尺度效應分析研究發(fā)現(xiàn),產(chǎn)匯流過程中重要的敏感參數(shù)飽和導水率T0與網(wǎng)格分辨率有關(guān),網(wǎng)格增大后如要保證模型的精度保持穩(wěn)定有效,必須給T0賦較大的值[3].本研究中,500m 的網(wǎng)格尺度下要提高模擬精度,應增加側(cè)向傳導率的值.結(jié)合前面分析的匯流演算方法,當側(cè)向傳導系數(shù)增大時,壤中流更容易流出網(wǎng)格,水量滯蓄作用將減弱,一定程度減小了網(wǎng)格增大帶來的尺度影響.但此時模型具有物理意義的參數(shù)可能會變成有效參數(shù),在一定程度影響了模型的物理基礎,故不推薦使用500m 的網(wǎng)格進行建模.
小網(wǎng)格在大流域的應用上受運算效率的限制,本研究流域為中等大小流域,100m 步長將流域劃分成了606×840個網(wǎng)格,在一臺配置為I3處理器/4G 內(nèi)存的計算機上,模擬一年平均需要25min,而200m網(wǎng)格步長平均只需要5分鐘,小網(wǎng)格明顯增加了模型參數(shù)率定的難度.模型在200m 網(wǎng)格大小下的模擬效果表明此網(wǎng)格大小分辨率已經(jīng)滿足水文模擬的要求,同時又兼顧了模型運算效率.
總而言之,本文揭示的尺度效應源于兩方面,一是由于網(wǎng)格增大導致地形數(shù)據(jù)均化,均化后的坡度、河長等參數(shù),輸入到模型中,間接影響了模擬結(jié)果.另一方面的原因是模型結(jié)構(gòu)具有尺度依賴性,使網(wǎng)格大小的選擇對匯流過程產(chǎn)生直接影響.如前所述,分布式水文模型匯流,有不少采用運動波或擴散波方法的匯流模型,需要利用有限差分求解微分方程,網(wǎng)格離散大小的選擇則不可避免地會對求解過程帶來誤差.DHSVM 模型匯流算法在一定程度避免了差分求解時網(wǎng)格離散帶來的尺度效應,但地表徑流匯流模塊仍面臨著尺度問題.已有相關(guān)的模型,如改進的TOPKAIPI模型,通過引進控制性方程的空間積分和參數(shù)的平均化處理,旨在減輕網(wǎng)格尺度效應的影響,模型在一個流域上應用于從幾米到幾千米的網(wǎng)格尺度,物理意義和模型的計算精度并未受影響[1].DHSVM 開發(fā)初衷是應用于30~200m 的空間尺度和1~3h的時間尺度上,屬于比較精細的分布式水文模型.以往應用中,考慮到模型模擬效果與運行效率,小流域(幾十平方千米)一般采用30m 網(wǎng)格尺度,中等流域(幾百平方千米至數(shù)千平方千米)則可以采用100~200 m 的網(wǎng)格尺度,而大于200m 的網(wǎng)格尺度應用較少[9,13].在以后的應用研究中,應特別注意網(wǎng)格大小對模擬結(jié)果的影響,選擇合適的模擬尺度.
本文利用3種不同網(wǎng)格分辨率的DEM(100m、200m、500m)建模,分別通過靜態(tài)參數(shù)法(即保持200 m 下率定好的參數(shù)不變)與等步長變化參數(shù)的方法進行多尺度模擬,分析了各個模擬尺度下的水文要素與模型參數(shù)的尺度效應,主要結(jié)論如下:
1)在靜態(tài)參數(shù)條件下,發(fā)現(xiàn)網(wǎng)格大小對模擬蒸發(fā)量、總徑流量影響較小,而對洪峰和徑流過程的模擬影響較大.網(wǎng)格越大,洪峰峰值越大,地表徑流占總徑流比例越大,網(wǎng)格水量滯蓄作用增大.
2)用變化參數(shù)法分析網(wǎng)格大小對參數(shù)敏感性的影響,發(fā)現(xiàn)側(cè)向傳導系數(shù)的敏感性隨著網(wǎng)格增大而增大,明顯具有網(wǎng)格尺度效應.大網(wǎng)格要提高模擬效果,需要增加側(cè)向傳導率的值,但要特別注意參數(shù)的有效性.
3)水文模型只有應用于一定尺度范圍內(nèi)才會得到預期的模擬效果,本文對DHSVM 模型尺度適用性初步分析可見,模型網(wǎng)格劃分時應充分考慮網(wǎng)格選擇對不同水文過程內(nèi)部機理的影響,本研究表明200 m 的分辨率已能達到較滿意的模擬效果.
[1] 徐宗學.水文模型[M].北京:科學出版社,2009.
[2] Kuo W L,Steenhuis T S,McCulloch C E,et al.Effect of Grid Size on Runoff and Soil Moisture for a Variable-Source-Area Hydrology Model[J].Water resources research,1999,35(11):3419-3428.
[3] 孫立群,胡成,陳 剛.TOPMODEL 模型中的DEM尺度效應[J].水科學進展,2008,19(5):699-706.
[4] Vieux B E.DEM Aggregation and Smoothing Effects on Surface Runoff Modeling[J].Journal of Computing in Civil Engineering,1993,7(3):310-338.
[5] Zhang W,Montgomery D R.Digital Elevation Model Grid Size,landscape Representation[J].Water resources research,1994,30(4):1019-1028.
[6] Molnar D,Julien P.Grid-size Effects on Surface Runoff Modeling[J].Journal of Hydrologic Engineering,2000,5(1):8-16.
[7] Ali M,R S,Ali R,et al.Simulations of Varying Grid Sizes on Catchment Yield by Using Calibrated and Validated MIKE SHE Models[M].18th World IMACS/MODSIM Congress.Cairns,Australia.2009.
[8] Hessel R.Effects of Grid Cell Size and Time Step Length on Simulation Results of the Limburg Soil Erosion Model (LISEM)[J].Hydrological Processes,2005,19(15):3037-3049.
[9] 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-1680.
[10]Dubin A M.Assessing the Influence of Digital Elevation Model Resolution in Hydrologic Modeling[M].University of Washington,1998.[11]Wolock D M,Price C V.Effects of Digital Elevation Model Map Scale and Data Resolution on a Topography-Based Watershed Model[J].Water resources research,1994,30(11):3041-3052.
[12]Tague C,Band L.Evaluating Explicit and Implicit Routing for Watershed Hydro-Ecological Models of Forest Hydrology at the Small Catchment Scale[J].Hydrological Processes,2001,15(8):1415-1439.
[13]Vanshaar J R,Haddeland I,Lettenmaier D P.Effects of Land-cover Changes on the Hydrological Response of Interior Columbia River Basin Forested Catchments[J].Hydrological Processes,2002:2499-2520.