三角洲是河流入湖(海)时,所携带的泥沙落淤沉积形成的沉积体.这个运动过程中挟沙水流经历了一个湍动射流过程,即无侧限水流的扩散过程,并在底部摩擦力和横向扩散的作用下向受水体方向减速[1-2],三角洲堆积体是沉积物沉降的结果.三角洲在地理位置上有独特的优势,加之三角洲沉积可以带来多种天然能源,故深入了解三角洲形成演变的过程和机理具有重要的意义.
对于三角洲形成演变过程的研究大致可以分为两类,一类为实验研究,另一类为数学模型研究.三角洲实验研究可以在实验室条件下塑造出和自然三角洲相似的特征过程,是研究三角洲形成演变的有效方式.Kim、Jerolmack[3]和Hoyal、Sheets[4]从实验角度突出了下游边界条件对三角洲动力演变的影响.Xin等[5]在实验过程中总结了上游来沙量和下游湖水位对入湖三角洲形成演变的影响.有学者引入三角洲自生行为来解释相似条件下三角洲演变差异的现象[6].三角洲数学模型研究是建立在水流泥沙运动基本理论基础上的,其机理明确.三角洲二维数学模型最早由Scheidegger[7]和Culling[8]提出.Price[9]阐述了一种基于随机游走理论的模型,并借此模拟了冲积扇的沉积过程.韩其为[10-11]根据不平衡输沙方程,推导得到了悬沙淤积的三角洲淤积百分数、前坡淤积百分数等三角洲形态特征的理论表达式.经典的扩散输运模型为模拟沉积体系中泥沙输移提供了可靠的方法[12-13].Rajeev等[14]研究了岸线运动的数学模型,提出了河流三角洲沉积过程中动边界问题的近似解法.Lorenzo-Trueba等[15]提出了一种几何模型来研究基准面变化在河控三角洲体系的岸线动态演变过程中的作用.
射流理论被广泛应用于描述河流进入静止状态水体所经历的过程[16-19].数学上,三角洲射流水动力可以由时均、非均匀、不可压缩的三维湍流方程来描述.而在某些特定的情况下,三维方程可以简化为二维方程并仍能精确地描述物理现象.浅水条件下的湍动射流可以由Navier-Stokes方程描述,方程为时均且深度平均[20].沿纵向速度断面自相似性是经典湍动射流的基础,实验证明自由射流情况下相似函数是成立的[21].Wang[16]将平面湍动射流理论应用于解释河口处的水沙现象,将连续方程和动量方程与泥沙扩散方程结合,利用相似函数得到了流速和含沙量断面的解析解.Muto和Steel[22]通过理论推导和现场数据验证,得出当下游相对水面以恒定速率持续上升,同时其他条件保持不变时,三角洲的岸线将不可避免的后退,三角洲前缘的运动速率将随时间变化的结论.Parker和Muto[23]通过一维数学模型计算了三角洲岸线在海平面上升时的变化.Edmonds等[24]建立了一个几何模型来计算以顶积层为主的三角洲在不同冲积和受水条件下顶积层和前积层厚度,以预测三角洲的类型.
尽管如此,我们对于入湖三角洲形成过程的机理和规律的认识尚未明晰,从理论角度探究这一过程有助于我们了解影响其发展演变的因素、深入认识入湖三角洲形成演变机制.
三角洲形成过程中,水流挟带泥沙运动至湖区,动力下降导致泥沙沉降淤积产生堆积体,逐渐形成入湖三角洲.这个过程相当于流体从面积较小的断面射入另一种流体的射流运动.在射出孔口进入下游环境后,初始阶段,射流的出流速度较高,主要依靠出射的初始动量来维持自身的继续运动,这一阶段动量对它的流动起支配作用;在之后的运动过程中射流将继续运动与扩散.以下将针对第一阶段的水沙运动特点,建立理论模型并求解.
假设水沙初始入射速度为u0,卷吸作用与掺混作用的流速小于u0的部分,即射流核心与静止液体之间的部分,称为射流边界层.纵向流速不等于零的射流区是以中心线为界的上下两个边界层的组合,如图1所示.
图1 三角洲初始段平面射流边界层示意图
Fig. 1 Schematic diagram of the plane jet boundary layer in the initial segment of a delta
浅水条件下,水深相对于水面区域尺度而言是小量,即垂向尺度远小于水平尺度,根据浅水问题特点,采用浑水积深浅水方程来描述水流运动过程:
连续方程
(1)
动量方程
(2)
(3)
其中,ux,uy分别为x,y向流速;h为水深;ρm,ρs,Δρ,ρ0分别为浑水密度、泥沙密度、泥沙与清水的密度差和床沙与饱和孔隙水的混合密度;y0为床沙底面高度;h0为淤积体高度;S为过水断面平均含沙量;εx,εy分别为x,y向湍流扩散系数;fx, fy分别为x,y向单位质量力;υ为动力黏滞系数;τsx,τsy分别为x,y向风应力;n为糙率.需要说明的是,在对方程时均化时,通常不考虑三阶关联项,在冲淤体积中不考虑对流运动和湍流脉动,脉动乘积项的时均值和可以分别看作是因湍流脉动在x,y轴方向引起的扩散,将其与时均参数相联系,设定
对上式进行如下假定:恒定流动条件下,忽略流速、水深对时间的变化率;浅水条件下,忽略水体中黏性应力的导数项和压力变化;忽略风应力和地转偏向力,只考虑床面上的黏性应力;认为浑水不可压缩ρm为常数.x方向的单位质量力为gsin θ,则连续方程、动量方程简化为
(4)
(5)
(6)
水深h是关于x的一维函数[16],三角洲形成初始阶段,当底床坡度较小时,水深h对x的导数项可忽略,且此时冲淤刚刚开始,水深h可视为常数,则上式可继续简化为
(7)
(8)
(9)
自相似现象是一种随时间发展的现象,通过彼此间的相似变换,可以得到系统性质在不同时刻、不同空间位置上的分布[26].自相似可以简化计算.对现象性质的描述,可以经所选自相似变量,将一群随机的经验数据点转化到一条曲线或一个面上.相似解法可以将偏微分方程简化为常微分方程.
相似解存在于某些具有特定物理对称性的流体问题中,对于所示的三角洲初始段形成过程,流速应存在相似解.实验资料亦表明,平面射流不同断面的流速分布具有相似性[12].对于图1,每一个x处都存在相应的速度分布,它们的共同特征是,射流边界层外缘速度为零,轴线处速度最大,射流各横截面上的速度分布具有相似性,即存在函数f,有
(10)
其中b为射流断面的特性半厚度,取为流速ux=um/e处的y值,根据射流边界沿直线扩散的特性,b=εx,前人实验结果得到ε=0.154[27].
根据实验资料和湍流随机性质的考虑,上式中的函数多采用Gauss正态分布的形式[28],其值为
(11)
设η=-y2/b2,令φ(η)=exp(η),则上式可表达为
ux=umφ(η).
(12)
引入流函数ψ,根据流函数定义有ux=∂ψ/∂y,将上式代入可得
(13)
令上式中由流函数定义可得
(14)
根据射流对称性,射流轴线处动量方程可简化为
(15)
代入上式可得
(16)
式(15)两边积分可得
(17)
即
(18)
其中C为常数.
由原点流速um=u0可得
代入上式得到轴线流速沿程变化关系式:
(19)
则流场分布公式为
ux=ume-y2/(ε2x2),
(20)
(21)
由以上两式可计算得
(22)
(23)
(24)
(25)
针对以推移质运动为主的三角洲,床面变形方程为
(26)
其中Pbl为第l组泥沙的组分, γs为泥沙的比重,Cm为浑水保持流体特性的最高含沙量,可按下式计算[29]:
Cm=0.755+0.222lg dl,
(27)
其中dl单位为mm;
qbl为第l组推移质输沙率,可按列维的方法进行计算[29]:
(28)
其中为第l组推移质泥沙的起动流速,其值为[29]
(29)
其中γ为水的比重, 列维计算公式的适用范围为dl=0.25~23 mm, h/dl=5~500, U/ul=1~3.5.则
(30)
由床面变形方程可得
(31)
将所求得第一阶段流场的分布代入,可求得淤积速度,即最初瞬间三角洲冲淤变化.在水力要素不变时,上式反映了淤积厚度沿程变化的特性.
本文采用实验室入湖三角洲实验数据进行验证.实验于天津大学内进行,实验装置包括流量和加沙速率都可调的水槽(长3.6 m,宽1.7 m,高0.2 m)、储水箱、实验区和循环系统,如图2所示.
通过调节水槽末端底部高度可以得到相应的实验坡度,本文中各组实验取坡度为1%.加沙系统向河道入口(河道宽度为6 cm)以预定的速率加沙,泥沙与水在河道入口处混合,挟沙水流将被输运到三角洲顶端处.水槽底部铺5 cm厚沙(本实验中的入口加沙和水槽底部铺沙都采用密度ρs=2 650 kg/m3,粒径大小d=0.3 mm的无黏性等粒径沙).
实验过程中,在设计实验流量下,泥沙主要以推移质形式运动,就像自然环境中以推移质运动为主的三角洲.实验模型参数参见表1.
图2 实验装置图
Fig. 2 The experimental device diagram
表1 实验模型参数设置
Table 1 Parameter setting of the experimental model
parametervaluedischarge Q/(cm3/s)100depth h/cm0.5slope J/%1particle size d/mm0.3roughness n0.02
实验通过固定精度为1 mm的激光测距仪到初始床面的距离来测量地形数据.本文中的三角洲实验用来观察时间和空间尺度上沉积过程的主要趋势.实验过程中实验条件的选择是基于前人实验、实验条件、自然三角洲特征以及预实验结果,保证可以展现三角洲演变过程中的沉积过程.
将实验数据和理论模型计算值直接进行比较是存在问题的,因为实验数据反映了水流泥沙的瞬时状态,即实验数据可以展现水沙的瞬间波动状态;而理论模型计算值模拟了水流运动和泥沙淤积在时间上的平均状态.同时,理论模型进行了合理简化后包括一些不能直接逐点测量的参数(例如糙率),实验过程中对于床面每一点随机因素不可避免(例如床面条件、来水来沙条件).因此,我们重点将模型预测的泥沙淤积主要趋势与实验数据所反映的趋势进行对比.为了保证实验过程中三角洲能比较充分的发展,在现有精度下采集的数据能全面反映地形变化趋势,同时将初期淤积形态对水流和泥沙运动的影响降低,选取实验1 h末冲淤实测值与模型计算值进行分析.1 h末的三角洲实测地形如图3所示.
实验进行1 h,堆积体在x方向推进了大约0.5 m,y方向上以x轴为对称轴上下各推进了大约0.2 m,故选定x方向0.5 m、y方向0.2 m的范围进行验证.x轴方向每0.1 m设定一个验证断面、y轴方向每0.05 m设定一个验证点,验证点布置示意图如图4所示.
图3 1 h末的三角洲实测地形
Fig. 3 Measured topography of the delta at the end of 1 h
注 为了解释图中的颜色,读者可以参考本文的电子网页版本.
图4 验证点布置示意图
Fig. 4 Layout of verification points
各断面冲淤z的计算值和实测值如图5所示.
从图5可以看出,各断面的冲淤趋势和特点不尽相同.由断面x=0.1 m的冲淤计算结果显示,此断面在x轴处的冲刷程度最大,在从y=0到y=±0.05 m的过程中zs/t的计算值有较快增长,之后向x轴两侧发展过程中zs/t的计算值增长缓慢,即从y=0到y=±0.2 m的过程中冲刷程度逐渐减小.由断面x=0.1 m的冲淤实测值可知,断面在y=0处的冲刷最剧烈,从y=0到y=±0.2 m的过程中冲刷程度逐渐减小,与计算结果吻合.
由断面x=0.2 m的冲淤计算结果可以看出, y=0处冲刷最剧烈,从y=0到y=±0.05 m的过程中zs/t值经历了快速增长, 从y=±0.05 m到y=±0.2 m的过程中地形变化速率增速逐渐放缓.由该断面的实测值可知, y=0处产生的冲刷最剧烈, 与计算结果一致,从y=0到y= ±0.2 m过程中zs值增速均匀,即冲刷程度均匀减小.虽然实测zs的极值点与计算值有所差异,但整体发展趋势基本相似.
断面x=0.3 m的计算结果显示,各验证点将受到冲刷,冲刷最剧烈发生在y=±0.05 m处,在从y=±0.05 m到y=±0.1 m的过程中冲刷程度迅速减弱.由断面实测结果可知,多数点处于冲刷状态,与计算结果比较一致.与计算结果相比,实测x=0.3 m断面冲刷最剧烈点向外延展,产生在y=±0.1 m处,计算冲刷最剧烈处(y=±0.05 m)成为了冲刷程度最小和冲刷程度最大点之间的过渡点.在y=±0.1 m到y=±0.2 m的过程中,冲刷程度逐渐减弱.计算值和实测值整体冲刷趋势都是 “先增强再减弱”.
从断面x=0.4 m处的计算结果可以看到, 各验证点将受到冲刷, 冲刷最剧烈发生在y= ±0.05 m处,在从y=±0.05 m到y=±0.2 m的过程中冲刷程度逐渐减弱.实测结果显示,在y=±0.05 m和y=0.1 m处产生冲刷,其余各点不同程度淤积,y=±0.05 m左右处冲刷最剧烈,与计算值一致.与计算值不同的是,实测值显示各验证点有冲刷也有少量淤积且断面不对称发育,从y=0到y=0.2 m的过程中,地形呈现“先淤后冲再淤”的趋势.计算值与实测值发展趋势基本一致.
(a) x=0.1 m
(b) x=0.2 m
(c) x=0.3 m
(d) x=0.4 m
(e) x=0.5 m
图5 各断面冲淤的计算值和实测值
Fig. 5 Calculated and measured values of erosion and deposition on each section
从断面x=0.5 m处的计算结果可知,各验证点将受冲刷,冲刷最剧烈处发生在y=±0.05 m处,在从y=±0.05 m到y=±0.2 m的过程中冲刷程度逐渐减弱.实测结果中,y=0处产生了于该断面来说较大的淤积.与x=0.4 m断面相比,该断面不对称发育程度加剧,zs的最小值产生在y=-0.05 m和y=0.2 m处.实测趋势与计算值发展趋势基本一致.
理论解和实验结果在趋势上较为一致,但在量级上有所差别,故有必要将理论解和实验结果与前人成果做对比研究.首先,利用上述实验数据,使用本文提出的计算方法与文献[27]中平面湍动射流的流场计算方法计算轴线上x=0.4~1 m范围内(文献[27]中的轴线流速针对主体段,实验验证中的初始段长约为0.3 m,故将计算范围设定为x=0.4~1 m)的理论冲淤速率,将计算结果画入图6.
图6 本文和文献[27]计算结果的对比 图7 本文与文献[30]中三角洲轴线地形数据对比
Fig. 6 The comparison of calculation results Fig. 7 The comparison of the delta axis topographic between this paper and ref. [27] data between this paper and ref. [30]
由图6可见,两种计算方法算得的x=0.4~1 m范围内轴线上的冲淤量结果属于同一量级且数值比较接近.
对比本文实验三角洲地形数据与文献[30]中实验三角洲对照组地形数据,二者水沙条件不同,但都是在实验室内模拟以推移质为主的入湖三角洲形成演变过程,对比两次实验x=0.1~0.5 m范围内的轴线淤积量,将结果画入图7.
由图7可见,两次实验的地形冲淤量在同一数量级,只是由于水沙条件不同,两组实验轴线上的淤积量数值有所不同.
基于以上对比,进一步分析可知,出现底床冲淤理论值变化幅度小,与实验值有量级上的区别,应是由于理论值较小,受限于实测精度,导致实测值量级较大.但将理论值和实测值的发展趋势加以对比,若较为一致,则能够说明理论方法的正确性.
综上,x=0.1~0.5 m断面zs变化率的计算值和实测值的发展趋势较为一致.其中,断面x=0.1~0.2 m范围内计算和实测趋势整体呈V形;x=0.3~0.5 m范围内计算和实测值整体呈W形;相较于计算值,实测地形发展较为平缓,x=0.3 m断面地形实测值出现了最深冲刷点外延现象,x=0.3~0.5 m断面出现了明显的不对称发育现象.
利用上述实验数据,分别使用本文提出的计算方法与文献[27]中的平面湍动射流的流场计算方法计算流场分布,使用相同的床面变形方程计算底床变形,比较两种方法计算的床面变形结果.文献[27]中的轴线流速针对主体段(x>0.3 m), 故分别使用两种计算方法计算x=0.4 m和x=0.5 m断面zs的变化速率,如图8所示.
由图8可知,两种算法下x=0.4 m和x=0.5 m断面zs/t值的发展趋势大体一致,主要区别在于冲刷最剧烈点的位置有所不同:文献[27]冲刷最剧烈点出现在x=±0.1 m处,本文的计算方法中冲刷最剧烈的点位于x=±0.05 m处.实验实测数据显示,x=0.4 m断面冲刷程度最大的点出现在y=0.05 m,其次为y=-0.05 m; x=0.5 m断面zs/t的最小值位于y=-0.05 m和y=0.2 m,其次为y=0.05 m.结合实验数据分析可知,本文的计算结果更加合理.
(a) x=0.4 m
(b) x=0.5 m
图8 本文与文献[27]方法计算的冲淤速率结果对比
Fig. 8 The comparison of calculated erosion and deposition rates between this paper and ref. [27]
本文根据三角洲初始段形成过程中水沙运动特点,以浑水浅水方程为基础,构建了三角洲初始段平面射流边界层理论模式.针对初始段形成过程水流运动特点,采用相似解解法,得到了三角洲形成初期的流场分布.利用河床演变通用数学模型推导得到了表征入湖三角洲初始段形态特征的理论表达式.通过实验验证,理论模型能够较好地反映三角洲形成初期的冲淤规律.验证过程中可以发现,相较于理论值,实测地形各点之间的过渡更平稳光滑,多个断面出现了不对称发育现象,这与入口附近的水流具有较大动能,水流湍动作用于不均匀分布的泥沙(具体到某一小撮的泥沙),放大了水沙随机性对三角洲形态表达的影响.
[1] JOSHI P B. Hydromechanics of tidal jets[J]. Journal of the Waterway, Port, Coastal and Ocean Division, 1982, 108(3): 239-253.
[2] ORTEGA-SNCHEZ M, LOSADA M A, BAQUERIZO A. A global model of a tidal jet including the effects of friction and bottom slope[J]. Journal of Hydraulic Research, 2008, 46(1): 80-86.
[3] KIM M, JEROLMACK D J. The pulse of calm fan deltas[J]. The Journal of Geology, 2008, 116(4): 315-330.
[4] HOYAL D C J D, SHEETS B A. Morphodynamic evolution of experimental cohesive deltas[J]. Journal of Geophysical Research Atmospheres, 2009, 114(F2). DOI: 10.1029/2007JF000882.
[5] XIN W Y, BAI Y C, XU H J. Experimental study on evolution of lacustrine shallow-water delta[J]. Catena, 2019, 182: 1-10.
[6] WHIPPLE K X, PARKER G, PAOLA C, et al. Channel dynamics, sediment transport, and the slope of alluvial fans: experimental study[J]. The Journal of Geology, 1998, 106(6): 677-694.
[7] SCHEIDEGGER A E. Hydraulic effects in geodynamics[J]. Geologie und Bauwesen, 1959, 25: 3-49.
[8] CULLING W E H. Analytical theory of erosion[J]. The Journal of Geology, 1960, 68(3): 336-344.
[9] PRICE W E. Simulation of alluvial fan deposition by a random walk model[J]. Water Resources Research, 1974, 10(2): 263-274.
[10] 韩其为. 论水库的三角洲淤积(一)[J]. 湖泊科学, 1995, 7(2): 107-118.(HAN Qiwei. Study on the deposition of reservoir delta (1)[J]. Journal of Lake Sciences, 1995, 7(2): 107-118.(in Chinese))
[11] 韩其为. 论水库的三角洲淤积(二)[J]. 湖泊科学, 1995, 7(3): 213-225.(HAN Qiwei. Study on the deposition of reservoir delta (2)[J]. Journal of Lake Sciences, 1995, 7(3): 213-225.(in Chinese))
[12] SWENSON J B, VOLLER V R, PAOLA C, et al. Fluviodeltaic sedimentation: a generalized Stefan problem[J]. European Journal of Applied Mathematics, 2000, 11: 433-452.
[13] CAPART H, BELLAL M, YOUNG D L. Self-similar evolution of semi-infinite alluvial channels with moving boundaries[J]. Journal of Sedimentary Research, 2007, 77: 13-22.
[14] RAJEEV, KUSHWAHA M S, KUMAR A. An approximate solution to a moving boundary problem with space-time fractional derivative in fluvio-deltaic sedimentation process[J]. Ain Shams Engineering Journal, 2013, 4(4): 889-895.
[15] LORENZO-TRUEBA J, VOLLER V R, PAOLA C. A geometric model for the dynamics of a fluvially dominated deltaic system under base-level change[J]. Computers & Geosciences, 2013, 53: 39-47.
[16] WANG F C. The dynamics of a river-bay delta system[J]. Journal of Geophysical Research: Oceans, 1984, 89(C5): 8054-8060.
[17] FALCINI F, JEROLMACK D J. A potential vorticity theory for the formation of elongate channels in river deltas and lakes[J]. Journal of Geophysical Research Atmospheres, 2010, 115(F4). DOI: 10.1029/2010JF001802.
[18] NARDIN W, MARIOTTI G, EDMONDS D A, et al. Growth of river mouth bars in sheltered bays in the presence of frontal waves[J]. Journal of Geophysical Research: Earth Surface, 2013, 118(2): 872-886.
[19] LEONARDI N, CANESTRELLI A, SUN T, et al. Effect of tides on mouth bar morphology and hydrodynamics[J]. Journal of Geophysical Research: Oceans, 2013, 118(9): 4169-4183.
[20] SCHLICHTING H, GERSTEN K. Boundary-Layer Theory[M]. Berlin: Springer Verlag, 1968.
[21] RAJARATNAM N, STALKER M J. Circular wall jets in coflowing streams[J]. Journal of the Hydraulics Division, 1982, 108(2): 187-198.
[22] MUTO T, STEEL R J. Retreat of the front in a prograding delta[J]. Geology, 1992, 20(11): 967-970.
[23] PARKER G, MUTO T. 1D numerical model of delta response to rising sea level[C]//3rd IAHR Symposium, River, Coastal and Estuarine Morphodynamics. Barcelona, Spain, 2003.
[24] EDMONDS D A, SHAW J B, MOHRIG D. Topset-dominated deltas: a new model for river delta stratigraphy[J]. Geology, 2011, 39(12): 1175-1178.
[25] 谢鉴衡. 河流模拟[M]. 北京: 水利电力出版社, 1990.(XIE Jianheng. The Simulation of River[M]. Beijing: Water Conservancy and Electric Power Press, 1990.(in Chinese))
[26] BARENBLATT G I, ISAAKOVICH B G. Scaling, Self-Similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics[M]. Cambridge: Cambridge University Press, 1996.
[27] 四川大学水力学与山区河流开发保护国家重点实验室. 水力学[M]. 北京: 高等教育出版社, 2016: 210-217.(State Key Laboratory of Hydraulics and Mountain River Engineering of Sichuan University. Hydraulics[M]. Beijing: Higher Education Press, 2016: 210-217.(in Chinese))
[28] 余常昭. 环境流体力学导论[M]. 北京: 清华大学出版社, 1992: 202-206.(YU Changzhao. Introduction to Environmental Hydrodynamics[M]. Beijing: Tsinghua University Press, 1992: 202-206.(in Chinese))
[29] 钱宁, 万兆惠. 泥沙运动力学[M]. 北京: 科学出版社, 1983.(QIAN Ning, WAN Zhaohui. Mechanics of Sediment Movement[M]. Beijing: Science Press, 1983.(in Chinese))
[30] 白玉川, 胡晓, 徐海珏. 入湖浅水三角洲形成过程实验模拟分析[J]. 水利学报, 2018, 49(5): 549-560.(BAI Yuchuan, HU Xiao, XU Haijue. Experimental analysis of the formation process of lacustrine shallow-water delta[J]. Journal of Hydraulic Engineering, 2018, 49(5): 549-560.(in Chinese))
白玉川(1967—),男,教授,博士,博士生导师(E-mail: ychbai@tju.edu.cn);
徐海珏(1977—),女,副教授,博士,硕士生导师(通讯作者. E-mail: xiaoxiaoxu_2004@163.com).