(中国水利水电科学研究院；清华大学 水电系；内蒙古水利学校)

1 皇甫川流域基本概况

　　皇甫川流域位于黄河中游龙门一河口区间的上段，是黄河粗泥沙产区及泥沙来源的主要地区，年平均向黄河输送约0.5亿吨的泥沙。主流长125km，流域面积3240km。皇甫川流域水系主要由干流纳林川和支流长川组成。按不同地表物质分布和侵蚀差异可将流域划分为三个类型区: (1)黄土丘陵沟壑区。主要分布于流域的东部和西南部，例如十里长川以东和中间部位，面积为1627km²，沟壑密度为5～9km/km²，本地除部分梁峁和缓坡地为耕地外，多为天然草场，植被度为20%左右。本区黄土较厚，约20～30m之间，是陕晋黄土的边缘，呈现较典型的黄土梁峁和黄土沟谷地貌。该区以水蚀为主，水蚀、风蚀和重力侵蚀交错,(2)沙土丘陵沟壑区。主要分布于纳林川中下游以东到长川以西地区和库布其沙漠边缘。面积为243km²，平均沟壑密度为4.2km/km²，该区地势较平缓，地形较完整，沟道较浅。表层原为黄土，由于长期开荒及过度放牧的结果，天然植被严重破坏，经风蚀，引起沙化。该区的主要植被有沙蒿、柠条和沙棘等，覆盖度为15%左右。该区水蚀较轻，风蚀为主要侵蚀方式。(3)砒砂岩丘陵沟壑区。主要分布于流域的西北部，例如干察板沟，圪秋沟，虎石沟和尔架麻沟。面积为1379km²，沟壑密度平均为7.42km/km²，该区水土流失极其严重，地形切割十分破碎，坡陡沟深，本区覆盖度很低，基岩大面积外露。该区侵蚀以水蚀为主复合重力侵蚀。
　　流域内降雨的年际和年内变率大，年内降雨集中于6～9月,占全年的80%以上。且常以暴雨出现，历时短，强度大，一次降雨往往集中在几分钟至几十分种分内。对流域面雨量的分析计算，多年平均雨量为379mm，1959年年雨量高达685mm，1965年只有192mm，两者相差三倍多。
　　皇甫川出口站皇甫站多年平均径流量为1.7亿m³，输沙量为5676万吨，多年平均含沙量320kg/m³。径流和输沙量在年内分配极不均匀，沙量尤甚，主要集中于汛期(6～9月)。汛期沙量占全年沙量的98.9%，汛期水量占全年水量的82.1%。流域历年最大洪峰流量高达1.16万m³/s，洪水均为暴雨形成。皇甫川经常出现高含沙洪水，最大含沙量为1270km/m³ 。洪水平均含沙量达422km/m³。悬沙粒径大于0.05mm，占54.5%，大于0.05mm粒径的沙量多年平均为3093万吨，皇甫站多年平均输沙模数为1.78万t/km²·a。

2 次暴雨产流产沙计算方法研究及模型建立

　　研究流域产流产沙规律，需要遵循成因分析的途径，在泥沙运动力学理论的基础上，结合水文学、气象学、土壤学和地质地貌学的基本原理，研究包括降雨、地形、土壤、植被、工程等诸多因素的影响和流域产流产沙的物理图形及定量关系。只有在此基础上建立流域次暴雨产流产沙的数学模型，才能实际地反映流域的产流产沙过程和水沙变化。流域产流产沙过程非常复杂，其影响因素很多，有气候因子（降雨量、降雨强度、蒸发等）、下垫面情况（土壤性质、地形、植被等）和人类活动影响等。
　　研究未经治理情况下的流域产流产沙规律，需要突出降雨和产流产沙之间的矛盾。作者把皇甫川流域根据具体情况分为三个计算时段：1966至1969年作为第一计算时段，这四年期间流域基本未经治理或水土保持措施尚未发挥作用，作为公式推导和验证；1970至1979年为第二计算时段；1980至1989年为第三计算时段。第二和第三计算时段作为公式应用。根据“黄河皇甫川流域土壤侵蚀系统模型和治理模式”一文中已有的资料，将皇甫川流域划分为三个小流域，即沙圪堵流域、长滩流域和皇甫区间（即沙圪堵流域和长滩流域以下区域）；在此基础上，又划分为14个子流域；每个子流域再划分为若干个类型不相同的单元，共510个单元，而每个单元作为模型的基本计算单元。
2.1 建立次暴雨产流模型
　　对本流域1963～1973年85场降雨径流分析研究，取Ie=0.2mm/min和Pe=6.0mm作为产流产沙的临界值，即认为降雨开始和结束时雨强小于0.2mm/min和雨量小于6.0mm时不产流产沙，忽略这两部分降雨，而取中间的一段时间作为产流历时，取产流历时内的降雨量和降雨强度分别作为该场暴雨的有效降雨量和有效降雨强度（以下简称降雨量和降雨强度）。

当I<10mm/hr时

当I>10mm/hr时

　　公式(2)中没有考虑前期雨量的影响，这是因为对于高强度降雨，前期雨量的影响较小。
2.2 产沙模型的建立
　　关于产沙公式，到目前为止已有许多成果，但仍存在许多问题，国外在这方面的发展较国内迅速，但直接用于我国黄土高原产沙计算是不合适的。国内产沙公式大多是从研究黄土高原上的实验数据或观测资料回归分析得到的经验公式，对于研究产沙的物理过程方面主要是在90年代，仍带有深刻经验烙印。为此，皇甫川流域产沙模型在谢树楠等的产沙模型基础上，结合水文、气象、土壤和地质地貌学的基本原理，着重解决了流域的不均匀性和K值率定等问题建立的。
2.2.1 基本假定
　　(1)暴雨产生的径流为坡面一维流动；(2)压强按静水压强公布；(3)流动中的动量系数为常量；(4)泥沙不考虑粘性；(5)沟道泥沙输移比为1。
2.2.2 公式推导
　　(1)挟沙力公式 根据爱因斯坦挟沙力的概念，水力参数Θ和输沙参数Φ的关系可以写为

将(3)式和(4)式代入(5)式，整理合并为

式中 β₁=αg^-β+1/2/{[(γ_s-γ)/γ]g}^β+1/2D^-β+3/2；q_s为输沙率；γ_s,γ为分别为泥沙和水的比重；g为重力加速度；τ₀为泥沙的中数粒径；α为常数；β为指数。
　　由式(6)可见，坡面输沙率主要与作用在坡面的水流剪切力成正比。
　　（2）泥沙输移的连续方程

式中P_s为单位面积上的起冲率；x为在流动方向上沿坡面的距离。
　　（3）水流的连续方程

式中 q为单宽流量；q₀为斜坡上段流入的单宽流量; q_*为单位时间、单位面积斜坡上由降雨产生的径流量。

　　（4）水流运动方程 对于明渠，具有侧向入流的水流运动方程

式中 y为水深；s₀,s_f为坡面坡度和能坡；β₂为动量因子；U为进入主流中的侧向水流的流速(即雨点的终极速度)；φ为x方向与速度U的夹角；θ为x方向与水平方向之间的夹角
　　（5）泥沙起冲率与边界剪切力 泥沙的起冲率可以表达为边界剪切力的函数。即

式中 a为描述坡面土壤侵蚀特性的系数；b为指数。
将(6)和(7)联立求导，整理合并得

与(10)式比较得出

又∵τ₀=γys_f

其中 γ为水流容量； 为平均摩阻流速。

将(14)式代入(12)式得

　　（6）径流深的求解 一般说来，对于坡面流，水深很小，可近似认为是均匀流，则(9)式可简化为

　　φ为x方向与速度U的夹角，可以考虑降雨与径流垂直，即cosφ=0，则上式可进一步简化为

则径流深为

对（16）式微分，得

将(17)式代入(15)式，得

　　（7）地表降雨下的边界剪切应力 明渠的阻力可表示为

　　在坡面上的水流由于降雨的加入，其阻力系数F已不同于明渠阻力，此时降雨强度对决定阻力系数起重要作用。为解决这个问题，70年代初先是南朝鲜的Y.N.Yoon，后是李日明和沈学汶，都先后做过室内的水槽实验。试验表明：当雷诺数Re<2000时，雨强I越大，F就越大；当Re>2000时，雨强I对F不起作用。
　　将阻力系数分为两部分：一部分是没有降雨时的阻力系数F₀，另一部分是由降雨所增加的阻力系数Fr，认为Fr在光滑和粗糙面上不变。当Re在2000～12600时

由于黄河中游各流域坡度陡，故在计算阻力时可以按(19)式计算。

               ∵y=τ₀/γs₀

　　（8）有关各参数的确定α,β值的确定 谢树楠等分析了团山沟和段川14个径流试验(1966～1967年)的645组资料，建立了水力参数Θ和输沙参数Φ的关系

将(22)与(3)式比较，得

α=1.81; β=2.16

　　α,β₂值的确定 Yoon把变流速的运动方程写成

　　Yoon采用实验资料对等式右边各项的影响程度进行了分析，发现s₀对s_f的影响最为显著，第二项及第三项对s_f的影响几乎相同，而且均相当于s₀的1/10，最后一项可以忽略不计，据此可得

　　另据沈学汶及李日明的研究，当雷诺数大于2000时，s=1.05。根据以上研究，Komura建议采用s=1.20，由降雨产生的径流动量因子β₂值，可从Yoon的实验结果确定。Yoon等人研究了降雨强度，雷诺数和渠道坡度对动量因子的影响，认为β₂值随雷诺数和雨强的增大而减小，但坡度影响不显著。据他们的研究成果，取β₂=1.1。
　　其它参数值的确定
　　据日本学者的研究，认为K₂=0.6为宜。

　　（9）产沙公式 单位面积、单位时间的产沙量Ws可由下式表示

将(8)、(10)、(13)、(18)和(21)式代入(24)式，积分并整理得

若忽略坡顶来流q₀，则(25)式可简化为

把以上各参数代入(26)式，整理得

　　暴雨产沙不仅仅与上式右端各因子（即径流量、坡长、坡降和中数粒径等）有关，而且与土壤性质(S_B)、土壤侵蚀特性(C_E)和裸露率(C_A)有关。土壤性质、土壤侵蚀特性和裸露率的不同导致产沙量大小差异。据研究，Ws与S_B、C_E和C_A成正比关系。因此在式(27)中右边乘以S_B、C_E和C_A值，就组成了次暴雨产沙公式。
式中 q_*表示由降雨产生的单位面积径流量，单位是m³/m²·s，实际为单位时间降雨产生的径流深，可表示为：q_*=fI
故次暴雨产沙公式为

次暴雨总产沙量为

上面各式中各因子的含义及单位：Wsi为次暴雨产沙量(10⁴t);Ws为次暴雨的总产沙量(10⁴t);Wsp为坡面产沙量(10⁴t);Wsg为沟道产沙量(10⁴t);S_B为土壤种类所占比例(%);S_Bp为某一土壤类型坡面所占比例;S_Bg为某一土壤类型沟道所占比例;C_A为地表裸露率;C_Ap为坡面地表裸露率;C_Ag为沟道地表裸露率;C_E为侵蚀因子;C_Ep为坡面侵蚀因子;C_Eg为沟道侵蚀因子;D₅₀为土壤的中数粒径(mm);f为径流系数;I为降雨强度(mm/hr);L为坡面长度(m);Lp为单元坡面长度(m);Lg为单元沟道长度(m);S₀为坡面比降;S_0p为单元坡面平均比降;S_0g为单元沟道平均比降;A为单元面积(km²);Δt_i为计算时间步长(hr)。

3 模型中各参数的率定

3.1 径流系数f
　　本模型采用的径流系数f是根据本流域降雨、径流资料，按有效降雨量、有效降雨强度和前期影响雨量进行回归分析所得到的公式：f=0.0006017P^0.91I^1.44P_a^0.14， (I<10mm/h_r)；f=0.0026930P^0.38I^0.90， (I>10mm/h_r)。
3.2 有效降雨强度I
　　流域产沙不仅取决于降雨量，而且更主要取决于降雨强度及分布。因此，对于产沙的降雨强度不能笼统地用次暴雨的有效降雨量除以有效降雨历时。
　　采用谢树楠等研究出的降雨强度分布模拟公式：I(t+Δt)={[(t+Δt)^0.366-t^0.366]/Δt}PT^-0.366
3.3 有关流域下垫面因子的确定
　　下垫面因子的确定依据地形图和实际调查资料。
　　地形图由中国测绘资料信息中心提供，由中国人民解放军总参谋部测绘局1971年11月航测、1977年5月调绘，1979年第一版。比例尺1∶50000，等高距10m。
　　流域总面积A及单元面积A_k: 采用工具是求积仪，并结合米格法确定面积。应用求积仪时，顺时针量测一圈，再逆时针量测一圈，取其平均值，最后平差，即为面积。
　　主沟长度L和单元坡长L_k: 主沟长度L是从流域出口沿主沟道一直延伸至接近分水岭。坡长是坡面的平均长度，对于单元坡长L_k来讲，L_k是定值。
　　单元沟坡坡度S: 对于流域内不同级别的沟道其沟坡平均坡度是不同的,即使在同一沟道内，上、中、下游各部位的沟坡坡度也有差别。计算产沙时，若采用平均坡度计算，势必带来极大的误差。因此，坡度要采用单元沟坡坡度，这样坡度是均匀的。
　　裸露率C_A: 裸露率C_A就是单元内无植被的面积和整个单元面积的比值，由实测资料知单元的植被覆盖率C_B，则C_A=1-C_B。
　　土壤中数粒径D₅₀:D₅₀由该流域表土的级配曲线率定，不同的土壤类型其D₅₀不同.
　　侵蚀因子C_E: 侵蚀因子C_E是反映流域侵蚀程度的综合系数,它主要决定于流域坡面土的缔固程度、粘土比、分散比及沟壑密度等因素。
　　本流域采用1970年以前(无水保工程前)的多年平均侵蚀模数W_s0作为确定各个单元的侵蚀因子的根据，即：C_E=KW_s0
　　研究皇甫川三个流域59组944场暴雨产沙的K值变化规律，分析结果表明，流域的K值系数为0.00175。

4 模型验证

　　为了验证产流产沙模型的适用性，采用1966～1969年的59组944车雨资料进行验证。三个流域逐次计算出每次暴雨的产流产沙量，再与实测值比较。 沙圪堵流域产流产沙量的计算值与实测值比较，相关系数分别为0.898和0.995；长滩流域和皇甫区间的计算值与实测值比较，相关系数分别为0.752和0.931；皇甫川流域的计算值与实测值比较，相关系数分别为0.967和0.962 (见图2和3)。这三个流域的共同特征是总径流量或总产沙量计算值和实测值相差很小，相关性较好(只有长滩流域和皇甫区间差一些，这是由于其实测值是全流域的径流量减去沙圪堵流域的径流量得出的原因)，说明模型精度较高，可以作为本流域产流模型应用。（注：图中Qc表示计算产流量值，Qm表示实测径流量值；W_sc表示计算产沙量值，W_sm表示实测输沙量值。

5 产流产沙模型应用及水沙变化原因分析

　　为了分析流域的水沙变化原因，采用1970～1989年共213组3408场暴雨，逐单元、逐次、分时段进行次暴雨产流产沙计算，然后求其总和，即某次降雨的产流产沙量。年产流产沙量由次暴雨累加，产沙量的计算按坡面、沟道产沙及砒砂岩、黄土和沙土等三种土壤类型分别计算。在分析水沙变化原因时，按三大区即沙圪堵站以上流域、长滩站以上流域及皇甫区间和全流域分别统计。

5.1 应用产流模型分析水量变化原因
　　沙圪堵流域在进入70年代以后：（1）降雨影响水量变化：70年代径流量增加了21857.1万m³，年均2185.7万m³，占39.09%；80年代径流量减少了5335.9万m³，年均533.6万m³，占9.54%；1970～1989年径流量增加了8260.6万m³，年均826.1万m³，占14.77%。（2）人类活动和其他因素影响水量变化：70年代径流量减少了184.8万m³，年均18.5万m³，占0。33%；80年代径流量减少了3179.9万m³，年均318.0万m³，占5.69%；1970～1989年径流量减少了1682.3万m³，年均168.2万m³，占3.01%。详见表1。说明该流域70年代和80年代人类活动和其他因素影响（主要是水保工程作用）水量减少甚微。
　　长滩流域和皇甫区间在进入70年代以后：（1）降雨影响水量变化：70年代径流量增加了5474.5万m³，占7.45%；80年代径流量减少了13603.6万m³，占18.51%；1970～1989年径流量减少了4064.5万m³，占5.53% 。（2）人类活动和其他因素影响水量变化：70年代径流量减少12037.8万m³，占16.38%；80年代径流量减少8658.0万m³，占11.78%；1970～1989年径流量减少了10347.9万m³，占14.08%。详见表1。说明该流域人类活动和其他因素影响（主要是水保工程）的减水作用较显著。
　　因此，整个流域进入70年代以后，有一定的减水效益：（1）降雨影响水量变化：70年代径流量增加了27331.6万m³，年均2733.2万m³，占21.12%；80年代径流量减少了18939。5万m³，年均1894万m³，占14.64%；1970～1989年径流量增加了4196.1万m³，年均419.6万m³，占3.24%。（2）人类活动和其他因素影响水沙量变化：70年代径流量减少了12222.6万m³，年均1222.3万m³，占9.45%；80年代径流量减少了11837.9万m³，年均1183.8万m³，占9.15%；1970～1989年径流量减少了11880.3万m³，年均1188.0万m³，占9.18%（详见表1）。说明皇甫川流域进入70年代以后，人类活动和其他因素影响（主要是水保工程）水量有一定的减水作用，但不明显。沙圪堵流域水保工程较少，减水作用很小；长滩和皇甫区间由于有较多的水土保持工程，使得减水效益较高。
5.2 应用产沙模型分析沙量变化原因
5.2.1 年降雨产沙量的计算方法
　　次暴雨的产沙量总和基本上代表了年产沙量。沙圪堵流域汛期实测输沙量占年输沙量的94.2%；皇甫川流域汛期实测输沙量占年输沙量的98.9%。故采用由次暴雨产沙累积求年产沙量的办法。

　　 注:（5）=（4*）-（3）,（6）=（5）/（4*）,（7）=（3）-（4*）,（8）=（7）/（4*）

5.2.2 结果分析
　　70年代，沙圪堵流域坡面产沙7773万吨，占总产沙量的29.4%，沟道产沙18683万吨，占总产沙量70.6%；而砒砂岩类产沙11383万吨，占总产沙量的43%，黄土类产沙13774万吨，占总产沙量的52.1%，沙土类产沙1297万吨，占4.9%(参见表2)。70年代实测总沙量比60年代增沙13969万吨，其中降雨因素增加6708万吨，占48.0%，而人类活动和其他因素增沙7262万吨，占52.0%(参见表3)。80年代，沙圪堵流域坡面产沙4483万吨，占总产沙量的28%，沟道产沙11527万吨，占72%；而砒砂岩类产沙6767万吨，占总产沙量的42%，黄土类产沙8333万吨，占总产沙量的52%，沙土类产沙910万吨，约占6%。80年代比60年代增沙390万吨，其中降雨因素减沙3738万吨, 人类活动和其他因素增沙4128万吨。80年代沙圪堵流域增沙主要是由于人类活动和其他因素引起的。说明沙圪堵流域增沙的趋势应引起足够的重视。
　　长滩流域及皇甫区间70年代该流域坡面产沙5028万吨，占总产沙量的15.1%, 沟道产沙28349万吨，占84.9%，而砒砂岩类产沙10432万吨，占31.3%，黄土类产沙20924万吨，占62.7%，沙土类产沙2020万吨，占6.1%(参见表2)。70年代比60年代减沙1224万吨，其中降雨因素增沙3779万吨，而水保工程减沙5002万吨，说明水保工程的作用是显著的(参见表3)。80年代该流域坡面产沙4088万吨，占总产沙量的14.6%，沟道产沙24003万吨，占85.4%。砒砂岩类产沙8784万吨，占31.3%，黄土类产沙17559万吨，占62.5%，沙土类产沙1747万吨，占6.2%。80年代比60年代减沙8350万吨，其中降雨减沙1507万吨，而水保工程减沙6843万吨。说明80年代水保工程的减沙作用比70年代更好。