自Leopold^[1]将冲积河流分为顺直、弯曲和辫状河流以来，国内外学者对河型分类开展了广泛的研究，国内具有代表性的应用方案是把冲积河流分为直流河、曲流河、分汊河、网状河及辫状河^[2-3]。

泥沙过多是辫状河道形成的重要原因^[4]，尤其是河流变化的推移质输沙量影响辫状河型^[5]的形成。汛期大流量条件下，上游来沙量远大于下游辫状河流水流挟沙力，河床淤积抬升，使其形成不稳定性沙洲，呈现出辫状交织的平面形态。如：雅鲁藏布江地形呈现出宽窄相间的特点，当大量挟带泥沙的冰雪融水从上游峡谷段进入辫状河流时，河床平缓，流速大幅减小，河水在河床漫滩，河床冲淤变形剧烈，挟带泥沙的水流在辫状流动的过程中大量泥沙沉降，辫状河流对泥沙产生滞留作用，上游来沙在辫状河流发育河段的输运发生变异，不能正常地随水流输运到下游。

Griffiths等^[6]认为辫状砾石河流的推移质输运能力随着下游水力学条件的变化而发生改变，受洪水期间推移质输运的影响，推移质输运能力的差异会导致河床物质的改变；反过来，河床物质的改变会使床面产生不规律的波动及泥沙输运波动，这些波动是辫状河流中推移质输运的作用机制。在辫状河流中，泥沙输运尤其是推移质输运对地形的动态变化至关重要^[4]，河床发生泥沙运输和地形变化的区域被定义为活动宽度，是辫状河形态动力学研究非常有用的要素，可通过辫状强度、浸湿宽度等参数求得活动宽度，并用于预测推移质输沙量^[7]。Yang等^[8]采用水动力和泥沙输运数值模型，模拟了在理想河流中悬移质控制的辫状河流的演变过程，模拟了天然辫状河中常见的河道节点和演变过程，如心滩；采用多级泥沙组分的方法模拟了泥沙粗化和细化对新塘和辫状河道演化的影响。

辫状河流发育河段宽窄、多流路、流速小的特征使泥沙输运发生变异，从而会使生源要素向下游的输运受到影响，从而影响下游的水生态系统。在泥沙输运方面，河流泥沙输运研究已经有了很多成果，但是基于辫状河流复杂的泥沙迁移研究甚少，而针对辫状河流及其输沙的研究基本致力于泥沙输运对辫状河流地形形成与变化的影响机制，而对于辫状河流的存在对泥沙输运所产生的冲淤影响研究不足。基于上述出发点，拟通过数值模拟对辫状河流输沙变异开展研究。主要研究辫状河流在不同来水、来沙条件下，挟沙水流流经辫状河流后的流量过程、含沙量过程和悬移质输沙率过程及其进出口过程的相位发生的变异规律。

1 数学模拟方法 1.1 模型简介

研究构建了可用于辫状河流输沙模拟的数学模型，该数学模型通过全沙输移理论计算悬移质和推移质输沙率。在计算泥沙输移的同时，考虑泥沙的启动、沉降淤积作用导致的地形更新，进而求解挟沙水流在辫状河流中运动时对地形变化的影响。

1.2 水动力控制方程

泥沙输运、冲淤的求解计算是建立在水动力学模型的基础之上，即求解连续性方程与雷诺时均N–S方程。其中，水动力学模型根据静水压力假定进行简化。假定流速沿水深保持不变，将不可压缩、取雷诺平均的Navier–Stokes控制方程沿水深进行积分，得到沿水深平均的2维浅水流动质量和动量守恒控制方程组：

1）连续性方程：

$ \frac{\partial h}{\partial t}+\frac{\partial h \bar{u}}{\partial x}+\frac{\partial h \bar{v}}{\partial y}=0 $

(1)

2）动量方程组：

$ \begin{aligned}[b] \frac{{\partial h\overline u }}{{\partial t}} +& \frac{{\partial h{{\overline u }^2}}}{{\partial x}} + \frac{{\partial h\overline {uv} }}{{\partial y}} = - gh\frac{{\partial \eta }}{{\partial x}}- \\ &\frac{h}{{{\rho _0}}}\frac{{\partial {p_{\rm{a}}}}}{{\partial x}} - \frac{{g{h^2}}}{{2{\rho _0}}}\frac{{\partial \rho }}{{\partial x}} - \frac{{{\tau _{{\rm{b}}{x}}}}}{{{\rho _0}}} + \\ &\frac{\partial }{{\partial x}}\left( {2h{\nu _{\rm{t}}}\frac{{\partial \overline u }}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {h{\nu _{\rm{t}}}\Bigg(\frac{{\partial \overline u }}{{\partial y}} + \frac{{\partial \overline v }}{{\partial x}}\Bigg)} \right) \end{aligned} $

(2)

$ \begin{aligned}[b] \frac{{\partial h\overline v }}{{\partial t}} + & \frac{{\partial h{{\overline v }^2}}}{{\partial y}} + \frac{{\partial h\overline {uv} }}{{\partial x}} = - gh\frac{{\partial \eta }}{{\partial y}}- \\ &\frac{h}{{{\rho _0}}}\frac{{\partial {p_{\rm{a}}}}}{{\partial y}} - \frac{{g{h^2}}}{{2{\rho _0}}}\frac{{\partial \rho }}{{\partial y}} - \frac{{{\tau _{{\rm{b}}y}}}}{{{\rho _0}}} + \\ &\frac{\partial }{{\partial y}}\left( {2h{\nu _{\rm{t}}}\frac{{\partial \overline v }}{{\partial y}}} \right) + \frac{\partial }{{\partial x}}\left( {h{\nu _{\rm{t}}}\Bigg(\frac{{\partial \overline u }}{{\partial y}} + \frac{{\partial \overline v }}{{\partial x}}\Bigg)} \right) \end{aligned} $

(3)

式（1）～（3）中：t为时间，s；η为水面相对于未扰动基面的高度，m；h为总水深，m； $\overline u $ 、 $\overline v $ 为基于水深的平均流速，m/s； $u$ 、 $v$ 为空间流速沿x、y方向的速度分量，m/s； $\;\rho $ 为水的密度，1 000 g/L； $\;\rho_0 $ 为水的相对密度，g/L；p_a为当地大气压强，Pa；g为重力加速度，9.81 m/s²； $\nu_{\rm{t}} $ 为水平方向涡粘系数，m²/s； $\tau_{{\rm{b}}x} $ 、 $\tau_{{\rm{b}}y} $ 为水面与河床的水流切应力在x、y方向的分量，N/m²。

1.3 泥沙输运模型及方程

非黏性泥沙模型的数值计算基于水动力条件进行，主要用于计算单一水流或波流共同作用下的非黏性泥沙输移能力，以及引入河床变形速率的反馈机制进行地貌演化计算，本文计算单一水流作用下的非黏性泥沙输运。

由对流和湍流扩散控制悬浮泥沙输运的方程^[9]为：

$ \begin{aligned}[b] \frac{{\partial h\overline C }}{{\partial t}} \!+\! \frac{{\partial h\overline u \overline C }}{{\partial x}} \!+\! \frac{{\partial h\overline v \overline C }}{{\partial y}} = h\bigg(\frac{\partial }{{\partial x}}\bigg(\frac{{{\nu _{\rm{t}}}}}{{{\sigma _{\rm{t}}}}}\frac{\partial }{{\partial x}}\bigg)+ \frac{\partial }{{\partial y}}\bigg(\frac{{{\nu _{\rm{t}}}}}{{{\sigma _{\rm{t}}}}}\frac{\partial }{{\partial y}}\bigg)\bigg)\overline C + \Delta S \\ \end{aligned} $

(4)

式中： ${\sigma _{\rm{t}}}$ 为普朗特数； $\overline C$ 为水深平均的悬移质泥沙浓度，g/m³； $\Delta S$ 为沉积汇， $\Delta S=\phi_{0}\left(\eta_{0}\right) \omega_{{\rm{s}}}\left(\overline{C}-\overline{C}_{{\rm{e}}}\right)$ ， $\phi_{0}$ 为单宽输沙率函数， $\eta_{0}$ 为床面以上标准高程， $\overline{C}_{{\rm{e}}}$ 为平衡浓度，g/m³， ${\omega _{\rm{s}}}$ 为悬移质泥沙的沉速。

Engelund–Hansen泥沙输移计算公式^[10]：

$ \phi = \frac{S}{{\sqrt {(s - 1)g{d^3}} }} $

(5)

式中：S为输沙率（推移质、悬移质和全沙），m²/s；d为泥沙粒径，m； $\phi $ 为无量纲输沙率。输沙率计算公式如下：

$ {S\!_{{\rm{tl}}}} = 0.05\frac{{{C^2}}}{g}{\theta ^{2.5}}\sqrt {(s - 1)gd_{50}^3} $

(6)

$ {S\!_{{\rm{bl}}}} = {k_{\rm{b}}}{S\!_{{\rm{tl}}}} $

(7)

$ S\!_{{\rm{sl}}}=k_{{\rm{s}}} S\!_{{\rm{tl}}} $

(8)

$ \theta=\frac{\tau}{({\rm{s}}-1) \rho g d_{50}} $

(9)

式中：d₅₀为泥沙中值粒径，m； $\tau $ 为水流剪切力；S_bl为推移质输沙率，m²/s；S_sl为悬移质输沙率，m²/s；S_tl为全沙输沙率，m²/s；k_b为推移质输沙率校正系数；k_s为悬移质输沙率校正系数；s为泥沙相对密度，g/L。

由于泥沙输移产生的河床变形主要取决于河床高程变化速率，其大小可通过Exner方程（泥沙输移连续性方程）求解：

$ -(1-n) \frac{\partial {\textit{z}}}{\partial t}=\frac{\partial S\!_{x}}{\partial x}+\frac{\partial S\!_{y}}{\partial y}-\Delta S $

(10)

$ {{\textit{Z}}_{{\rm{new}}}} = {{\textit{Z}}_{{\rm{old}}}} + \frac{1}{{1 - n}}\frac{{\partial {\textit{Z}}}}{{\partial t}}\Delta {t_{{\rm{HD}}}} $

(11)

式中，n为床沙孔隙率，z为河床高程，S_x、S_y为全沙在x、y方向上的输移量， $\Delta $ t_HD为水动力学模块中的时间步长，Z为高程，Z_new、Z_old分别为新河床高程和旧河床高程。

1.4 模型参数设置

模型中水动力模块主要有水平涡黏系数和糙率，非黏性泥沙输运模块主要参数有床沙孔隙率、床沙相对密度和糙率，非黏性泥沙输运模块中的糙率取值与水动力模块一致。各参数值具体设置如表1所示。

2 工况设置及地形模型 2.1 工况设置

在雅鲁藏布江中游河段采取了28个水样样品，并做级配分析，如图1所示。图1中，P为小于某颗粒的累计质量百分比，Y–6–L(R)、Y–7–L(R)、Y–8–L(R)、Y–9–L(R)为水样采样点，分布在曲水至桑日河段。并通过对雅鲁藏布江中辫状河流典型河段上下游水文站流量、含沙量和悬移质输沙率的年内过程的分析，得到两个站点各自的流量、含沙量和悬移质输沙率年内变化过程线之间出现不同程度的相位差，尤其是上下游含沙量过程之间产生较大的相位差，下游含沙量峰值滞后于上游含沙量峰值至少2 d。

泥沙过多时辫状河流形成的原因也是结果，在地形条件复杂多变的辫状河流中，势必对泥沙输运过程造成影响。根据对雅鲁藏布江各段辫状河流的形态特征分析，辫状河流多出现在山区河流的宽谷段，挟沙水流从峡谷段进入宽谷段的辫状河流后，流速迅速减小，泥沙难以被水流带走而出现泥沙滞留和冲刷现象。本文针对辫状河流中的水沙运动，采用数值模拟结合室内水槽铺沙试验的方法对其进行了研究。

设置6个工况，研究模拟来流洪峰流量10 L/s时，中值粒径分别为0.3、0.2、0.1、0.05、0.02和0.01 mm的泥沙在辫状河流中的输运变异规律。

2.2 辫状河流地形模型

辫状河流地形形态如图2所示，地形由滩地和辫状水流通道组成，计算范围宽3 m，长17 m，左边为进口，右边为出口。为模拟天然辫状河流在宽窄相间河道的宽谷段发育的特征，进口宽度设置为0.3 m，位于进口断面正中央，滩地和辫状水流通道的地形坡度均为1.07%。辫状河流水具有心滩、分汊和流线3个方面的辫状强度特征，共有5个指标来定量描述：心滩辫状指标BI_B^[11]和BI_B*^[12]分别为3.1和4.3，分汊辫状强度指标BI_T^[13]为3.6，流线弯曲度辫状强度指标P_T^[13]和P_T*^[14]分别为4.7和4.1。

对地形变化较大或者水流进口区域进行了加密，均采用非结构三角形网格，网格单元的最大面积为0.008 m²，加密区域的为0.002 5 m²，网格数目为11 778个。为保持模型计算的稳定性，计算时间步长设置为5 s。初始流场为静止状态。进口条件为流量、泥沙浓度边界，且在断面上均匀分布；出口条件为水位边界，其他岸坡设为陆地边界，即法向速度为0。

3 辫状河流输沙特征及变异分析 3.1 泥沙输运特征

辫状河流上游来流洪峰过程时长为9 h，峰值分别为10.0和1.5 g/L。当来流泥沙粒径为0.2 mm时，洪峰过后有49.2 kg泥沙在辫状河流中淤积，洪峰过程时间段记为0：00至9：00，4：00达到洪峰值。选取河床地形高程有显著变化的6个时刻（0：06、3：20、3：40、4：00、4：40、9：00）分析河段入口及下游辫状河流的冲淤变化，图3为计算河段入口及下游辫状河流冲淤变化情况，辫状河流地形在洪峰来临前后产生的冲淤变化如图4所示。由图4可知，洪峰来临前2.5 h，入口和下游的最大冲刷厚度分别是2.0 cm和0.7 cm，下游最大淤积厚度0.96 cm。3：20至4：00之间，来流流量和含沙量迅速增大到洪峰值，辫状河流河床高程随之迅速增大，最大淤积厚度达到1.5 cm。4：00洪峰过后，来流流量和含沙量迅速减小，冲刷和淤积厚度仍在增加，入口冲刷厚度直到5：20达到最大值4.3 cm，然后开始减小，直到9：00计算结束，冲刷厚度减小至3.8 cm，表明入口处不再冲刷，而形成淤积，最大淤积厚度为0.5 cm。4：40时，下游的辫状河流冲刷厚度达到最大值1.9cm，之后一直减小；9：00时减小至1.35 cm，表明随着洪峰过后，来流流量迅速减小，冲刷区域不再继续冲刷，而开始淤积，最大淤积厚度为0.55 cm。下游辫状河流淤积厚度在5：30时达到最大值1.8 cm，之后趋于稳定，表明随着来流量和含沙量减小，辫状河流中淤积趋于稳定。

泥沙输运是水流和泥沙相互作用的过程，河床高程的变化导致水深和流场产生相应变化，如图5和6所示。随着来流流量增大，直至达到洪峰流量，流动区越来越大，流速和水深也越来越大，4：00洪峰来临，入口流速和下游辫状河流中流速及水深均达到最大值，入口最大流速0.61 m/s，下游最大流速0.34 m/s，最大水深4.4 cm，最大流速和最大水深均出现在河床高程变化最大区域，即冲刷最严重区域；洪峰过后，入口和下游流道中流速越来越小，但是入口处由于持续冲刷，水深仍在增大，直到4：40达到最大值5.8 cm。

与0：06时相比，9：00时与之有着相同的流量，但流场和水深却产生了很大变化，0：06时淤积还未对流场和水深产生影响，入口流速大，辫状分汊河流中流速和水深均匀，流动区面积大。9：00时，入口处由于冲刷导致水深增大3.7 cm，流速减小0.12 m/s，下游辫状分汊河段中最大流速增大0.08 m/s，最大水深增大1.9 cm，但只有在局部冲刷区域出现流速增大，其他区域流速减小，且流场不再均匀，流动区的面积大幅度减小；虽然在局部冲刷区域水深增大，但冲刷区域很小，且在其他区域水深减小或不变，静水区面积变大，表明淤积使辫状分汊河段中水流在辫状河流中漫滩。

3.2 输沙变异分析

各粒径泥沙出口流量过程线见图7。由图7可知，水流流经辫状河流后，出口流量过程与进口流量在时间上完全同步，没有产生相位差。但是由于泥沙孔隙率的存在，使其能够储存一部分水，故在没有分流的情况下，河段出口流量相对于进口流量有所减小。泥沙粒径为0.01、0.02 mm的挟沙水流出口流量相对进口只减少了3%，挟带大颗粒泥沙水流出口流量减小20%以上。

挟沙水流经过辫状河流后流量减小，含沙量过程也产生了进出口过程不同步和含沙量大小的变化。相位差定义为出口含沙量和悬移质输沙率峰值滞后或者提前于进口含沙量峰值的时间与总时间历时之比，若提前，则相位差为负值；若滞后，则相位差为正值。

总体上，出口含沙量过程沙峰滞后于入流含沙量过程，随着泥沙粒径减小，滞后程度随之减小。泥沙中值粒径为0.02、0.01 mm的水流出口含沙量过程基本与入流含沙量同相位，相对于来流含沙量沙峰时间点4：00，分别提前6 min和滞后8 min，如表2所示，其相位差分别为–1%和1%。中值粒径为0.02、0.01 mm的挟沙水流进入辫状河流时，出口含沙量在2 min之内增大到较大值，中值粒径为0.01 mm水流甚至增大到4.153 g/L，大于4：00沙峰和洪峰流量时的含沙量值3.907 g/L。中值粒径为0.3、0.2、0.1 mm水流出口含沙量峰值滞后于4：00，均超过了1 h，相位差超过了15%，粒径越大，相位差越大。泥沙中值粒径大于0.05 mm的水流出流含沙量过程出现了双峰或者多峰，0.1和0.2 mm泥沙则在第一个峰值之后波动较大，但整体上变化规律仍与来流含沙量过程相对应。

辫状河流挟沙水流出流含沙量过程与入流含沙量过程相比，除产生相位差以外，其含沙量大小也产生了相应的变化。由于水流挟带细颗粒泥沙的能力比挟带粗颗粒泥沙的能力大的缘故，不同中值粒径泥沙的冲淤程度有所不同，如表3所示。

由表3可知，9 h内，当入流来沙量均为57.0 kg时，0.30、0.20、0.10 mm泥沙产生淤积，0.05、0.02、0.01 mm泥沙产生冲刷，随着泥沙粒径增大，辫状河流出流输沙量随之减小，表明泥沙粒径越大，越容易淤积，但是0.01 mm泥沙出流输沙量却小于0.02 mm泥沙。水流出流含沙量值的大小随所挟带泥沙粒径的变化与冲淤量的变化一致，中值粒径为0.30、0.20和0.10 mm的泥沙出流含沙量峰值相对于入流含沙量减小，但减小程度随着粒径减小依次变弱，表明这3种粒径泥沙整体上沿程淤积，辫状河流中泥沙滞留，随着泥沙中值粒径减小，淤积量变小。中值粒径为0.05、0.02以及0.01 mm的泥沙出流含沙量相对于入流含沙量增大，表明这3种粒径的泥沙整体上沿程产生了冲刷，辫状河流中泥沙释放，随着粒径减小，释放程度增大，但是中值颗粒0.01 mm泥沙释放程度小于0.02 mm泥沙。

各工况入流加沙量峰值均为0.005 kg/s，入流加沙量峰值点均为4：00，辫状河流出流悬移质输沙率与入流来沙量相比产生的变化，与出流含沙量过程相对于来流含沙量过程所产生的变化相似。如表4和图8所示，悬移质输沙率过程同样产生了相位差，但是比含沙量过程产生的相位差稍小，出口悬移质输沙率峰值滞后于入流来沙量，随着泥沙粒径减小，相位差随之减小，泥沙中值粒径为0.02 、0.01 mm的水流出口悬移质输沙率过程均滞后2 min，其相位差很小，可认为是同相位。但是中值粒径为0.30、0.20 mm泥沙出流悬移质输沙率峰值滞后均超过了1 h，相位差超过了18%。图8（a）～（c）中中值粒径大于0.05 mm的水流出流含沙量过程出现了双峰且峰值大小出现增大、减少、相当3种情况，但整体上变化规律仍与来流含沙量过程相对应，这是由于不同程度的淤积产生的。

3.3 试验模型验证

通过室内水槽铺沙试验，对比流量分析出口含沙量相关性来验证模型计算的合理性。模型计算采用与试验相同的边界条件，进口采用相同的流量边界，出口采用相同的水位边界，试验中记录了水槽出口位置流量过程。

在水槽试验中，由于水槽铺床沙的多孔性，能够存储一部分的水量，因此，流出水槽中的水流流量略微小于进入水槽中的流量。模型计算中水动力模块的计算出流量是只考虑水动力的计算结果，其大小与进口流量相同，而加上泥沙输运模块之后的计算出流量比进口流量小，这与试验中出现的进口流量大于出口流量的现象一致，如图9所示。

模型计算泥沙输运模块出口流量与试验出口流量能较好的吻合，绝对相对误差范围为0.6%～11.7%。

各中值粒径泥沙出口悬移质含沙量试验值与模拟值相关性如图10所示。由图10中各中值粒径泥沙出口悬移质含沙量试验值与模拟值对比结果可知，数值模拟结果可信。

4 结论与展望

作者研究了不同来水来沙条件下辫状河流输沙变异度，即出口含沙量和悬移质输沙率相对进口的变化，出口含沙量和悬移质输沙率峰值相对进口的相位差，以及出口输沙量变化。

辫状河流本身的特征使其输沙产生上述变异，随着泥沙在辫状河流中淤积，静水区增大，流动区减小，表明淤积使得水深变浅，水流在辫状河流中漫滩。辫状分汊河段中局部冲刷区域水深增大，流速增大，其他区域水深和流速减小。

当模拟来流洪峰流量10 L/s时，辫状河流出口含沙量过程和悬移质输沙率过程峰值滞后于入流过程，形成相位差。泥沙粒径越小，相位差越小，中值粒径0.01和0.02 mm泥沙相位差约为0，出口悬移质输沙率产生的相位差小于含沙量；中值粒径0.30、0.20和0.10 mm泥沙在辫状河流中的淤积量超过50%。

综上，辫状河流输沙变异度可用出口含沙量变化百分比、出口悬移质输沙率变化百分比和两者的相位差来描述，这对类似的数值研究和物理试验具有一定的借鉴意义。

辫状河流水沙输运模型的建立、验证和实际应用是一项系统而又复杂的工作，诸多研究都未能完满地解决这一问题。因为辫状河流中的泥沙运动牵涉因素繁多，加上时间限制，该研究尚存在许多有待完善和深入的地方，现整理如下，为后续研究工作的深入开展提供借鉴。

1）水槽试验中，本文采用中值粒径较大的泥沙进行试验，得到较好的试验结果，但是更小颗粒泥沙在辫状河流中的输运仍需探究；

2）模型中只能计算单一粒径泥沙，与天然辫状包含各级粒径泥沙的混合沙有所差别，下一步计算中可考虑混合沙的计算。