分布式水沙物理模型CASC2D-SED的应用研究

李致家，沈 洁，张鹏程，李 娟，姚 成，郭 元

(1.河海大学水文水资源学院，江苏 南京 210098;2.宁夏瑞沃水资源工程研究院，宁夏银川 750021;3.宁夏大学土木与水利工程学院，宁夏银川 750021)

流域的水沙物理过程非常复杂。对流域侵蚀产沙与输沙过程的模拟，国外进行了大量研究，开发了通用土壤流失方程(USLE)［1］、产流产沙预报模型(WEPP)［2］、荷兰模型(LISEM)［3］、欧洲产流产沙预报模型(EROSEM)［4］等。最初的模型大多属于集总式模型，使用均一化处理的参数对水文泥沙过程进行平均化模拟，其结果无法真实地体现流域水文泥沙过程的时空分布特性［5］。随着计算机技术、GIS以及RS技术的发展，分布式模型成为当今水文预报、土壤侵蚀领域的研究重点与发展方向。

本文以分布式水沙物理模型CASC2D-SED的产流模块在国内湿润流域——舒家流域适用性研究为基础［6］，选取宁夏三关口流域及好水川流域为典型研究流域，进一步研究其产流及产沙模块在国内流域的适用性，以及产流模块在无资料地区径流模拟中的适用性。

1 CASC2D-SED模型简述

CASC2D模型最初的结构源于Julien教授对二维地面径流计算方法的发展［6-7］，他采用APL语言编写计算程序，之后模型被逐渐加入Green-Ampt下渗计算［8-9］和一维显式扩散波河道计算［10］、二维土壤侵蚀算法［11］、地下径流模拟等，使得模型更加完善、系统化，模型名称也变更为CASC2D-SED。

1.1 水流演算部分

1.1.1 降雨计算

当流域上只有1个雨量站时，认为降雨强度均匀分布，即每个栅格单元上的降雨强度与此雨量站的降雨强度一致;当流域上有多个雨量站时，采用距离平方倒数法插值估计每个栅格单元上的降雨强度，降雨强度在流域上是分布式的:

式中: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——雨量站总数。

1.1.2 截留计算

在CASC2D-SED中，计算下渗之前，先从降雨中扣除截留深，只有当累积的降雨深等于植物截留深I时，其后的降雨才不再扣除截留损失。I根据栅格对应的土地利用类型求得。

1.1.3 产流计算

在CASC2D-SED模型中，采用Green-Ampt产流模式。用Green-Ampt方程估算流域上每个栅格单元相应的下渗率:

式中:f——下渗率，cm/h;Ks——饱和土壤水力传导度，cm/h;Hc——毛管水头，cm;Md——土壤缺水量，cm3/cm3;F——累积下渗深度，cm。

1.1.4 坡面汇流计算

在CASC2D-SED中，坡面径流采用扩散波的二维显式有限差分来计算。

二维连续方程:

式中:ho——地面径流的深度;qx，qy——x，y 方向上的单宽流量;e——超渗降雨量(p － f);p——降雨强度。

二维扩散波水流方程(简化动量方程):

式中:Sox，Soy——x和y方向的坡底比降，分别由数字高程模型计算;Sfx，Sfy——x和y方向的坡底摩阻比降。

曼宁阻力方程:

式中n是曼宁糙率系数。

1.1.5 河道汇流计算

CASC2D-SED模型认为河道水流为一维明渠流，河道汇流计算采用一维显式有限差分扩散波方法。

一维连续方程:

式中:A——水流断面面积;Q——河道流量;q0——单位长度河道上的旁侧入流或出流。

一维扩散波方程:

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

1.2 泥沙演算部分

1.2.1 地表侵蚀和泥沙输运

对于一个大小为Δx的网格，在一个时间间隔Δt内，网格内能产生的总泥沙体积∀SKR由改进的Kilinc＆Richardson输运能力方程(KR方程)计算，公式如下:

式中:S0——地表坡度;q——单宽流量;K——土壤侵蚀因数;C——覆盖管理因数;P——实践因数。

对流过程输运的粒级i中悬浮质运移量为

式中:V——平均流速，m/s;SSUSi——粒级i的悬浮颗粒量。

从源网格输移到接收网格的∀SSUSi是对流输运量和KR方程计算的最大值:

从源网格输移到接收网格的粒级i床沙质运移量∀SBMi为

式中:SBMi——源网格粒级i床沙质体积;CET——水流的剩余运移能力，由下式计算:

如果悬浮质和床沙质都被输运的情况下还有运移能力剩余，土壤母质会根据粒级i相应比例地被侵蚀:

式中:∀SEROSi——粒级i母质侵蚀容量;Pi——粒级i在母质土壤中的侵蚀比例。

该算法提供了较好的解决流域产沙和水流运移能力有限条件下的泥沙输运方法。

1.2.2 河道泥沙输运

流域坡面侵蚀泥沙是通过河道输移到出口断面的。河道目前不允许侵蚀，但允许泥沙沉降。采用Engelund＆Hansen方程计算河道泥沙运移能力。

单位时间Δt内粒级i泥沙颗粒的输移体积被估算为

式中CWi为按质量计算的粒级i泥沙平均质量浓度。

河道中通过对流作用输运的粒级i悬浮质运移量为

用来搬运河道中粒级i床沙质的剩余运移能力为

河道输运的粒级i床沙质运移量是剩余运移能力和粒级i对流运移床沙质的最小值:

如果还有输移能力剩余，这部分剩余输移能力将不使用。

1.2.3 悬浮泥沙沉降

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

悬浮颗粒部分在坡地和河道网格中都允许沉降，CASC2D-SED假定一个离散的颗粒沉降模式，颗粒彼此独立沉降，互不影响，粒级i悬浮泥沙沉降比率PSi由下式确定:

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

式中:wi——中值粒径估算的沉降速度，m/s;h——网格水深，m。

模型计算过程中一个时间间隔内网格悬浮泥沙沉降量从悬浮泥沙中减去，加到沉积沉降量中。该算法提高了回水区域悬移质颗粒连续沉降算法的计算效率。

1.3 单元流域划分与参数

采用正交网格划分单元流域。模型参数有:植物截留深;土壤饱和水力传导度，毛管水头，土壤缺水量;坡面的曼宁糙率系数;河道的宽度、深度、糙率等。这些参数都是栅格式空间分布的，其中植物截留深与栅格的土地利用相关，下渗参数与土壤类型相关。河道宽度和深度等参数的率定应以实际资料为参考。本文采用人工试错法进行参数率定。

2 典型研究流域概况

2.1 流域Ⅰ概况

选取宁夏三关口流域为典型研究流域Ⅰ，流域有产期的降雨径流资料。流域位于六盘山地区的固原市泾源县境内，即东经106°～107°、北纬35°～36°之间，是泾河上游支流颉河的上游，总面积为218 km2。气候属温带湿润、半湿润气候区，具有春寒、无夏、秋短、冬长的特点，年平均气温5.7℃。

流域多年平均降水量547 mm，降水量年内分布有明显的典型大陆性特点，主要集中在6—9月。流域多年平均径流深111 mm，多年平均径流量2777万m3。

2.2 流域Ⅱ概况

选取宁夏好水川流域为典型研究流域Ⅱ。好水川流域缺乏径流观测资料，但有长期的降雨资料。流域位于六盘山西侧，宁夏回族自治区隆德县境内，是渭河上游支流葫芦河上游左岸一级支流的上游，地理坐标为东经105°56'～106°15'，北纬35°38'～35°45'，平均气温5.1℃，流域总面积为102km2。流域多年平均降水量517mm，多年平均水面蒸发量900 mm，干旱指数为1.5，流域内由于土壤沙化，植被覆盖率低，流水侵蚀比较严重。

3 模型在典型研究流域Ⅰ的应用

图1 三关口流域河道及站点分布Fig.1 Channels and stations of Sanguankou Watershed

3.1 模型输入及参数率定

模型的输入数据包括降雨量、实测流量、实测沙量、流域90 m×90 m DEM、河道数据(图1)、土地利用(图2)以及土壤类型资料。DEM资料来自于美国地质调查局(USGS)公共域免费提供的90 m×90 m的DEM数据，利用arcGIS软件处理DEM数据。流域内有4个雨量站:大湾、什字、清水沟、三关口，其中清水沟及三关口也是流量站，本文采用三关口的流量资料，其站点分布见图1。

模拟时段步长选取2 s，降雨资料的输入时段长是0.5 h，在一个时段内雨强是不变的。实测流量、沙量数据的时间间隔不均匀，采用线性插值，将实测流量、沙量资料的时段长整理为0.5 h。模型的输出时段长为计算步长的倍数，输出时段可以是1 min，5 min。为了对应于整理的实测流量、沙量时段长，将模型的模拟流量、沙量输出时段长也定为0.5 h。

首先根据实测数据确定参数，然后根据流域DEM数据、土地利用类型、土壤类型及河道特征估算模型参数，最后对一些敏感的参数如土壤饱和水力传导度、毛管水头、土壤侵蚀因数等再进行率定。产流参数率定、验证完成后，固定产流参数再进行泥沙参数的率定及验证。模型参数见表2和表3。

图2 三关口流域土地利用情况Fig.2 Land use in Sanguankou Watershed

表2 土地利用参数Table 2 Land use parameters

表3 土壤参数值及组成Table 3 Soil parameters and percentage composition

3.2 模型模拟结果

对三关口流域1983—1987年的8场洪水进行产流、产沙模拟。其中前5场洪水对模型参数进行率定，后3场洪水对模型进行验证。首先进行产流参数的率定及验证，产流模拟结果的相关特征值见表4。

表4 CASC2D-SED模型产流模拟结果特征值Table 4 Runoff statistical values simulated by CASC2D-SED model in Sanguankou Watershed

固定产流参数，进行泥沙参数的率定及验证。泥沙模拟结果的相关特征值见表5。

19860624场洪水的流量、沙量实测模拟对比见图3。从模拟结果看，模型在三关口流域的产流模拟效果较好，8场洪水的确定性系数均在0.70以上，且均值为0.78。洪峰相对误差合格率为100%，洪量相对误差合格率为88%，只有19860624场洪水洪量相对误差不合格。产沙模拟效果良好，8场洪水中有6场洪水的确定性系数在0.70以上，且均值为0.76，沙峰相对误差合格率为75%，洪水呈陡涨陡落情况。

表5 CASC2-SED模型产沙模拟结果特征值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.产沙模拟精度略低于产流模拟精度。由于未能融入流域下垫面［12－13］变化信息，得到的土地利用类型及土壤类型分布不精确，这必然会影响模型的模拟精度，同时会影响模型参数率定，而产沙模块中的参数较产流模块中多，故产沙模拟精度低于产流模拟精度。

b.洪峰及沙峰普遍偏低。这是由于每个栅格上的降水量均由研究区域内雨量站的实测降水资料整理得到，降雨的时空分布与实际情形有一定的出入，给产流、产沙计算带来一定的误差。另外，由于实测资料时段长取0.5 h，原始记录降雨量、流量和沙量资料有一定程度的均化，可能会错过降雨、流量和沙量的峰值。

c.从模拟的洪水过程线来看，洪水退水过程快，这是由于CASC2D-SED模型在产流计算上采用的是基于霍顿产流机制的Green-Ampt下渗，没有壤中流和地下径流的支持，洪水过程呈现陡涨陡落情况。

4 模型在典型研究流域Ⅱ的应用

4.1 研究方法

对于有实测径流资料的流域，一般是运用实际的观测资料对模型参数进行率定、验证之后再用于径流模拟;对于无实测径流资料的流域，一般采用区域化的思想，包括:(a)以包括所研究的无资料子流域在内的一个大的区域作为研究对象，应用该区域内其他有实测资料的流域优选模型参数，并对优选的参数值求均值，之后再移用于该无资料流域;(b)通过地形特征推求模型参数;(c)应用区域回归分析推求模型参数;(d)直接移用临近的相似流域资料。本文选用第(d)方法，将黄河上游三关口流域的参数移用于临近无资料流域——好水川流域，检验CASC2D-SED模型的产流模块在无资料地区的适用性。

4.2 模型输入及参数率定

模型的输入数据包括降雨、实测流量、流域90 m×90m DEM、河道数据(图4)、土地利用(图5)以及土壤类型资料。土地利用类型有:林地、水面、耕地、草地。流域内土壤为黑垆土。流域内有2个雨量站:张银、杨河，15个计算点:阴山、兰家湾、范湾、上岔、张家台子、后沟、下老庄、团结、后海子、张银、岔口、下岔、何家岔、阳山庄、老张沟。

模型时段步长选取2 s，降雨资料的输入时段长0.5 h，实测流量资料的时段长及模拟流量输出时段长均定为0.5 h。

由于好水川流域与三关口流域多年平均降水量均超过500 mm，降雨年内分布均有明显的典型大陆性特点，主要集中在6—9月，且两流域均位于六盘山地区，植被相近，流域内土壤均为黑垆土，所以好水川流域的参数移用三关口流域的参数。

4.3 模型模拟结果

好水川流域缺乏径流观测值，但有长期的降雨资料。根据杨河站及张银站2010—2011年4—9月的实测降雨资料，采用CASC2D-SED中的产流模块进行降雨径流模拟，得到15个计算点的流量过程，并计算得到日来水量。由于好水川流域无实测径流资料进行对比，故将模拟计算得到的日来水量交由相关部门验证。经验证，模拟计算得到的日来水量是合理的，说明将三关口流域的参数移用于好水川流域后，CASC2D-SED模型的计算结果是可靠的。

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

从以上结果看，CASC2D-SED模型在三关口流域的产流模拟中计算精度良好，将其移用于相似流域好水川流域时，模拟结果也是合理的，说明CASC2D-SED模型可用于无资料地区的径流模拟。

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

5 结 语

随着地理信息系统、遥感技术以及气象卫星的进一步发展，资料的空间分辨率将会更高，可以更准确地反映流域的地形、土地利用、土壤和河道资料，分布式水沙物理模型将会得到更好的数据支持;在解决资料问题后，随着计算机技术的进一步发展，模型自身结构的不断优化以及模型算法的不断完善，分布式水沙物理模型CASC2D-SED的运行速度将不断提高，模拟精度将会更加优良且稳定。

［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，FOSTER 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］曹文洪，祁伟，郭庆超，等.小流域产汇流分布式模型［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］李致家，胡伟升，丁杰，等.基于物理基础与基于栅格的分布式水文模型研究［J］.水力发电学报，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.