物理成因产沙模型研究中亟待解决的几个问题

摘要：剖析了HUM-1、THU和CSU三个国内外有代表性的物理成因产沙模型的结构、特点和存在的问题，分析了物理成因产沙模型应该考虑和解决的主要水沙物理过程，指出该类模型当前急需解决的六个问题和基本思路。

关键词：物理成因产沙模型 模型结构 泥沙输移 参数区域规律 水土保持

1 引言

流域侵蚀产沙模型是对流域内土壤侵蚀和泥沙输移现象的概化和近似表达，包括统计模型和物理成因模型两大类。前者基于大量的实测气象、水文、泥沙观测资料，经统计分析得到。又可分单位线法、经验方程和随机模型等，其中，以经验方程形式多样，研究和应用最多；后者以流域上实际发生的水沙物理过程为基础，用一个或多个数学方程加以描述并用一定的数值方法加以求解，参数具有物理意义。又可分为集中性模型和分散性模型。 经验性模型是根据实际资料应用统计相关的方法建立起来的，其结构简单，使用方便，在制定公式使用的资料范围内有足够的精度，在实际生产中得到了广泛的应用。但它们缺乏充分的物理基础，外延效果差，在作地区移用和向设计条件延伸时精度难以控制，不能反映侵蚀产沙的时空变化过程和人类活动影响后所发生的变化。物理过程模拟模型试图对流域内发生的侵蚀产沙过程进行概化和近似，并用数学方程描述流域上侵蚀产沙的主要物理过程，再用比较严格的数值解法计算水沙运动过程。模型的物理基聪强，外延精度高，有利于地区移用和向设计条件延伸，可以模拟侵蚀产沙的时空变化，并可通过参数反映人类活动影响后水沙的变化。 鉴于此，物理成因产沙模型的研究与应用为大家所重视，国内外有关学者都在探索研制当中，近20年来发展很快，代表性的模型有美国的CSU、ANSWERS、FESHM、ARM和中国的HUM-1、THU等模型。本文着重剖析了HUM-1、THU和CSU三个典型的分散性物理成因模型，通过分析比较，指出当前物理成因模型研究中亟待解决的若干问题。

2 典型模型的剖析

根据反映模型结构的参数不同，物理成因模型又分为集中性模型和分散性模型。如美国的Stanford、ARM等模型都可归结为用许多模型参数来处理每一个单元上有意义的那一部分面积，考虑的因素较多，较全面，参数多，成因性更强。而CSU、HUM-1、THU和ANSWERS等模型是分散性的，它们把流域划分成许多较小的均匀元素，在每一个均匀子流域上都有反映影响水文和侵蚀过程诸如土壤、植被、土地管理等特征的参数，它可以分过程考虑，既有物理成因，概念明确，又有适用性，对于大而复杂的流域，其应用显示出独特的优越性，是物理成因模型研究的发展趋势。 2.1 HUM-1模型[１] HUM-1是河海大学1号模型的英文缩写，于80年代中期由陈国祥和汤立群提出，以后经在不同尺度的流域上应用，并不断加以补充与完善，逐渐发展成为一个物理过程清晰，成因性强的小流域产流产沙动力学模型。 2.1.1 模型结构 由水文模型和泥沙模型两部分组成。水文模型包括产汇流计算；泥沙模型包括产输沙计算。 水文模型把各单元分成透水面积（1-IM）和不透水面积IM，在不透水面积上降雨扣除蒸发后无其它损失，得径流R2；在透水面积上，通过模拟蒸发和下渗等降雨损失，计算产流R1，两者相加即为总径流量R。用运动波方程分别描述各单元的坡面和沟道汇流过程，采用四点隐式有限差分格式进行数值计算，可得单元和总出口断面的水文过程。 泥沙模型用坡面土壤侵蚀率公式计算坡面上的面蚀（包括细沟侵蚀）量和沟蚀（包括重力侵蚀）量，用沟槽侵蚀率公式计算沟道侵蚀量，两者之和便得单元流域总的土壤侵蚀量。然后，将各单元流域的土壤侵蚀量演算到出口断面得流域产沙量。 2.1.2 水文模型基本方程 (1) 产流计算 用Horton下渗公式计算超渗产流量，考虑下渗能力在各单元透水面积上分布的不均匀性，采用随机分布的B次方抛物线加以概化，则可按（1）式计算产流

(1)

式中 I为雨强；A是单元面积上点的最大下渗能力；为时段平均下渗率；R为时段产流量。

(2) 坡面及沟道汇流计算 黄土地区坡面很陡，洪水波的传播速度快，沿程坦化小，具有运动波的传播特征。黄土地区小流域的沟道比降较大，沟道水流过程仍可用运动波方程组来描述，这样，坡面及沟道水流运动都用运动波形式的圣维南方程组来描述。动量方程与连续方程联立求解，采用Preissmann四点隐式差分格式，可推求出水流运动差分方程

(2)

式中θ为权重系数，在[0.6，1]取值；K为系数；α为指数。对坡面流Q为单宽流量，R(x，t)为净雨过程；对于沟道断面水流，Q为流量，R(x，t)为旁侧来流。 根据（2）式可推求出任意时空不均匀降雨的坡面单宽流量过程和沟道任意点的流量过程。 2.1.3 泥沙模型基本方程 该模型将单元流域概化为“一本打开的书”，“书面”的交线代表单元流域的沟道，左右两面对称，每面上由不同坡度的梁峁坡和沟谷坡组成。从上游至下游，梁峁坡、沟谷坡及沟道的土壤侵蚀率Er、Eg、Ec可由下列各式表示

（3）

（4）

（5）

式中 Ar、Ag、Bc为系数；br、bg分别为梁峁坡和沟谷坡宽度；γm、γs分别为浑水容重和泥沙密实干容重；g为重力加速度；τ0、V分别为各区坡面或沟道水流的切应力和平均速度；为斜坡上泥沙起动切应力。

由此可得到单元流域土壤侵蚀率ET为

ET=2(Er+Eg)+Ec

（6）

2.1.4 全流域产沙量 坡面及沟道的侵蚀物质必须输移到出口断面才能成为流域的产沙量，侵蚀物质的输送受到水流条件、泥沙条件、沟道边界条件和流域尺度等多种因素的制约，尤其以流域尺度的影响反映最直观。一般来讲，在不考虑其它因素条件下，流域越小，泥沙输移比就越大。模型研究的流域尺度均属小流域，假定泥沙输移比为1.0，完全能够满足计算精度要求，并且可使计算大大简化。这样，流域出口断面的产沙量可由各单元土壤侵蚀量错开若干个传播时段叠加求得。 2.2 THU模型[2] THU于90年代初由清华大学谢树楠等研制，主要应用于黄土丘陵沟壑区的暴雨产沙模拟。 2.2.1 模型结构 由产流模型与产沙模型组成。产流模型用径流系数法，考虑降雨量、降雨强度和前期影响雨量对它的作用，通过回归分析确定径流系数；产沙模型假定坡面流为一维流动，压强按静水压强分布，流动中的动量系数为常量，不考虑泥沙的粘性，沟道泥沙输移比为1，根据水流连续方程、动量方程、泥沙连续方程和挟沙力公式，结合水文、气象、土壤和地质地貌学的基本原理，导得流域侵蚀产沙量计算公式。 2.2.2 产流计算公式 次暴雨过程中，单位面积单位时段内产生的径流深R用下式计算

R=f·P

(7)

式中 P为时段降雨量；f为径流系数，与降雨量P、降雨强度I和前期影响雨量Pa有关

{

f=0.002693P0.3814I0.9023

I>10mm/h

(8)

f=0.0006017P0.9051I1.44Pa0.1442

I<10mm/h

2.2.3 产沙计算公式 坡面产沙由（9）计算

(9)

沟道产沙由（10）计算

(10)

时段产沙量为

Wsi=Wsgi+Wsbi

（11）

单元面积次暴雨总产沙量为

（12）

以上各式中 Wsi为时段暴雨产沙量；Ws为次暴雨总产沙量；CA为地表裸露率；CE为侵蚀因子；D50为土壤中值粒径；L为单元坡面长度；S0为单元坡面比降；Sob为坡面平均比降Sog为沟道平均比降；SBb为某一土壤类型坡面所占比例；SBg为某一土壤类型沟道所占比例；Δti为计算时段长；A为单元面积。 2.2.4 全流域产沙量 THU模型根据研究流域的地理地貌特别是地质特性，按自然水系将其分为若干个区，每个区又分成多个单元，假定泥沙输移比为1，则全流域产沙量只要把所有单元产沙量累加即可。 2.3 CSU模型[3] CSU模型全称叫柯罗拉多州立大学模型，于70年代末由Ruh-Ming Li等人提出。该模型可以模拟流域表面水文过程、产沙过程和小流域上的水沙运动和时空变化，是一个物理成因较强的分散性模拟模型。 2.3.1 模型结构 模型分成坡面漫流部分和河道系统部分。坡面漫流部分可以模拟截留、蒸发、填洼、下渗等损失及雨滴溅蚀和面流冲刷等过程，并将坡面水沙演算到最近河槽；河槽系统部分模拟坡面水沙在河槽内的运动，确定槽泥沙的冲淤。 2.3.2 产流计算 时段平均下渗率用（13）式表示

（13）

ΔF=-[(2F-KΔt)]/ 2+{[(2F-KΔt)2+δKΔt(δ+F)]/1/ 2}/2

（14）

式中 Δt为时段长；ΔF为下渗改变量；F为下渗量；K为水力传导度；δ为势头参数。 产流量计算如下

(15)

2.3.3 汇流计算 坡面及沟道汇流均用一维运动波方程组来描述

(16)

式中 Q为流量；A为断面面积；ql为超渗雨或旁侧来流；S0为坡面或沟道坡度；f为阻力系数；R为水力半径；g为重力加速度。 方程组（16）采用有限差分方法联立求解，即可求得坡面或沟道内任意位置的流量过程及其它水力要素。 2.3.4 产沙计算 流域产沙量计算遵循以下产沙模式：假定雨滴溅蚀量为Za，供沙量为Sa，Sa=Za，水流挟沙能力为GT，径流冲刷量为D，则

{

WT = GT

GT≤Sa

（17）

WT = Sa+D

GT＞Sa

Za=DrΔt

(18)

(19)

D = - Df(ΔZP+Za)

(20)

GT = B(gb+gs)

(21)

Gb =a(τo-τc)b

(22)

gs =(gb/11.6)[AZ-1/(1-A)Z ]{[(V/u*)+2.5)]I1+2.5I2}

(23)

式中 Dr为雨滴溅蚀率；Δt为时段；I为雨强；a5、b5分别为系数和指数；Zw是水层厚度与松散土层厚度之和；Zm为3倍的雨滴中值粒径；Cg为地面植被覆盖度；Cc为树冠覆盖度；ΔZ p是时段内径流冲刷量；Df为与土质有关的系数；B为坡面或沟道宽度；gb为推移质挟沙能力；gs为悬移质挟沙能力；a、b分别为系数和指数；A为无量纲系数；Z为悬浮指数；V为坡面或沟槽平均流速；u*为摩阻流速；I1、I2为Einstein积分。 2.3.5 全流域产沙量 CSU模型按自然水系划分流域，每个子流域被概化成具有一个坡面的对称的“open book”，各子流域产沙量按（17）式计算，所有子流域产沙量叠加得全流域产沙量。 2.4 模型评价 通过对上述各模型结构及计算方法的介绍，可以看出各模型有以下特点：（1）模型都考虑分单元并对各单元进行概化；（2）考虑降雨分布的不均匀性，各单元都有其代表性的雨量过程；（3）产流计算方法相近，模拟效果均较满意；（4）HUM-1和CSU模型考虑了坡面及沟道汇流计算，所得结果为时空变化过程。THU模型未考虑汇流计算，所得结果为时间变化过程；（5）CSU模型单独模拟降雨溅蚀过程，HUM-1模型将雨滴溅蚀考虑在梁峁坡土壤侵蚀计算中，THU模型将之列于流域总侵蚀量计算中；（6）HUM-1、CSU模型都考虑了坡面径流侵蚀产沙和沟道侵蚀产沙计算，THU模型将之都考虑在总量计算中；（7）三个模型在实测资料验证的基础上，都较好地应用于产流产沙量预报计算。其中，HUM-1、THU模型进一步用来计算水土保持减水减沙效益，分析水沙变化之原因。 同时，各模型存在的主要问题有：（1）成因性越强的模型，对基本资料的要求也越高，如HUM-1、CSU模型要求有连续记录的降雨过程资料，THU模型要求有详细的地形地貌资料，而现有的水文资料和下垫面资料均较粗，给模型的应用带来困难；（2）有些重要的侵蚀产沙过程还没有单独描述，独立模拟，如HUM-1模型将雨滴溅蚀过程隐含于梁峁坡坡面侵蚀之中，重力侵蚀隐含于沟谷坡坡面侵蚀之中，CSU模型没有专门描述并模拟重力侵蚀过程，THU模型从整体上描述侵蚀产沙过程，不区分各种侵蚀过程，因此，模型的微结构还不尽合理；（3）三个模型都没有很好地考虑沟道内泥沙的沉积与再输移问题。其中，HUM-1和THU模型都直接假定泥沙输移比为1，以简化计算。这种简化计算带来的误差基本上靠参数的拟合来弥补；（4）模型在各种侵蚀产沙过程没有合理描述前提下之所以都有较好的模拟结果，主要是因为产沙模型各方程都带有若干系数，而这些系数又是通过流域实测产沙量资料率定的，这在研究方法和技术途径上都是允许的，但带来的问题是不同地区不同的研究对象，其参数可能是不一样的，这就给模型的移用带来不便；（5）各模型的应用面积均还很小，除THU模型已用于3240km2的皇甫川流域外，HUM-1模型最大的应用面积为41.2km2，而CSU模型只有1.27km2，如何推广应用于较大流域的水沙预测，仍然有许多问题值得探讨；（6）HUM-1、THU模型都能用来还原计算人类活动影响的产流产沙量，进而计算出水土保持总的减水减沙效益，这是比较宏观的分析结果。但在小流域规划设计中，需要了解各项水土保持措施的减水减沙效果，现有的模型似乎都不能进行水保效益的分离评估。此外，进一步提高模型的计算精度，克服累计误差也是值得探讨与改进的问题。

3 亟待解决的几个问题

鉴于以上分析，归纳起来物理成因模型研究中有6个方面的问题需要深入分析研究，进一步加以完善：(1)建模所需资料的获取问题；（2）重力侵蚀过程的模拟；（3）沟道泥沙输移计算模式；（4）参数区域规律；（5）模型的推广应用；（6）水土保持措施效益的分离评估。以下就这6个方面问题分别进行讨论。 3.1 建模所需水沙资料的获取 由于种种原因，现有的流域水文泥沙基本资料都比较粗，给模型参数的率定和应用带来困难。解决的方法之一是对模型结构加以改进，使之适用于现状资料的产流产沙计算；另一办法是将现状资料进行加工处理，重新形成水文泥沙资料，使之适应于现有的模型，但有一定难度。 3.2 重力侵蚀产沙的模拟 黄土地区的坡面尤其是沟谷坡，重力侵蚀发生频繁，目前还缺乏有效的解决方法。主要原因是很少有描述重力侵蚀发生过程与数量的观测资料，缺乏描述重力侵蚀的基本条件，因而不能把过程分离出来进行定量描述。改进的方法应加强资料的观测。传统的方法可以用打桩法进行分析研究，也可用高速摄影法或同位素测定法加以研究，甚至可利用卫星或航片资料进行分析，获得第一手的重力侵蚀过程资料，再探索用数学描述的可能性。重力侵蚀的数学表达应将水动力学、土壤学、土力学等基本理论与方法结合起来加以研究。同时也可借鉴滑坡、泥石流方面的研究成果。 3.3 河道泥沙输移计算 河道泥沙输移计算的关键是确定泥沙输移比。目前国内开发的各类模型包括物理成因模型，都假定泥沙输移比为1，使问题处理得到简化。但从理论上讲，任何一个流域的泥沙输移比不可能为1，它受到多种因素的影响，如降雨径流条件，土壤类型及泥沙特性，地形地貌，地质构造运动和流域尺度等。对于特定的黄河中游地区，只有那些支毛沟的泥沙输移比才接近于1，面积越大，泥沙输移比就越小。具体表现为来自坡面的侵蚀泥沙进入沟道后存在淤积与再输移现象。 所谓泥沙输移比Drr是指流域某一断面实测输沙量Wr与断面以上流域侵蚀产沙总量T之比，可用下式表示

Drr=Wr/T

(24)

文献[4]认为，确定某一个区域的泥沙输移比有三个前提条件：（1）泥沙的粒径。河道中泥沙的最大粒级与最小粒级要相差百万倍，输移不同粒级的泥沙需要不同的环境条件，小粒径的物质可以在一般径流条件下运动，大粒径的物质一般水流条件就搬运不了；（2）时间系列。对于一个确定的流域，泥沙的输移能力从长系列看是稳定的，但短时间又是不稳定的。主要原因是输移泥沙的动力条件为径流，而决定径流的直接因素是降雨。降雨在年内分配是不均匀的，年际间也存在较大的差异，具有周期性。因此，泥沙输移过程也随水流动力条件具有周期性，只有在一个稳定的较长系列内分析泥沙输移比，才能得出符合实际的Drr值，经常用一年甚至多年；（3）空间范围。指的是流域的大小，流域越大其环境条件越复杂。一般情况下，中小流域的空间内环境因素相似性较大，而大流域环境因素相似性要差。从这一点上讲，国内外都得出了基本结论，即泥沙输移比与流域面积成反比。 解决这个问题的途径，一是根据沟道水流泥沙的动力条件，建立沟道泥沙输送函数，将来自坡面上的泥沙演算到流域出口，得出流域侵蚀产沙量，这是一种简单实用的处理方法。二是也可将一维水流泥沙冲淤数学模型引进到流域产沙模型，将来自坡面的水流泥沙作为一维数模的旁侧分散来流来沙，根据沟道边界条件、地形条件、初始条件和河床物质组成用有限差分求解，把水流泥沙运动过程演算到出口断面即得流域侵蚀产沙量。 3.4 参数的区域规律 当前，模型应用的流域还不够多，类型单一，若应用于同一类型区，模拟结果都能达到较高的精度，但移用于不同类型区时，参数都要重新率定，因此，存在着参数区域规律问题。要解决这个问题，须做大量耐心细致的研究工作，将模型在不同类型区的各种不同流域上对水沙进行模拟，并对参数分布进行分析研究，找出其规律性，这样，当模型要用来解决实际问题时，可立即得到参数值，不需要再进行参数率定，提高使用效率，扩大适用面。 3.5 模型的推广应用 成因模型大多是在小流域上建立起来的。由于该类模型对基本资料的要求相对较高，而小流域又易获得所需资料。因此，目前普遍应用于小于100km2的流域，效果比较好。对于中等或大流域，由于缺乏降雨、径流、泥沙过程资料，控制密度又较疏，应用还不够理想，必须在现有模型基础上加以改进，使其适合于现状动力条件下的中大流域产流产沙量的计算。另一方面，实用型的物理成因产沙模型大多是分散性的，如本文所述各模型，它们的最大特点是除分过程分别描述水沙过程外，流域按自然水系划分单元，因此，从理论上讲，单元流域可划分到任意小，实用起见，单元流域通常划分到能获得水沙资料的那一级支流，各单元通过“搭积木”的方法进行演算，便可得任意大小流域的产流产沙量，也就是说，物理成因模型的分散性结构为推广应用奠定了基础。 3.6 水土保持措施减水减沙效益分离评估 水保效益的评价方法主要有水文法和水保法。近年来作者用模型法评价水保效益取得了较好的效果，但仍属总量评价[5]，没有对每种措施加以分离评估，这也是水保效益评价的难点。 笔者近期通过国家自然科学基金重大项目的研究，试图找到用模型对水土保持措施减水减沙效益加以分离评估的方法，其主要思想是通过参数的变化来反映人类活动影响下的产流产沙过程，从而分别计算各项措施的效益，模型参数可通过坡面小区或微小流域加以率定。此问题解决后，可直接指导小流域的水土保持规划设计，为流域治理提供决策依据。

4 结论

分散性物理成因产沙数学模型是当前流域泥沙模拟研究的发展趋势，具有过程明确，参数少，概念清晰，使用方便，适应面广的特点。本文列举的三个分散性模型都得到了较好的应用，具有实用价值和较广阔的推广前景。但在基本资料的获取、模型微结构、泥沙输移计算、参数区域规律、模型推广应用、水土保持效益的分离评估等方面还需要加以改进完善。

参 考 文 献

[1] 汤立群，陈国祥。小流域产流产沙动力学模型。水动力学研究与进展，A辑，Vol.12，No.2，1997.

[2] 陈国祥，谢树楠，汤立群。黄土高原地区流域侵蚀产沙模型研究。黄河水利出版社，1996.

[3] Ruh-Ming Li， Water and Sediment Routing from Watersheds ，Published by John Wilely & Sons ， Inc. 1979.

[4] 唐克丽，熊贵枢，等。黄河流域的侵蚀与径流泥沙变化。中国科学技术出版社，1993.

[5] 汤立群，陈国祥。物理概念模型在水保效益评价中的应用。水利学报，1998，(9).