近年来，山洪灾害已成为全球范围内致死率最高的一类洪涝灾害，其损失排在各类自然灾害前两位的国家超过百个^[1-2]。中国山区面积广阔，山洪致灾问题尤为突出，2011—2017年间的山洪致死人数超过洪涝灾害死亡总人数的70%^[3]。山洪致灾机理与防灾减灾技术研究已成为国内外热点。

山洪不简单是山丘溪沟中的暴涨洪水，其汇集演进过程时常还伴随着剧烈的泥沙输移与河床冲淤。文[4-8]分别在北美、欧洲、非洲、南亚地区的山洪过程中观察到了明显的泥沙运动与河道变化。北京“2012.7.21”房山山洪(简称“7.21”山洪)、“2017.6.18”门头沟山洪、“2018.7.16”密云山洪过程中，同样发生了强烈的泥沙输移、边岸侵蚀现象。Rickenmann等^[9]则进一步揭示，山洪过程通常有“源头支沟产沙-中游剧烈冲淤-下游侧向侵蚀”的空间分布特征。

一般认为，泥沙输移在山洪冲击掏蚀、漫滩淹没、泥沙淤埋、阻塞沟道-级联溃决等山洪致灾过程中发挥着重要作用^[10]，近年的灾后分析愈加揭示出其重要性^[10-11]。Badoux等^[12]统计分析瑞士40年间的山洪灾害数据，发现推移质输移致灾损失可占到灾害总损失的35%以上。文[13-14]发现，泥沙输移及其导致的沟道形态变化对灾害结果的影响与水文要素同等重要。与之相对应，试验研究逐步明确了泥沙作用的力学机制。文[15-17]试验探索推移质的增阻壅水作用，先后建立了一系列经验公式以描述推移质作用下的水深、阻力特征。

当前的山洪灾害风险分析，通常仍以水文模型或水力学模型为工具开展洪水演进计算；而山区河流泥沙运动模拟的关注焦点则多为输沙量、河床形貌等泥沙特征。如Chiari等^[18]考虑床面形态阻力的影响，以一维模型SETRAC复演洪水过程中的沟道输沙特征；Viparelli等^[19]以活动层概念为基础，研究河床垂向沉积结构的变化；文[20-21]研究耦合边岸-沟道侵蚀过程，模拟拆坝后河床的长期演变过程。近几年，山洪演进与泥沙输移模拟相结合的洪水致灾模拟研究逐渐出现。如Radice等^[22]用一维模型BASEMENT计算山区沟道设计洪水水面线，表明考虑泥沙输移有助于提升山洪模拟准确性；Chen等^[23]得到了波罗电站山洪灾害的二维模拟结果，指出了泥沙输移对洪水位的增值作用。但是，关于泥沙输移与河床冲淤影响山洪最高水位、最大流速、峰现时间等洪水特征的系统研究还很少见。

本文以北京房山“7.21”山洪为例，模拟计算山洪过程中的床沙级配、床面高程、水流阻力等变化，探索泥沙输移对山洪水位、流速等洪水特征值空间分布、变化过程的影响机制。这对于深化认识山洪致灾机理、揭示泥沙输移效应、改进山洪防治技术等均具有重要价值。

1 红螺谷沟“7.21”洪水概况

京郊土石山区地势陡峻，地表风化强烈，沟道松散堆积物丰富，山前迎风区易出现局地性暴雨，陡涨陡落的山洪频繁发生。2012年“7.21”暴雨期间，北京多处站点降水超过建站以来的历史极值，引发了严重山洪灾害。房山处于暴雨中心地带，局地日最大降雨量460 mm，受灾人口超过80万人，经济损失超过50亿元，灾情尤为严重^[24]。灾后调查发现，区内多数沟道发生了显著冲淤变化，所调查的沟道有26 km主沟冲刷损毁、76.8 km沟道明显淤积抬升，估测冲淤量分别达7.89和87.3万m³。

红螺谷沟位于房山区周口店镇，是“7.21”洪水的重灾区，流域形态如图 1所示，基本参数见表 1。该流域内分布多处工厂、村落，沟道沿线人口密集，沟源区及干沟内积累了大量砂石堆积物。

“7.21”山洪中，该流域降雨过程持续约20 h，日降雨量高达370 mm，灾后调查估算推得新黄路与宝涞路交汇处的河道断面洪峰149.4 m³/s，洪水总量510.9万m³。沟道内部冲淤强烈，支流入汇处、河宽束窄处上游、桥梁上游、坡度平缓段有大量泥沙淤积，桥梁下游明显冲刷，下游河段岸坡坍塌损毁，河道明显拓宽，多处洪水淹没上岸^[25]。

综合调查结果及相关资料，选取干沟中下游为对象，分析泥沙输移对山洪特性的影响。计算河段自新黄路与宝涞路交汇处至下游流域出口，全长8.4 km。沟道沿线有3处村落及水泥厂、养殖场、风景区等生产服务设施，6条支流依次汇入其中。计算河段纵剖面及洪前河宽如图 2所示，沿程局部比降0.04%~3.2%，平均比降1.04%；岸高为2.5~3.6 m；河宽为8.4~44.4 m，部分河段局部“突扩突缩”，河宽变幅较大。洪水过后，在计算河段上游及中下游选取3处河宽变幅较小、沟道相对顺直的区域，每处取3点样方，量测河床表面及下部0.4 m深处泥沙粒径级配并取均值。测得粒径床面泥沙中值粒径19.8 mm，下部床沙中值粒径约4.2 mm，级配曲线如图 3所示。

2 “水文-水动力学-输沙”模型

山洪通常具有实时量测困难、历史资料少、易灾区域沿沟道狭长分布等特征。为探索泥沙运动对洪水特征的影响，以水文模型描述山洪流域的产汇流过程、提供山洪沟道的入流边界条件，以一维水沙动力学模型描述狭长沟道内的水沙输移与河床冲淤过程。

2.1 水文模块——HEC-HMS模型

HEC-HMS模型被广泛应用于山区小流域的降雨-产汇流过程模拟，在干旱区的城市化区域有山洪应用的成功先例^[26]。根据红螺谷流域的气候、下垫面条件选取合适参数，利用HEC-GeoHMS软件和Arcgis平台对地形数据进行填洼、子流域划分、水系提取，得到模型所需的流域信息；分别用SCS(soil conservation service)单位线法、Muskingum-Cunge法计算产汇流，得到各支流入汇干流的洪水过程线(如图 4所示)。算得干沟入口处洪峰流量149.9 m³/s，洪水总量503.3万m³，符合灾后调查估算结果。应该指出，HEC-HMS模型模拟结果具有一定的不确定性，但本文探索一定洪水过程下的泥沙输移效应，洪水入流过程的不确定性对本文结论没有显著影响。

2.2 水动力学、输沙模块

考虑红螺谷沟洪水季节性强、陡涨陡落特征以及沟道沙石粒径范围宽、床面粗化特征，采用基于浅水方程的山区河流推移质输沙模型。该模型已被用于描述震后山区泥石流阵性补给对河床演变的影响^[27]。

1) 水流控制方程

采用一维Saint Venant方程组描述水流运动。考虑支流垂直入汇，水流方程为：

$ {\frac{{\partial A}}{{\partial t}} + \frac{{\partial Q}}{{\partial x}} = {q_{\rm{L}}}, } $

(1)

$ {\frac{{\partial Q}}{{\partial t}} + \frac{\partial }{{\partial x}}\left( {\frac{{{Q^2}}}{A}} \right) = - gA\frac{{\partial z}}{{\partial x}} - g\frac{{{n^2}{Q^2}}}{{A{h^{4/3}}}}, } $

(2)

其中：t为时间，x为流向坐标，A为过流断面面积，Q为流量，q_L为单位长度上侧向入汇流量，h为水深，g为重力加速度，z为水位，n为综合Manning阻力系数。n表示为

$ n=\gamma n_{\mathrm{c}}+\Delta n. $

(3)

其中：n_c为清水条件下的阻力值，γ为推移质运动阻力增值系数，Δn为跨河建筑导致的阻力增量。模型中，n_c由适于山区河流定床条件的Rickenmann^[28]模型计算：

$ {n_{\rm{c}}} = \frac{{{S^{0.33}}D_{90}^{0.45}}}{{0.56{g^{0.44}}{Q^{0.11}}}}\quad (S > 0.8\% ), $

(4a)

$ {n_{\rm{c}}} = \frac{{{S^{0.08}}D_{90}^{0.24}}}{{2.73{g^{0.49}}{Q^{0.03}}}}\quad (S < 0.8\% ). $

(4b)

其中：S为床面坡降，D₉₀为床面泥沙中累积粒度分布达90%时对应的泥沙粒径。

推移质运动下的水流阻力与定床条件有所不同，引入Song等^[15]的结论反映这一影响：

$ \gamma=\left[30.1\left(\frac{\rho_{\mathrm{s}}-\rho}{\rho} \frac{g}{\nu^{2}}\right)^{1 / 6} D_{50}^{1 / 2} \frac{Q_{\mathrm{bT}}}{Q}+1\right]^{0.51}. $

(5)

其中：ρ_s与ρ分别为沙、水密度，υ为运动黏度，Q_bT为体积输沙率，定床条件下γ取1。

建有跨河桥梁的区段过水断面缩减，水流阻力增加，在郑国栋等^[29]成果基础上，以阻力系数增值Δn概化其增阻壅水作用，推得Δn可表示为

$ \Delta n=\frac{\xi h^{4 / 3}}{4 g \Delta x n_{\mathrm{c}}}. $

(6)

其中：Δx为河段网格长度，ξ为局部损失系数，有桥河段取0.5，无桥河段取0^[29]。本文计算条件下，局部阻力增值的影响幅度、影响范围较为有限。

2) 河床变形方程

假定泥沙运动以推移质为主。将河床分成表面活动层和下部底床层，活动层与推移质直接进行泥沙交换，描述河床形貌调整过程的非均匀沙Exner方程为：

$ \left(1-\lambda_{\mathrm{p}}\right) B \frac{\partial z_{\mathrm{b}}}{\partial t}=-\frac{\partial\left(B q_{\mathrm{bT}}\right)}{\partial x}+q_{\mathrm{sT}}, $

(7a)

$ \begin{aligned}\left(1-\lambda_{\mathrm{p}}\right) B\left[f_{\mathrm{Ii}}\right.&\left.\frac{\partial}{\partial t}\left(z_{\mathrm{b}}-L_{\mathrm{a}}\right)+\frac{\partial}{\partial t}\left(F_{\mathrm{i}} L_{\mathrm{a}}\right)\right]=\\ &-\frac{\partial\left(B q_{\mathrm{bi}}\right)}{\partial x}+q_{\mathrm{si}}. \end{aligned} $

(7b)

其中：q_bi为总单宽输沙率，λ_p为床沙孔隙率，z_b为河床高程，B为河宽，L_a为活动层厚度，q_bi为分组单宽输沙率，i为粒径组编号(本文计算中共分9组)，q_sT、q_bi为单位长度上因边岸拓宽或支流入汇带来的总补沙量和分组补沙量，F_i为第i组级配的泥沙占活动层的体积百分比，f_Ii为活动层与底床层的泥沙交换量。f_Ii表达式为

$ {f_{{\rm{Ii}}}} = \left\{ {\begin{array}{*{20}{l}} {{{\left. {{f_{\rm{i}}}} \right|}_{{z_{\rm{b}}} - {L_{\rm{a}}}}}, }&{\left( {\frac{{\partial \left( {{z_{\rm{b}}} - {L_{\rm{a}}}} \right)}}{{\partial t}} < 0} \right);}\\ {\alpha {F_{\rm{i}}} + (1 - \alpha ){p_{{\rm{bi}}}}, }&{\left( {\frac{{\partial \left( {{z_{\rm{b}}} - {L_{\rm{a}}}} \right)}}{{\partial t}} > 0} \right).} \end{array}} \right. $

(8)

其中：f_i为底床层中i组级配沙的体积分数，p_bi为推移质中i组级配沙的体积分数，α为表征掺混比例的系数，计算中通常假设二者等比例混合，α取0.5^[27]。

使用Wilcock & Crowe^[30]公式计算推移质的单宽体积输沙率：

$ q_{\mathrm{bi}}=\frac{F_{\mathrm{i}} u_{*}^{3}}{g R} W_{\mathrm{i}}^{*}. $

(9)

其中：u_*为摩阻流速，R=(ρ_s-ρ)/ρ，W_i^*为无量纲推移输沙强度。

$ W_{\rm{i}}^* = \left\{ {\begin{array}{*{20}{l}} {0.002{\varphi ^{7.5}}, }&{(\varphi <1.35);}\\ {14{{\left( {1 - \frac{{0.894}}{{{\varphi ^{0.5}}}}} \right)}^{4.5}}, }&{(\varphi > 1.35).} \end{array}} \right. $

(10)

其中：φ=τ_b/τ_ri为水流强度参数，τ_b为水流切应力，τ_ri为第i组级配推移质泥沙的参照剪切力，输沙所用的肤面阻力仅为式(3)总阻力的一部分，采用Manning-Strickler公式^[31]计算。

灾后考察得知，山洪过程中沟道中下游多呈现边岸侵蚀坍塌的侧向泥沙补给现象^[12]，选取工程界常用的Osman & Thorne模型^[32]，计算单位时间内边岸侵蚀补沙量：

$ q_{\mathrm{sT}}=\frac{C_{1}\left(\tau_{\mathrm{f}}-\tau_{\mathrm{c}}\right)}{\gamma_{\mathrm{bk}}} e^{\left(-1.3 \pi_{\mathrm{c}}\right)}\left(z_{\mathrm{d}}-z_{\mathrm{b}}\right) \Delta t. $

(11)

其中：τ_f=ρgn²u³h^-1/3为水流作用在边岸上的切应力，τ_c=0.047(ρ_s-ρ)D₅₀为边岸抵抗冲刷的临界剪切力，γ_bk为边岸容重，z_d为河岸高程，C_l为边岸侵蚀系数，计算中取常用推荐值0.09^[33]。

2.3 计算条件

以红螺谷沟“7.21”山洪过程为典型案例。应用上述模型描述沟道形态随底部冲淤、边岸侵蚀动态变化。为揭示泥沙输移与河床冲淤对最高水位、最大流速、淹没时间等山洪特征值空间变化的影响，也采用定床模型(沟道形态固定、无泥沙输移)结果进行对比。模型计算条件见表 2，考虑到流域产汇流、洪水波演进所需时间，计算时段取30 h。2种模型均以HEC-HMS模型在“7.21”暴雨条件下的计算结果作为入流过程，以均匀流的水位流量关系作为出口条件，阻力系数取出口网格阻力n值。灾后调查发现，流域沟道内泥沙补给充足，干支流交汇处淤沙较多，模型入口、支流来沙量以式(9)计算。针对局部损失系数ξ、掺混比例系数α、边岸侵蚀系数C_l等关键参数的敏感性分析发现，各参数在推荐取值范围内的变化对计算结果影响较小。

3 泥沙输移影响下洪水要素的时空分布特征

山洪演进过程中，沿程断面的水位、流速等山洪要素决定了灾害的影响范围与程度。山洪沟道的河床冲淤对洪水位和洪水流速有直接影响。

3.1 沿程泥沙冲淤输移

山洪含沙量的时空变化情况如图 5所示。洪水演进过程中，断面含沙量受河宽、坡降、桥梁等因素影响沿程分布不均，并随水流、地形、级配条件的变化不断调整。同一时刻，含沙量沿流向递减区域发生落淤(如3 000 m附近)，递增河段发生冲刷(如1 000 m附近)，洪峰过境、回落时段含沙量较高。图 6给出了断面河底冲淤变化值在洪前、5 h、10 h、20 h、30 h等典型时刻的沿程分布，从而反映出山洪过程中沟道纵剖面的调整变化。可以看出，多数断面冲刷、淤积程度随时间推移持续加深，部分断面经历了“先冲后淤”(如1 800、3 300、5 900 m附近)或“先淤后冲”(如1 600、4 300 m附近)的过程。15 h之后，各断面河底高程基本不再变化。对照图 2中原始地形可得，冲淤变化与沟道形态关系密切。坡降由缓变陡处多出现冲刷，由陡渐缓处多出现淤积，河宽束窄段、跨河桥梁处呈现“上淤下冲”的特征。计算结果表明：30 h后8.4 km的干沟中有5.5 km长发生淤积，淤积总量8.17万m³，断面最高淤抬1.55 m；2.9 km长处于冲刷状态，侵蚀总量3.06万m³，断面最深冲刷1.88 m。该结果与灾后调查所得房山区山洪沟冲淤统计规律相符，说明模型在描述山洪过程中泥沙输移时具有一定的可靠性。

3.2 沿程最大表观水深

山洪演进过程中，洪水漫滩淹没是造成财产损失、人员伤亡的重要原因。定义沟道各断面处表观水深Z为水面高程与原河底高程之差，则最大表观水深Z_max是确定洪水淹没范围、评估灾害损失的重要参数。灾后调查发现，接近下游沟道出口的娄子水小学断面岸上洪痕超过原河底高程4.5 m。动床、定床条件下，算得该断面最大表观水深随时间的变化过程如图 7所示。动床模型所得Z_max为4.57 m，与洪痕调查结果基本相符。

图 8给出了计算所得的最大表观水深沿程分布及其与定床结果的对比。因底坡、河宽变化及桥梁的增阻壅水作用，河道束窄段水位普遍较高。

结合图 6可知，与定床计算结果比较，淤积段Z_max多有抬升，冲刷段Z_max普遍下降。选取间隔50 m的代表性断面，受到断面淤积和输沙增阻的影响，74.6%的断面数的Z_max较定床值抬升，32.5%的断面数的Z_max抬升幅度超过定床值的30%，增幅最大处较定床值增加103%，全河段Z_max较定床条件平均抬升14.1%，与王海周等^[34]通过试验验证的水沙变化响应关系一致。由此可见，输沙耗能增阻、河底边岸冲淤、表层级配变化致使泥沙输移影响下沿程最大表观水深与定床条件出现明显差异，坡降由陡渐缓段、河宽束窄段、跨河桥梁处壅水淹没风险加大。

3.3 沿程最大流速

同等水位下，流速是决定洪水危险度的另一关键参数，流速越高，洪水冲击掏蚀、席卷挟带的能力越强，对沿线财产、人员的威胁越大。当淹没水深大于0.3 m时，3 m/s的洪流便可构成严重威胁^[35-36]。图 9给出了计算得到的“7.21”洪水过程中沟道沿线最高流速U_max。

可以看出，各断面U_max普遍高于3 m/s，且因地形变化和桥梁的存在沿程上下波动。相对于定床条件，动床U_max沿程波动幅度较小，体现了河床的自动调整作用。受输沙耗能增阻作用的影响，动床U_max普遍低于定床假设下的计算值，这与大部分断面的最大表观水深增加一致；其中，72.8%的断面动床U_max低于定床值，最多减小25.1%，平均降低6.3%。

3.4 演进过程特征

定床假设下，通常认为洪峰自上游至下游传播演进，断面流速、淹没水位均在流量最大时刻T_Qmax(峰现时间)前后达到最大值，单峰洪水流经简单一维河道时，下游断面达到最大流量、最高水位、最高流速的时间一般应晚于上游断面。故在工程应用中，通常以少数监测点位观察到的洪水涨落现象作为判断流域整体洪水变化趋势的依据。

图 10a和10b分别给出了各断面达到流量最大值、水位最大值的时刻T_Qmax、T_Zmax沿程分布及其与定床结果的对比。为了表征洪水特征值维持在高数值的持续时间，将流量、最大表观水深高于90%最大值的时间区间P_Q、P_Z也绘制在相应图中。

可以看到在定床条件下，洪峰于第479 min进入沟道，历经21 min传播至出口断面；各断面峰现时间T_Qmax总体上随流向距离x增加而滞后；由于存在沿程支流入汇情形，在支流入汇处T_Qmax存在小幅波动。动床条件下，T_Qmax沿程有类似的波动。2种情形下各断面实时流量大于断面最大流量90%的时间区间P_Q基本一致，洪峰基本上仍符合自上游向下游传播的特征。

对于各断面达最大表观水深的时刻T_Zmax而言，动床与定床情形的结果呈现较大差异。定床条件下，T_Zmax沿流向均匀推后，与洪峰传播时间T_Qmax的沿程变化规律基本一致。各断面水位高于90%Z_max的持续时长P_Z差异较小。而在实际情况下，考虑泥沙运动后的断面T_Zmax、P_Z受冲淤的影响均发生了显著的变化，与T_Qmax、P_Q的变化不再保持同步。T_Zmax相较于定床条件最早提前121 min，最迟滞后143 min，P_Z值在30~430 min之间大幅波动，部分河段(如x=500 m处，见图 11)水位雍高的时间显著增加。

图 12给出了定、动床条件下典型河段(1 500~2 500 m)河底高程、水深、淹没水位随时间的变化过程。由图可知，定床条件下各断面达到最高水位的时刻基本一致。动床条件下，1 800、2 400 m断面附近河段侵蚀下切，水位较早开始回落，T_Zmax较上游部分断面有所提前；而在1 500、2 200 m断面附近，河床淤积抬升，水位持续提高，T_Zmax较下游一定距离内的断面有所滞后。

结合沿程冲淤情况(见图 6)可知，实际动床情形下T_Zmax沿程分布与断面冲淤关系密切，淤积区段达到最高水位的时刻相对滞后，维持在高水位的时段有所延长；冲刷区段各断面达到最高水位的时刻有所提前，维持在高水位的时段相对缩减。单峰洪水演进过程中，有可能因剧烈冲淤导致“下游A断面水位已在回落，而其上游B断面水位仍在抬升”的现象。因此，山洪预警管理中，淤积段淹没时间可能较定床结果显著延长，单一站点观测到的水位回落现象无法作为上游河段水位均已开始下降的判据，不能用有限个点位的水位变化结果判断全程各处的变化趋势。布设监测点位时，应充分考虑泥沙冲淤输移可能带来的影响。

4 结论

本文采用“水文-水动力学-输沙”模型，基于野外调查数据，复演了北京红螺谷沟“7.21”山洪过程，探索了泥沙输移对山洪特征值及其时空分布的影响。泥沙输移显著改变沟床形态、泥沙级配与水流阻力。沟道纵坡连续减小段多出现淤积，连续增加段多出现冲刷；河宽束窄段、跨河桥梁处多出现“上淤下冲”现象。泥沙输移显著影响最高洪水位、最大流速等洪水特征值的空间分布，大幅提升局部河段山洪危险。在阻力变化与河床冲淤的共同作用下，各断面最高洪水位较定床情况整体抬升，局部有所下降。动床条件下，断面流速最大值稍低于定床值。泥沙输移显著影响洪水过程，水位变化与流量变化过程不再同步。定床假设下各断面达到最高水位的时间T_Zmax与峰现时间T_Qmax基本一致。而实际动床条件下，T_Zmax与T_Qmax最多存在±2 h的时间差，淤积受淹段高水位淹没时间显著延长。当前防灾实践中常用的“以少数点位监测结果判断全程各处洪水变化趋势”的预警方法可靠性有限，这一现象的认识对于山洪预警工作具有一定参考价值。受实测山洪水沙数据较少的限制，本文结论的定量准确性仍可提升完善。后续将结合现场观测，进一步深入研究。