李林林,夏军强,邓珊珊,周美蓉,李志威

(武汉大学水资源与水电工程科学国家重点实验室,湖北 武汉 430072)

三峡工程运行后,大量泥沙在坝前淤积,下泄沙量大幅减少,使得水流经常处于严重次饱和状态,坝下游河段往往会发生长距离长历时的河床冲刷[1-4]。河床剧烈冲刷变形会对河道防洪、航运安全、两岸居民生活及工农业用水等带来一系列不利影响。当前对坝下游河床长距离长时段冲刷变形的预测精度还有待提高,对冲积河流不平衡输沙规律的认识也需要进一步深入,研究三峡大坝下游典型河段的泥沙不平衡输移问题具有重要的理论意义及实际应用价值。

在计算悬移质不平衡输沙时,恢复饱和系数是泥沙数学模型中的重要参数,其取值的合理性直接影响数学模型预测结果的准确性。国内外学者针对泥沙恢复饱和系数开展了大量的实验分析及理论研究,但一般都是基于泥沙连续方程和河床变形方程[5-8]。基于泥沙连续方程和河床变形方程,窦国仁[5]将恢复饱和系数解释为泥沙沉降概率,其计算结果恒小于1；韩其为[9]假定不平衡输沙和平衡输沙的河底含沙量梯度相同,对二维扩散方程积分后得出恢复饱和系数为河床近底含沙量与垂线平均含沙量的比值,但其计算值大于1,并建议在淤积时取0.25,冲刷时取1；Zhou和Lin[10]试图沿横向积分以降低其数值,但是由于只考虑了流速分布的影响,并未反映扩散及“恢复”的作用。为了使数模计算结果与实际情况更加接近,部分学者对非均匀沙恢复饱和系数进行研究,从不同角度提出了恢复饱和系数的经验或半经验半理论公式[6,11-15]。王新宏等[7]采用概率统计方法指出分组沙的恢复饱和系数与该粒径组泥沙的沉速成反比,提出了半经验半理论的分组沙恢复饱和系数计算公式；韦直林等[12]认为恢复饱和系数只能通过模型验证来率定,提出了分组沙恢复饱和系数与其沉速成反比的经验关系；葛华等[8]基于三峡水库蓄水后荆江河段实测水沙资料,用蓄水前多年平均输沙水平代替不平衡输沙方程中的挟沙力,计算了该河段的恢复饱和系数,并通过分析指出恢复饱和系数通常随着泥沙粒径(沉速)的增大而减小,随着河床冲刷历时的增加和床沙粗化程度提高而呈递减趋势。另外,Garcia和Parker[16]也提出了恢复饱和系数的经验计算公式。但是,上述经验公式结构较为复杂,应用范围有限,不具有普适性。因此,韩其为和陈绪坚[6]根据泥沙运动统计理论,建立了不平衡输沙的边界条件方程,得到恢复饱和系数的理论表达式，在此基础上,从基于悬移质扩散理论的河床变形方程出发,提出了底部恢复饱和系数的概念；但在应用该公式计算恢复饱和系数时,需先确定与河床冲淤状态有关的非饱和调整系数,主观性较强。

综上所述,虽然对恢复饱和系数进行了大量研究,但对于坝下游河道恢复饱和系数的确定、应用范围等方面还缺乏理论共识，有必要从理论公式推导出发,研究坝下游河段恢复饱和系数的沿程变化特点,尤其是分组悬沙恢复饱和系数的取值,从而提高数学模型的计算精度和可靠性。因此,本文基于Markov随机过程及泥沙起动理论,提出非均匀悬沙恢复饱和系数计算式,可以同时考虑床沙交换过程中推移质与床沙组成的影响。

1 研究河段概况

1.1 荆江河段概况

如图1所示,荆江河段上起枝城,下迄城陵矶,全程长约347 km,其中枝城距离三峡大坝约102 km；以藕池口为界分为上、下荆江,河长分别约为172 km、175 km。上荆江河段由枝城、江口、沙市等6个河湾及过渡段组成,且以杨家垴和观音寺为界又分为枝江、沙市及公安3个子河段,属于微弯分汊型河道；其中沙市距离三峡大坝约192 km。下荆江河段由石首、调关、监利等10个弯曲河段组成,且以塔市驿为界分石首和监利2个子河段,属于典型的蜿蜒型河道,局部河段主流摆动频繁，其中监利距离三峡大坝约345 km。本文计算断面为荆42(沙市站)和荆144(监利站),分别位于沙市(荆25—荆52)、监利(荆136—荆185)河段,且两河段均为沙质河床,床沙组成多为中细沙,其中值粒径分别介于0.20～0.27 mm、0.17～0.21 mm[17]。

图1 荆江河段示意Fig.1 Sketch of the Jingjiang reach

1.2 水沙变化及河床冲淤情况

分别点绘沙市、监利河段1998—2017年水量、沙量逐年(水文年)变化过程曲线(图2)。可以看出,2003—2017年,沙市河段年平均水量相比于三峡工程运行之前减少10.2%,沙量减少近85.9%；监利河段年平均水量减少7.3%,沙量减少近77.9%。总之,三峡工程运行后,下泄沙量大幅减少,使得水流经常处于严重次饱和状态,荆江河段往往会发生长距离长历时的河床冲刷,从而使含沙量沿程逐渐恢复。

图2 1998—2017年沙市及监利断面水量、沙量变化过程Fig.2 Temporal variations in water volume and sediment load at Shashi and Jianli during the period 1998—2017图3 沙市及监利河段累计冲淤量变化过程Fig.3 Accumulated volumes of channel evolution in subreaches of Shashi and Jianli

根据夏军强等[18]、Zhou等[19]的研究,荆江河段蓄水前(1975—2002年)平滩河槽年均冲刷量为0.14亿m3/a,蓄水后(2002—2012年)增加为0.37亿 m3/a,远大于蓄水前冲刷量。如图3所示(冲刷为正),三峡工程运行后,沙市、监利河段累计冲刷量(平滩河槽)逐年增加,其中沙市河段年均冲刷量为0.167亿m3,监利河段为0.131亿 m3,河床处于持续冲刷状态；尤其2008年以来,沙市、监利河段河床冲刷均有所加剧,其中沙市河段较监利河段明显。

2 分组悬沙恢复饱和系数

2.1 状态转移概率

在次饱和水流条件下,上游来水不断冲刷河床,为了形成稳定的抗冲层(粗化层),床沙组成不断调整。在床沙调整过程中,床沙(状态1)、推移质(状态2)、悬移质(状态3)泥沙的运动状态可以相互转换,这一状态转移过程称为Markov随机过程[20],而由某一运动状态转移到另一运动状态的概率即为状态转移概率。在复杂流场中,床沙交换过程具有随机性,其内在机理需要深入探究[21]。一些学者考虑了床沙与推移质或者悬移质之间的交换过程,提出了适用于均匀沙或输沙平衡状态的两态模型[22-25],但该模型在天然河流应用中存在局限性。为了弥补已有研究的不足,Kuai和Tsai[26]提出了同时考虑床沙、推移质及悬移质的三态模型,其转移概率矩阵εi可以表示为

(1)

在此基础上,Li等[20]对Kuai和Tsai[26]提出的Markov非齐次离散模型进行修正,得到基于非均匀沙隐蔽效应的三态转移概率矩阵为

(2)

式中:i为非均匀沙分组粒径；εxyi为状态转移概率,表示泥沙颗粒由x状态转移到y状态的概率；PTi为总概率,且PTi=PRi+PSi,PSi、PRi分别为起悬、滚动概率(与非均匀沙隐蔽效应有关)；1-P2i、1-P3i分别为止动概率、止悬概率,可按Li等[20]提出的方法进行计算。

2.2 紊流中泥沙的落距

如图4所示,假定运动泥沙颗粒可由悬移状态直接转变为床沙(忽略推移质运动的影响),紊流中泥沙落距(x0)定义为悬移质泥沙在沉降过程中的运动步长或距离,即运动泥沙起始点至床面接触点的水平距离。因此,运动泥沙的落距可由式(3)确定:

图4 泥沙颗粒在紊流中的沉降轨迹Fig.4 Settling track of a sediment particle

(3)

式中:V(y)为紊流中纵向水流流速；Ud(y)为颗粒沉降速度,即为(ωi-uy)的数学期望E(ωi-uy),ωi为沉速,uy为水流垂向脉动流速。

已知紊流中流速服从对数分布[7],泥沙落距可以表示为

(4)

式中:h′i为泥沙在紊流中的平均悬浮高；q(h′i)为自底部河床至平均悬浮高处的单宽流量。

参考韩其为和陈绪坚[6]的研究成果,q(h′i)服从对数分布,从而得到:

(5)

式中:η′i=h′i/h为泥沙在紊流中的相对悬浮高；u*为水流摩阻流速；V为垂线平均流速；κ为卡门常数,一般取κ=0.4。

(6)

将(6)式代入(5)式可得:

(7)

在式(7)中,用K表示泥沙落距影响系数,则式(7)可以进一步表示为

(8)

式(8)中,由于η′i=h′i/h≤1,所以K<1；C为非饱和调整系数,C>1表示冲刷,C<1表示淤积；C=1表示冲淤平衡。

2.3 分组悬沙恢复饱和系数

根据泥沙运动统计理论,可以得到第i粒径组悬沙落淤到床面的泥沙颗粒数N1i,以及从床面冲起的泥沙颗粒N2i；当N1i=N2i时,得到含沙量沿程变化方程为[6]

dSi/dx=-μ1ε21i(1-ε13i)(Si-S*i)

(9)

式中,Si为含沙量；S*i为水流挟沙力；ε21i为止滚概率,且ε21i=1-P2i；1-ε13i为止悬概率[21]；μ1=1/l1i,l1i为悬移质单步运动距离[7],可由式(10)确定。

(10)

式中:Uu(y)为悬移质颗粒上升的平均速度；Ud(y)为悬移质颗粒下降的平均速度；在紊流中,q(h′i)=q(η′i)。

根据韩其为和陈绪坚假定[6],在式(9)中:

αi/x0=μ1ε21i(1-ε13i)

(11)

式中:αi为恢复饱和系数；层流条件下x0=q/ωi[6],紊流条件下可由(5)式计算。

将式(11)代入式(9),可得泥沙连续方程为

(12)

可以看出,αi越大,含沙量沿程变化速率dS/dx越大,且含沙量恢复到水流挟沙力的速度越快。

联立式(8)—式(11),可得紊流条件下分组泥沙恢复饱和系数的计算关系式为

(13)

同时,也可将改进的恢复饱和系数计算方法应用于河床变形方程:

(14)

式中:ρ′为床沙干密度；y0为床面高程。

在计算天然河流悬沙恢复饱和系数时要结合实际水沙条件。如图5所示,长江中游荆江河段推移质较少,在输沙过程中所占比重不足10%；其中，沙市站推移质输沙量约占总输沙量的8.45%,监利站推移质输沙量占总输沙量的6.97%左右,泥沙输移的主要形式为悬移质。另外,床沙组成也会很大程度上影响悬沙恢复过程,即分组悬沙在沿程恢复过程中能否从床沙得到补给(床沙中是否含有相应粒径组的泥沙)。

图5 2003—2019年输沙率变化过程Fig.5 Temporal variations of sediment discharges at Shashi and Jianli in 2003—2019

由状态转移概率的定义可知,三态(床沙、推移质、悬移质)转移概率矩阵εi中,ε13i为泥沙起悬概率、ε21i为止滚概率,泥沙止动概率(既不起悬也不起滚)为ε21i(1-ε13i)。不考虑推移质的影响时,三态转移概率矩阵变为二维矩阵(床沙、悬移质),泥沙止动概率变为1-ε13i；也就是说,该情况下床沙只能通过悬移质方式进行输移,而床面落淤泥沙的来源也仅有悬移质泥沙。因此,当不考虑床沙交换过程中推移质的影响,即忽略床沙转变为推移质的概率(止滚概率ε21i),且考虑床沙组成是否具备补给悬沙的条件时,式(13)可以进一步表示为

(15)

式中:Pai为床沙中某粒径组泥沙所占比重。

既不考虑床沙交换过程中推移质的影响,也不考虑床沙组成是否具备补给悬沙的条件时,式(15)可以表示为

(16)

由式(13)可以看出,恢复饱和系数与止动概率、止悬概率、沉速及悬移质颗粒上升或下降的速度等因素有关；且与止滚概率及止悬概率成正比,与悬移质颗粒上升或下降速度亦成正比。此外,式(13)不用考虑非饱和调整系数(冲刷、淤积、平衡)的影响,只与一般的水流或泥沙参数有关,便于计算和应用,这也是本文式(13)与文献[6]公式的主要区别。本文通过式(13)计算了不同粒径组αi与泥沙粒径(d)、悬浮指标(Z)及摩阻流速(u*)之间的变化关系,如图6、图7所示。

图6表示在不同ω/u*(沉速与摩阻流速之比)取值情况下,αi随泥沙粒径的变化曲线。由图6可以看出,当ω/u*≤c时(0.001≤c≤0.01),αi随着泥沙粒径先减小,然后趋于某一定值(≈0.001)；当ω/u*>c时,αi先趋于某一定值,然后随着泥沙粒径的增大逐渐减小,最后亦趋于定值。当ω/u*一定的情况下,图6之所以出现拐点,是因为泥沙颗粒大于或者小于某粒径时,泥沙粒径已不是影响αi的主要因素；这时,αi主要受悬浮指标的影响,图7也可以解释这一现象。

图6 恢复饱和系数与泥沙粒径的变化关系Fig.6 Relationship between αi and d图7 恢复饱和系数与悬浮指标的变化关系Fig.7 Relationship between αi and Z

图7表示在不同泥沙粒径时,αi随Z的变化曲线。由图7可以看出,当泥沙粒径d<0.1 mm时,随着Z的增大,αi先逐渐增大,再趋于定值(≈0.3),最后当悬浮指标大于2时逐渐减小；当泥沙粒径d>0.1 mm时,αi随着Z的增大先增大后减小(Z=5为拐点)。在图7中,αi之所以会随着随着Z的增大先增大后减小,是因为本文考虑了床面泥沙颗粒间的隐蔽效应对泥沙起动、起悬的影响(隐蔽效应及Z同时影响αi)；在摩阻流速一定时,泥沙沉速越小,说明对应的泥沙粒径越小,该粒径组床沙所受到的隐蔽效应越大,泥沙颗粒就越难起动或起悬。

3 公式应用

为说明荆江河段各水文年内不同时期分组悬沙恢复饱和系数的取值情况,尤其是αi的沿程变化特点,将沙市、监利河段床沙、悬沙组成按粒径范围分为细沙(d<0.031 mm)、中沙(0.031 mm0.125 mm)3组；将沙市、监利站某一水文年按水位变化过程划分为枯水期、涨水期、洪水期及落水期等4个计算时段,其中2017年各时段划分结果如图8所示。结合河床平均高程变化过程可知,沙市站在4个时段的冲淤情况为冲淤交替、先淤积后平衡、冲淤交替及冲刷,监利站冲淤情况大致为先淤后冲、先淤积后平衡、冲淤交替及冲刷。由图8可以看出,沙市、监利站4个计算时段对应的时间点分别为:枯水期(1—2月)、涨水期(3—5月)、洪水期(6—10月)、落水期(11—12月)。

图8 沙市站及监利站2017水文年时段划分Fig.8 Division of the 2017 hydrological year at Shashi and Jianli

结合沙市、监利站2017年实测流量过程、床沙级配、含沙量及悬移质级配资料,运用式(11)分别计算荆42(沙市站)和荆144(监利站)断面分组悬沙在枯水期、涨水期、洪水期及落水期对应的恢复饱和系数取值及其变化特点。为验证公式的合理性及计算精度,在运用式(15)和式(16)计算αi时,分为2种情况:① 不考虑床沙组成的影响,忽略推移质的影响,即式(16)中不包含床沙组成项(Pai)；② 考虑床沙组成的影响,忽略推移质的影响,即式(15),计算结果如表1所示。其中,流速(U)、水深(h)为对应计算时段的水流参数平均值。

表1 沙市、监利站2017年分组悬沙恢复饱和系数计算结果Table 1 Calculation of saturation recovery coefficient of the grouped suspended sediment at Shashi and Jianli in 2017

由表1可以看出,在情况①下,荆42和荆144断面计算的分组悬沙恢复饱和系数有如下规律:

(1)细沙>中沙>粗沙,即恢复饱和系数随着泥沙粒径的增大而减小。

(2)αi在各计算时段的变化过程对悬沙分组具有敏感性；细沙对应的恢复饱和系数在各计算时段变化很小,即恢复饱和系数与划分的计算时段关系较小,并趋于某一常数；中沙和粗沙对应的恢复饱和系数在各计算时段的变化较细沙明显,即恢复饱和系数与划分的计算时段关系较好,荆42断面各时段恢复饱和系数表现为:洪水期>涨水期>枯水期>落水期,荆144断面大致表现为:洪水期>落水期>涨水期>枯水期,但粗沙对应的恢复饱和系数变化较中沙明显,同时也体现出冲刷时段的恢复饱和系数较淤积时大的规律。

(3)除落水期外,荆42断面在各计算时段得到的αi大于荆144断面,这是因为在不考虑床沙组成影响时,距离三峡工程越近,水流的不饱和程度越高,水流对河床的冲刷越剧烈,床沙对悬沙的补给就越多。此时,水流条件及泥沙粒径是决定悬沙沿程恢复快慢的主要因素。

(4)在不考虑床沙组成时,荆42和荆144断面αi的计算值都在0.12～0.27范围内变化。

本文还运用河床变形方程(式(14))计算了荆42和荆144断面分组悬沙在枯水期、涨水期、洪水期及落水期对应的恢复饱和系数,并与情况② 的计算结果进行对比,如表1所示。可以看出,除中沙对应的个别计算值外,运用式(15)和式(14)计算的细沙和中沙各时段的恢复饱和系数近似相等,且2种方法计算的粗沙各时段的恢复饱和系数符合较好,但式(14)个别计算值为负,与实际不符(式(14)计算值应大于或等于0)。式(15)计算的分组悬沙恢复饱和系数随泥沙粒径的增大而减小,除细沙外,中沙和粗沙在各时段的恢复饱和系数变化规律与式(15)计算结果相同；而式(14)的个别计算结果为负值,无明显变化规律,随机性较强。另外,由于在计算恢复饱和系数时,没有考虑冲泻质的影响,所以用式(15)和式(14)2种方法计算的细沙恢复饱和系数都为0,但这并不代表该河段在此时刻不能得到细沙的补给。造成这种计算结果的原因,可能是测量床沙级配时取样位置及取样深度的选取存在随机误差或者判定冲泻质的代表粒径偏大。

在考虑床沙组成时,荆42和荆144断面计算的分组悬沙恢复饱和系数有如下规律:

(1)荆42和荆144断面αi计算值分别在0.000 3～0.171 8、0.003 5～0.157 9范围内变化(细沙除外)；且细沙、中沙、粗沙对应的αi存在关系:细沙<中沙<粗沙,即恢复饱和系数随着泥沙粒径的增大而增大,与情况①得出的结论相反；这里恢复饱和系数随着泥沙粒径的增大而增大的结果只是一种表面现象,不是因为泥沙粒径变大导致恢复饱和系数增大,而是因为在考虑床沙组成时,床沙会随着次饱和水流的不断冲刷发生粗化,造成荆42和荆144断面2017年实测床沙组成中细沙、中沙所占比重很小,粗沙的比重较大(沙市站约为99%,监利站大于93%),进而导致式(15)中细沙和中沙对应的Pai很小,粗沙对应的床沙组成项Pai接近1；所以在其他条件不变的情况下,荆42和荆144断面中细沙、中沙对应的αi计算值很小(细沙为0),而粗沙对应的αi计算值较大,且与式(16)的计算结果相比变化很小。

(2)除落水期外,荆42断面粗沙在各计算时段得到的αi同样大于荆144断面,与情况①结论相同；但中沙的计算结果与粗沙相反,即荆42断面中沙在各计算时段得到的αi小于荆144断面。之所以出现这种现象,是因为考虑床沙组成的影响后,床沙组成(粗化程度)成为决定悬沙沿程恢复快慢的主要因素,而荆42的床沙粗化程度大于荆144断面，即荆42断面粗沙所占比重大于荆144断面,同时荆42断面中沙所占比重却小于荆144断面。

(3)αi在各计算时段的变化过程对悬沙分组同样具有敏感性；细沙对应的恢复饱和系数在各计算时段近似为0；中沙和粗沙对应的恢复饱和系数在各计算时段的变化较为明显,即恢复饱和系数与划分的计算时段关系较好,得出的结论与情况①相同,即荆42断面各时段恢复饱和系数表现为:洪水期>涨水期>枯水期>落水期,荆144断面大致表现为:洪水期>落水期>涨水期>枯水期。另外,结合实测水沙资料分析发现,三峡水库运行后荆江河段各粒径组悬沙含沙量在该河段都得到了不同程度的恢复,但不同粒径组悬沙恢复特点不同,其中细沙和中沙沿程恢复缓慢,而粗沙沿程恢复较快,并且荆42断面的恢复速度大于荆144断面,这也主要与该河段床沙组成有关。

另外,对比情况①和情况②可知,在情况②下,荆42断面细沙和中沙在各计算时段的恢复饱和系数近似为0,与情况①计算结果相比变化较大,粗沙计算结果与情况①近似相等；荆144断面细沙在各计算时段的恢复饱和系数近似为0,中沙很小,两者计算结果与情况①计算结果相比也发生了显著变化,而粗沙计算结果与情况①相差不大。造成这种现象的原因有2个方面:① 在不考虑床沙组成时,影响分组恢复饱和系数的主要因素为泥沙粒径、悬浮指标及止动概率、止悬概率等,而荆42和荆144断面的水流计算参数差别不大。② 当考虑了床沙组成的影响后,Pai成为决定分组恢复饱和系数大小的主要因素,而荆42和荆144断面实测床沙级配中细沙、中沙成分较少(主要由粗沙组成),但悬移质泥沙主要由细沙和中沙组成[27]。床沙只能补给悬沙中的粗沙成分,而不能为悬沙中细沙和中沙的沿程恢复提供补给作用,其中沙市、监利站2017年部分计算时段对应的床沙及悬沙级配如图9所示。

图9 沙市、监利站2017年床沙及悬沙级配Fig.9 Compositions of bed material and suspended load during different stages in 2017

4 结 论

基于Markov随机过程及泥沙运动理论,推导出非均匀悬沙的泥沙落距表达式,进而修正了分组悬沙恢复饱和系数公式,该方法可以同时考虑推移质及床沙组成对悬沙恢复过程的影响。得到的主要结论有:

(1)本文提出的恢复饱和系数公式主要与泥沙粒径、悬浮指标、床沙组成及止动概率、止悬概率等因素有关,无需考虑非饱和调整系数的影响；不考虑床沙组成时,沙市站和监利站分组悬沙恢复饱和系数计算值均为0.12～0.27,考虑床沙组成时分别为0.000 3～0.171 8和0.003 5～0.157 9,床沙组成对分组悬沙恢复饱和系数的影响较大。

(2)三峡水库运行后荆江河段各粒径组悬沙含沙量在该河段都得到了不同程度的恢复,但不同粒径组悬沙恢复特点不同,由于沿程床沙不段粗化,床沙组成中细沙和中沙成分不断减少,从而造成细沙、中沙对应的恢复饱和系数较小,而粗沙对应的恢复饱和系数较大；除落水期外,沙市站悬沙沿程恢复速度小于监利站。