分布式水沙物理模型CASC2D-SED的應(yīng)用研究

李致家，沈 潔，張鵬程，李 娟，姚 成，郭 元

(1.河海大學水文水資源學院，江蘇 南京 210098;2.寧夏瑞沃水資源工程研究院，寧夏銀川 750021;3.寧夏大學土木與水利工程學院，寧夏銀川 750021)

流域的水沙物理過程非常復雜。對流域侵蝕產(chǎn)沙與輸沙過程的模擬，國外進行了大量研究，開發(fā)了通用土壤流失方程(USLE)［1］、產(chǎn)流產(chǎn)沙預報模型(WEPP)［2］、荷蘭模型(LISEM)［3］、歐洲產(chǎn)流產(chǎn)沙預報模型(EROSEM)［4］等。最初的模型大多屬于集總式模型，使用均一化處理的參數(shù)對水文泥沙過程進行平均化模擬，其結(jié)果無法真實地體現(xiàn)流域水文泥沙過程的時空分布特性［5］。隨著計算機技術(shù)、GIS以及RS技術(shù)的發(fā)展，分布式模型成為當今水文預報、土壤侵蝕領(lǐng)域的研究重點與發(fā)展方向。

本文以分布式水沙物理模型CASC2D-SED的產(chǎn)流模塊在國內(nèi)濕潤流域——舒家流域適用性研究為基礎(chǔ)［6］，選取寧夏三關(guān)口流域及好水川流域為典型研究流域，進一步研究其產(chǎn)流及產(chǎn)沙模塊在國內(nèi)流域的適用性，以及產(chǎn)流模塊在無資料地區(qū)徑流模擬中的適用性。

1 CASC2D-SED模型簡述

CASC2D模型最初的結(jié)構(gòu)源于Julien教授對二維地面徑流計算方法的發(fā)展［6-7］，他采用APL語言編寫計算程序，之后模型被逐漸加入Green-Ampt下滲計算［8-9］和一維顯式擴散波河道計算［10］、二維土壤侵蝕算法［11］、地下徑流模擬等，使得模型更加完善、系統(tǒng)化，模型名稱也變更為CASC2D-SED。

1.1 水流演算部分

1.1.1 降雨計算

當流域上只有1個雨量站時，認為降雨強度均勻分布，即每個柵格單元上的降雨強度與此雨量站的降雨強度一致;當流域上有多個雨量站時，采用距離平方倒數(shù)法插值估計每個柵格單元上的降雨強度，降雨強度在流域上是分布式的:

式中:rt(j，k)——在 t時刻柵格(j，k)處的降雨強度，mm/h;rm，t(jm，km)——位于柵格(jm，km)處的 m 雨量站在t時刻的實測降雨強度，mm/h;dm——柵格(j，k)與m雨量站所在柵格(jm，km)之間的距離，m;M——雨量站總數(shù)。

1.1.2 截留計算

在CASC2D-SED中，計算下滲之前，先從降雨中扣除截留深，只有當累積的降雨深等于植物截留深I(lǐng)時，其后的降雨才不再扣除截留損失。I根據(jù)柵格對應(yīng)的土地利用類型求得。

1.1.3 產(chǎn)流計算

在CASC2D-SED模型中，采用Green-Ampt產(chǎn)流模式。用Green-Ampt方程估算流域上每個柵格單元相應(yīng)的下滲率:

式中:f——下滲率，cm/h;Ks——飽和土壤水力傳導度，cm/h;Hc——毛管水頭，cm;Md——土壤缺水量，cm3/cm3;F——累積下滲深度，cm。

1.1.4 坡面匯流計算

在CASC2D-SED中，坡面徑流采用擴散波的二維顯式有限差分來計算。

二維連續(xù)方程:

式中:ho——地面徑流的深度;qx，qy——x，y 方向上的單寬流量;e——超滲降雨量(p － f);p——降雨強度。

二維擴散波水流方程(簡化動量方程):

式中:Sox，Soy——x和y方向的坡底比降，分別由數(shù)字高程模型計算;Sfx，Sfy——x和y方向的坡底摩阻比降。

曼寧阻力方程:

式中n是曼寧糙率系數(shù)。

1.1.5 河道匯流計算

CASC2D-SED模型認為河道水流為一維明渠流，河道匯流計算采用一維顯式有限差分擴散波方法。

一維連續(xù)方程:

式中:A——水流斷面面積;Q——河道流量;q0——單位長度河道上的旁側(cè)入流或出流。

一維擴散波方程:

式中:Sf—— 河道摩擦比降;Sc—— 河床底坡;hc—— 河道水深。

1.2 泥沙演算部分

1.2.1 地表侵蝕和泥沙輸運

對于一個大小為Δx的網(wǎng)格，在一個時間間隔Δt內(nèi)，網(wǎng)格內(nèi)能產(chǎn)生的總泥沙體積?SKR由改進的Kilinc＆Richardson輸運能力方程(KR方程)計算，公式如下:

式中:S0——地表坡度;q——單寬流量;K——土壤侵蝕因數(shù);C——覆蓋管理因數(shù);P——實踐因數(shù)。

對流過程輸運的粒級i中懸浮質(zhì)運移量為

式中:V——平均流速，m/s;SSUSi——粒級i的懸浮顆粒量。

從源網(wǎng)格輸移到接收網(wǎng)格的?SSUSi是對流輸運量和KR方程計算的最大值:

從源網(wǎng)格輸移到接收網(wǎng)格的粒級i床沙質(zhì)運移量?SBMi為

式中:SBMi——源網(wǎng)格粒級i床沙質(zhì)體積;CET——水流的剩余運移能力，由下式計算:

如果懸浮質(zhì)和床沙質(zhì)都被輸運的情況下還有運移能力剩余，土壤母質(zhì)會根據(jù)粒級i相應(yīng)比例地被侵蝕:

式中:?SEROSi——粒級i母質(zhì)侵蝕容量;Pi——粒級i在母質(zhì)土壤中的侵蝕比例。

該算法提供了較好的解決流域產(chǎn)沙和水流運移能力有限條件下的泥沙輸運方法。

1.2.2 河道泥沙輸運

流域坡面侵蝕泥沙是通過河道輸移到出口斷面的。河道目前不允許侵蝕，但允許泥沙沉降。采用Engelund＆Hansen方程計算河道泥沙運移能力。

單位時間Δt內(nèi)粒級i泥沙顆粒的輸移體積被估算為

式中CWi為按質(zhì)量計算的粒級i泥沙平均質(zhì)量濃度。

河道中通過對流作用輸運的粒級i懸浮質(zhì)運移量為

用來搬運河道中粒級i床沙質(zhì)的剩余運移能力為

河道輸運的粒級i床沙質(zhì)運移量是剩余運移能力和粒級i對流運移床沙質(zhì)的最小值:

如果還有輸移能力剩余，這部分剩余輸移能力將不使用。

1.2.3 懸浮泥沙沉降

假定粗沙、粉沙和黏土顆粒的中值粒徑d已知，表1為其沉降速度的估算。

懸浮顆粒部分在坡地和河道網(wǎng)格中都允許沉降，CASC2D-SED假定一個離散的顆粒沉降模式，顆粒彼此獨立沉降，互不影響，粒級i懸浮泥沙沉降比率PSi由下式確定:

表1 顆粒中值直徑和沉降速度Table 1 Particle’s median diameter and fall velocity

式中:wi——中值粒徑估算的沉降速度，m/s;h——網(wǎng)格水深，m。

模型計算過程中一個時間間隔內(nèi)網(wǎng)格懸浮泥沙沉降量從懸浮泥沙中減去，加到沉積沉降量中。該算法提高了回水區(qū)域懸移質(zhì)顆粒連續(xù)沉降算法的計算效率。

1.3 單元流域劃分與參數(shù)

采用正交網(wǎng)格劃分單元流域。模型參數(shù)有:植物截留深;土壤飽和水力傳導度，毛管水頭，土壤缺水量;坡面的曼寧糙率系數(shù);河道的寬度、深度、糙率等。這些參數(shù)都是柵格式空間分布的，其中植物截留深與柵格的土地利用相關(guān)，下滲參數(shù)與土壤類型相關(guān)。河道寬度和深度等參數(shù)的率定應(yīng)以實際資料為參考。本文采用人工試錯法進行參數(shù)率定。

2 典型研究流域概況

2.1 流域Ⅰ概況

選取寧夏三關(guān)口流域為典型研究流域Ⅰ，流域有產(chǎn)期的降雨徑流資料。流域位于六盤山地區(qū)的固原市涇源縣境內(nèi)，即東經(jīng)106°～107°、北緯35°～36°之間，是涇河上游支流頡河的上游，總面積為218 km2。氣候?qū)贉貛駶?、半濕潤氣候區(qū)，具有春寒、無夏、秋短、冬長的特點，年平均氣溫5.7℃。

流域多年平均降水量547 mm，降水量年內(nèi)分布有明顯的典型大陸性特點，主要集中在6—9月。流域多年平均徑流深111 mm，多年平均徑流量2777萬m3。

2.2 流域Ⅱ概況

選取寧夏好水川流域為典型研究流域Ⅱ。好水川流域缺乏徑流觀測資料，但有長期的降雨資料。流域位于六盤山西側(cè)，寧夏回族自治區(qū)隆德縣境內(nèi)，是渭河上游支流葫蘆河上游左岸一級支流的上游，地理坐標為東經(jīng)105°56'～106°15'，北緯35°38'～35°45'，平均氣溫5.1℃，流域總面積為102km2。流域多年平均降水量517mm，多年平均水面蒸發(fā)量900 mm，干旱指數(shù)為1.5，流域內(nèi)由于土壤沙化，植被覆蓋率低，流水侵蝕比較嚴重。

3 模型在典型研究流域Ⅰ的應(yīng)用

圖1 三關(guān)口流域河道及站點分布Fig.1 Channels and stations of Sanguankou Watershed

3.1 模型輸入及參數(shù)率定

模型的輸入數(shù)據(jù)包括降雨量、實測流量、實測沙量、流域90 m×90 m DEM、河道數(shù)據(jù)(圖1)、土地利用(圖2)以及土壤類型資料。DEM資料來自于美國地質(zhì)調(diào)查局(USGS)公共域免費提供的90 m×90 m的DEM數(shù)據(jù)，利用arcGIS軟件處理DEM數(shù)據(jù)。流域內(nèi)有4個雨量站:大灣、什字、清水溝、三關(guān)口，其中清水溝及三關(guān)口也是流量站，本文采用三關(guān)口的流量資料，其站點分布見圖1。

模擬時段步長選取2 s，降雨資料的輸入時段長是0.5 h，在一個時段內(nèi)雨強是不變的。實測流量、沙量數(shù)據(jù)的時間間隔不均勻，采用線性插值，將實測流量、沙量資料的時段長整理為0.5 h。模型的輸出時段長為計算步長的倍數(shù)，輸出時段可以是1 min，5 min。為了對應(yīng)于整理的實測流量、沙量時段長，將模型的模擬流量、沙量輸出時段長也定為0.5 h。

首先根據(jù)實測數(shù)據(jù)確定參數(shù)，然后根據(jù)流域DEM數(shù)據(jù)、土地利用類型、土壤類型及河道特征估算模型參數(shù)，最后對一些敏感的參數(shù)如土壤飽和水力傳導度、毛管水頭、土壤侵蝕因數(shù)等再進行率定。產(chǎn)流參數(shù)率定、驗證完成后，固定產(chǎn)流參數(shù)再進行泥沙參數(shù)的率定及驗證。模型參數(shù)見表2和表3。

圖2 三關(guān)口流域土地利用情況Fig.2 Land use in Sanguankou Watershed

表2 土地利用參數(shù)Table 2 Land use parameters

表3 土壤參數(shù)值及組成Table 3 Soil parameters and percentage composition

3.2 模型模擬結(jié)果

對三關(guān)口流域1983—1987年的8場洪水進行產(chǎn)流、產(chǎn)沙模擬。其中前5場洪水對模型參數(shù)進行率定，后3場洪水對模型進行驗證。首先進行產(chǎn)流參數(shù)的率定及驗證，產(chǎn)流模擬結(jié)果的相關(guān)特征值見表4。

表4 CASC2D-SED模型產(chǎn)流模擬結(jié)果特征值Table 4 Runoff statistical values simulated by CASC2D-SED model in Sanguankou Watershed

固定產(chǎn)流參數(shù)，進行泥沙參數(shù)的率定及驗證。泥沙模擬結(jié)果的相關(guān)特征值見表5。

19860624場洪水的流量、沙量實測模擬對比見圖3。從模擬結(jié)果看，模型在三關(guān)口流域的產(chǎn)流模擬效果較好，8場洪水的確定性系數(shù)均在0.70以上，且均值為0.78。洪峰相對誤差合格率為100%，洪量相對誤差合格率為88%，只有19860624場洪水洪量相對誤差不合格。產(chǎn)沙模擬效果良好，8場洪水中有6場洪水的確定性系數(shù)在0.70以上，且均值為0.76，沙峰相對誤差合格率為75%，洪水呈陡漲陡落情況。

表5 CASC2-SED模型產(chǎn)沙模擬結(jié)果特征值Table 5 Sediment statistical values simulated by CASC2D-SED model

圖3 19860624場洪水實測模擬對比過程Fig.3 Hydrographs of observed and simulated data of flood No.19860624

a.產(chǎn)沙模擬精度略低于產(chǎn)流模擬精度。由于未能融入流域下墊面［12－13］變化信息，得到的土地利用類型及土壤類型分布不精確，這必然會影響模型的模擬精度，同時會影響模型參數(shù)率定，而產(chǎn)沙模塊中的參數(shù)較產(chǎn)流模塊中多，故產(chǎn)沙模擬精度低于產(chǎn)流模擬精度。

b.洪峰及沙峰普遍偏低。這是由于每個柵格上的降水量均由研究區(qū)域內(nèi)雨量站的實測降水資料整理得到，降雨的時空分布與實際情形有一定的出入，給產(chǎn)流、產(chǎn)沙計算帶來一定的誤差。另外，由于實測資料時段長取0.5 h，原始記錄降雨量、流量和沙量資料有一定程度的均化，可能會錯過降雨、流量和沙量的峰值。

c.從模擬的洪水過程線來看，洪水退水過程快，這是由于CASC2D-SED模型在產(chǎn)流計算上采用的是基于霍頓產(chǎn)流機制的Green-Ampt下滲，沒有壤中流和地下徑流的支持，洪水過程呈現(xiàn)陡漲陡落情況。

4 模型在典型研究流域Ⅱ的應(yīng)用

4.1 研究方法

對于有實測徑流資料的流域，一般是運用實際的觀測資料對模型參數(shù)進行率定、驗證之后再用于徑流模擬;對于無實測徑流資料的流域，一般采用區(qū)域化的思想，包括:(a)以包括所研究的無資料子流域在內(nèi)的一個大的區(qū)域作為研究對象，應(yīng)用該區(qū)域內(nèi)其他有實測資料的流域優(yōu)選模型參數(shù)，并對優(yōu)選的參數(shù)值求均值，之后再移用于該無資料流域;(b)通過地形特征推求模型參數(shù);(c)應(yīng)用區(qū)域回歸分析推求模型參數(shù);(d)直接移用臨近的相似流域資料。本文選用第(d)方法，將黃河上游三關(guān)口流域的參數(shù)移用于臨近無資料流域——好水川流域，檢驗CASC2D-SED模型的產(chǎn)流模塊在無資料地區(qū)的適用性。

4.2 模型輸入及參數(shù)率定

模型的輸入數(shù)據(jù)包括降雨、實測流量、流域90 m×90m DEM、河道數(shù)據(jù)(圖4)、土地利用(圖5)以及土壤類型資料。土地利用類型有:林地、水面、耕地、草地。流域內(nèi)土壤為黑壚土。流域內(nèi)有2個雨量站:張銀、楊河，15個計算點:陰山、蘭家灣、范灣、上岔、張家臺子、后溝、下老莊、團結(jié)、后海子、張銀、岔口、下岔、何家岔、陽山莊、老張溝。

模型時段步長選取2 s，降雨資料的輸入時段長0.5 h，實測流量資料的時段長及模擬流量輸出時段長均定為0.5 h。

由于好水川流域與三關(guān)口流域多年平均降水量均超過500 mm，降雨年內(nèi)分布均有明顯的典型大陸性特點，主要集中在6—9月，且兩流域均位于六盤山地區(qū)，植被相近，流域內(nèi)土壤均為黑壚土，所以好水川流域的參數(shù)移用三關(guān)口流域的參數(shù)。

4.3 模型模擬結(jié)果

好水川流域缺乏徑流觀測值，但有長期的降雨資料。根據(jù)楊河站及張銀站2010—2011年4—9月的實測降雨資料，采用CASC2D-SED中的產(chǎn)流模塊進行降雨徑流模擬，得到15個計算點的流量過程，并計算得到日來水量。由于好水川流域無實測徑流資料進行對比，故將模擬計算得到的日來水量交由相關(guān)部門驗證。經(jīng)驗證，模擬計算得到的日來水量是合理的，說明將三關(guān)口流域的參數(shù)移用于好水川流域后，CASC2D-SED模型的計算結(jié)果是可靠的。

圖4 好水川流域河道Fig.4 Channels of Haoshuichuan Watershed

從以上結(jié)果看，CASC2D-SED模型在三關(guān)口流域的產(chǎn)流模擬中計算精度良好，將其移用于相似流域好水川流域時，模擬結(jié)果也是合理的，說明CASC2D-SED模型可用于無資料地區(qū)的徑流模擬。

圖5 好水川流域的土地利用Fig.5 Land use in Haoshuichuan Watershed

5 結(jié) 語

隨著地理信息系統(tǒng)、遙感技術(shù)以及氣象衛(wèi)星的進一步發(fā)展，資料的空間分辨率將會更高，可以更準確地反映流域的地形、土地利用、土壤和河道資料，分布式水沙物理模型將會得到更好的數(shù)據(jù)支持;在解決資料問題后，隨著計算機技術(shù)的進一步發(fā)展，模型自身結(jié)構(gòu)的不斷優(yōu)化以及模型算法的不斷完善，分布式水沙物理模型CASC2D-SED的運行速度將不斷提高，模擬精度將會更加優(yōu)良且穩(wěn)定。

［1］WISCHMEIER W H，SMITH D D.Predicting rainfall erosion losses:a guide to conservation planning［M］.Washington D.C.:Agriculture Handbook，1978.

［2］NEARING M A，F(xiàn)OSTER G R，LANE L J，et al.A process-based soil erosion model for USDA-water erosion prediction project technology［J］.Transactions of the ASAE，1989，32(5):1587-1593.

［3］de ROO A PJ，WESSELING CG，RITSEMA CJ.LISEM:a single－event physically based hydrological and soil erosion model for drainage basins.I:theory，input and output［J］.Hydrological Processes，1996，10(8):1107-1118.

［4］MORGAN R P C，QUINTON J N，SMITH R E，et al.The European soil erosion model(EUROSEM):a dynamic approach for predicting sediment transport from fields and small catchments［J］.Earth Surface Processes and Landforms，1998，23:527-544.

［5］曹文洪，祁偉，郭慶超，等.小流域產(chǎn)匯流分布式模型［J］.水利學報，2003，34(9):48-54.(CAO Wenhong，QI Wei，GUO Qingchao，et al.Distributed model for simulating runoff yield in small watershed［J］.Journal of Hydraulic Engineering，2003，34(9):48-54.(in Chinese))

［6］李致家，胡偉升，丁杰，等.基于物理基礎(chǔ)與基于柵格的分布式水文模型研究［J］.水力發(fā)電學報，2012(2):5-13.(LI Zhijia，HU Weisheng，DING Jie，et al.Study on distributed hydrological model of solving physical equation on grids［J］.Journal of Hydroelectric Engineering，2012(2):5-13.(in Chinese))

［7］JULIEN P Y，SAGHAFIAN B.CASC2D user manual-a two dimensional watershed rainfall-runoff Model［M］.Fort Collins:Colorado State University，1991.

［8］GREEN W H，AMPT G A.The flow of air and water through soils［J］.The Journal of Agricultural Science，1911，4:1-24.

［9］RAY K，LINSLEY J，MAX A，et al.Hydrology for engineers［M］.3rd ed.New York:McGraw-Hill，Inc，1982.

［10］OGDEN F L.De-St Venant channel routing in distributed hydrologic modeling［J］.Hydraulic Engineering，1994，1:492-496.

［11］JULIEN P，ROJASR.Upland erosion modeling with CASC2D-SED［J］.International Journal of Sediment Research，2002，17(4):265-274.

［12］MAIDMENT D R.Handbook of hydrology［M］.New York:McGraw-Hill，Inc，1993.

［13］ANDERSON M G，MCDONNELL J J.Encyclopedia of hydrologic science［M］.New York:John Wiley ＆ Sons Ltd，2005.