从20世纪50年代起我国沿海河口修建了大量挡潮闸以挡潮御卤、排洪蓄淡^[1].由于河口建闸后潮波变形加剧，河道内涨落潮过程的明显不对称形成了河道涨潮大于落潮的不平衡输沙，长时期关闭闸门使得闸下河道发生持续性淤积^[2-4]，严重影响河道防洪安全，降低其通航能力，限制河口资源开发和利用.近年来，随着我国经济的持续发展，淤积型河口与港口、航道的开发之间的矛盾渐显，因此实现河口、港口和水运航道一体发展的理想格局是河口建设发展的主要方向.

在现阶段，河口航道建设及维护主要依靠修筑导堤、岸线整治和机械清淤等工程措施.此类工程措施可通过改变河口地形打破原有水动力和输沙规律，改变河道冲淤特性^[5-6].一般而言，河道边滩的宽度为过流断面的50%时，可维持河道涨落潮相对平衡，一旦大量潮间带被浚深会使涨潮流增强，海相细颗粒泥沙将造成河道回淤^[7].Winterwerp等^[8]通过对欧洲河道内潮差分析发现，20世纪许多河道内的潮差在人类工程影响下呈现增大的态势，同时指出潮差增加与水深增加、河道长度变短及边滩面积减小有关.Maren等^[9]通过数值模拟的方法研究了河道疏浚和港口建设对Ems河口泥沙输运的影响，发现河道浚深是迫使泥沙向河口输运的主要因素.此外，不同的工程方案(如挖槽深度、宽度、长度)也使回淤率存在较大差异.刘杰等^[10]利用长江口深水航道回淤资料研究了北槽航道回淤时空变化特征及其对河槽地形的影响，发现北槽主槽地形冲淤调整引起的高差变化，是航道回淤重心年际间上提的主要原因.李大鸣等^[11]通过建立泥质河口二维水沙数学模型对宽浅型和深窄型两类局部清淤方案效果进行了分析，发现清淤宽度越小、深度越大，回淤率越高；清淤长度越长，回淤率越小.在相同的清淤条件下，宽浅型清淤方案与窄深型清淤方案相比，具有回淤量小、回淤率低的特点.综上可以看出，河道的冲淤特征与工程建设方案密切相关.然而，对于原本潮波反射剧烈、涨落潮不平衡的建闸河口，地形和岸线的改变引起的河道冲淤响应更为复杂和剧烈^[12-15]，因此，建闸河口整治工程对河床演变的影响需要进一步研究.

双龙河口位于河北省北部沿海，在嘴东工业园区与曹妃甸工业区交界处，是该海域重要的入海航道.双龙河口建闸后，河床10年淤高2~4 m，河道束窄变浅，严重降低了原有的河道通航能力，行洪能力降为原来的32%^[16].随着曹妃甸围海造陆工程和嘴东工业园区的发展，双龙河口战略位置日益凸显.为改善河口的水动力环境，提升河道通航等级，自2011年开展双龙河口整治工程.本文以双龙河口整治工程为例，利用数值模拟手段对一期、二期整治工程后挡潮闸下河道水动力调整及工程后河床演变进行计算分析，揭示建闸河口在整治工程影响下的演变规律，提出建闸河口整治及维护建议.

1 河口及工程概况

双龙河(图 1)属季节性河流，河长65 km，集水面积为443 km²，河道坡降(0.005%~0.02%)较平缓，底质以0.005~0.008 mm的粉砂黏土为主.河口潮流主要受渤海潮波系统控制，属不正规半日混合潮，平均潮差2.24 m.双龙河挡潮闸建于1978年，为6孔直升式闸门，闸口净宽48 m，闸门仅在汛期排洪或需排污水时提起，平时为常闭闸.建闸后，潮流成为双龙河口的控制性动力，涨潮历时缩短，落潮历时延长，涨潮流速大于落潮流速.由于当地海相来沙丰富，且涨潮输沙量大于落潮输沙量，闸下河道发生严重淤积，至1988年河床淤高2~4 m^[16]，河床-0.3 m等深线宽度由90 m缩窄至30 m(图 2a).为改善河口的水动力环境和提升河道通航等级，双龙河口整治工程自2011年开始实施，主要对挡潮闸下至入海口长度为4.09 km河段进行岸线整治、河道清淤及岸坡防护.一期工程(图 2b)于2013年竣工，完成闸下1.834 km范围内的整治.河道按0.05%的坡降进行清淤，河宽拓展至180 m，河口海域面积增加5.17 km².为营造亲水、绿色友好的海岸环境，河道右岸采用沉井搭建亲水平台，地面回填至4.65 m，亲水平台与地面采用1:10缓坡相连，作为景观用地；左岸采用1:4坡度与现状地面相连，岸坡采用植草护坡形式控制滩面泥沙启动，促使泥沙就地落淤.二期工程(图 2c)计划对一期工程下游至河口(距闸约4.09 km)范围内河段以0.05%的坡降继续清淤和岸线整治.

2 数学模型的建立及验证

采用丹麦水力学研究所(DHI)研发的MIKE 21 Flow Model模型^[17]中的水动力模块(HD)和黏性泥沙(MT)模块建立双龙河口潮流、泥沙输运及河床演变数学模型.模型通过有限体积法进行离散，显式欧拉法求解，具备计算速度快、易收敛、计算精度高等优点^[18].

2.1 基本控制方程 2.1.1 潮流模型

潮流模型建立在Navier-Stokes方程的基础上.在笛卡尔坐标系下，通过对三维水平动量方程和连续方程沿水深积分得到如下二维浅水方程：

$\frac{{\partial h}}{{\partial t}} + \frac{{\partial hu}}{{\partial x}} + \frac{{\partial hv}}{{\partial y}} = 0$

(1)

$\begin{array}{l} \frac{{\partial hu}}{{\partial t}} + \frac{{\partial h{u^2}}}{{\partial x}} + \frac{{\partial hvu}}{{\partial y}} = fvh - gh\frac{{\partial \eta }}{{\partial x}} + \frac{{{\tau _{{\rm{s}}x}} - {\tau _{{\rm{b}}x}}}}{\rho } + \\ \quad \frac{\partial }{x}\left( {h{T_{xx}}} \right) + \frac{\partial }{{\partial y}}\left( {h{T_{xy}}} \right) \end{array}$

(2)

$\begin{array}{l} \frac{{\partial hv}}{{\partial t}} + \frac{{\partial huv}}{{\partial x}} + \frac{{\partial h{v^2}}}{{\partial y}} = - fuh - gh\frac{{\partial \eta }}{{\partial y}} + \frac{{{\tau _{{\rm{s}}y}} - {\tau _{{\rm{b}}y}}}}{\rho } + \\ \quad \frac{\partial }{{\partial x}}\left( {h{T_{xy}}} \right) + \frac{\partial }{{\partial y}}\left( {h{T_{yy}}} \right) \end{array}$

(3)

式中：t为时间；x，y为笛卡尔坐标下纵横坐标；h为总水深，h=η+d，η为水位，d为静止水深；u，v分别为x，y方向的水深平均流速；f为科氏力系数；ρ为水的密度；τ_sx，τ_sy分别为x，y方向上的表面风应力；τ_bx，τ_by分别为x，y方向的底床切应力，${\tau _{{\rm{b}}x}} = \rho g\frac{{\sqrt {{u^2} + {v^2}} }}{{{M^2}{h^{\frac{1}{3}}}}}u,{\tau _{{\rm{b}}y}} = \rho g\frac{{\sqrt {{u^2} + {v^2}} }}{{{M^2}{h^{\frac{1}{3}}}}}v$; M为曼宁数；侧应力T_ij包括黏性摩擦和紊动摩擦等，通过水深平均流速沿梯度方向的变化计算如下：

$\left\{ \begin{array}{l} {T_{xx}} = 2A\left( {\frac{{\partial u}}{{\partial x}}} \right)\\ {T_{xy}} = A\left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}} \right)\\ {T_{yy}} = 2A\left( {\frac{{\partial v}}{{\partial y}}} \right) \end{array} \right.$

(4)

亚网格尺度涡流黏度A由Samagorinsky公式(5) 确定：

$A = c_{\rm{s}}^2{l^2}\sqrt {2{S_{ij}}{S_{ij}}} $

(5)

式中：c_s为常数，在本文中取0.28；l为特征长度；变形率S_ij根据下式计算：

${S_{ij}} = \frac{1}{2}\left( {\frac{{\partial {u_{\rm{i}}}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_{\rm{i}}}}}} \right)$

(6)

2.1.2 泥沙输运模型

双龙河口泥沙运动以悬移质为主，因此采用如下考虑淤积和侵蚀的对流扩散方程描述泥沙输运：

$\begin{array}{l} \frac{{\partial hc}}{{\partial t}} + \frac{{\partial huc}}{{\partial x}} + \frac{{\partial hvc}}{{\partial y}} = \frac{\partial }{{\partial x}}\left( {h{\varepsilon _x}\frac{{\partial c}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {h{\varepsilon _y}\frac{{\partial c}}{{\partial y}}} \right) + \\ {Q_{\rm{L}}}{C_{\rm{L}}} - S \end{array}$

(7)

$S = \left\{ {\begin{array}{*{20}{c}} {{\omega _{\rm{s}}}{c_{\rm{b}}}\left( {1 - \frac{{{\tau _{\rm{b}}}}}{{{\tau _{{\rm{cd}}}}}}} \right)}&{{\tau _{\rm{b}}} ＜ {\tau _{{\rm{cd}}}},{\rm{ }}河床淤积}\\ 0&{{\tau _{{\rm{cd}}}} \le {\tau _{\rm{b}}} \le {\tau _{{\rm{ce}}}},河床不冲不淤}\\ { - E{{\left( {\frac{{{\tau _{\rm{b}}}}}{{{\tau _{{\rm{ce}}}}}} - 1} \right)}^n}}&{{\tau _{\rm{b}}} ＞ {\tau _{{\rm{ce}}}},{\rm{ }}河床冲刷} \end{array}} \right.$

(8)

式中：c为沿水深平均含沙量；ε_x，ε_y分别为x，y向的泥沙扩散系数；S为淤积/冲刷项；Q_L为水平方向单位面积的源流量；C_L为源的含沙量；ω_s为泥沙沉速；c_b为近底含沙量；τ_b为底床切应力；τ_cd为临界淤积切应力；τ_ce为临界侵蚀切应力；E为侵蚀系数；n为侵蚀指数.

2.1.3 河床演变模型

河床演变方程为

$\rho _{\text{S}}^{'}\frac{\partial {{\eta }_{\text{s}}}}{\partial t}=S$

(9)

式中：ρ′_S为悬沙干密度；η_s为悬沙造成的河床淤积厚度，正值表示淤积，负值表示冲刷.

2.2 计算范围及网格

为了保证计算的准确性及连续性，且计算边界不受曹妃甸围海造陆工程的影响，模型计算范围自双龙河口挡潮闸向海延伸40 km，包括整个双龙河口及曹妃甸在内的4 500 km²海域(图 1).模型采用非结构网格对计算区域进行渐变剖分，网格结点总数为15 337个，网格单元总数为28 473个，河道横截面方向网格数约为20个，最高分辨率在闸下工程区域，为10 m，至外海网格逐渐扩至2 800 m.

2.3 边界条件、初始条件及参数选取 2.3.1 边界条件和初始条件

模型开边界包括东、西、南外海开边界和闸控边界.东、西外海开边界的走向基本与主流方向垂直，从而增强模型的稳定性.潮流模型外海开边界条件为渤海潮流模型^[19](边界为大连至烟台)提供的潮位过程，由于挡潮闸常年关闭，闸控边界流量取零.东、西、南外海开边界悬沙浓度根据渤海多年平均含沙量分布进行线性插值给定，分别为0.005~0.010 kg·m^-3, 0.010~0.020 kg·m^-3, 0.005~0.010 kg·m^-3.潮位初始条件取开始时刻的边界潮位平均值，初始流速为零，初始含沙量为0.01 kg·m^-3.

2.3.2 关键参数选取

随着潮位涨落变化，部分区域露滩，模型通过干湿网格法界定动边界，干水深h_dry、淹没水深h_flood和湿水深h_wet分别取0.005 m, 0.05 m和0.1 m，M根据底部泥沙粒径和水深分布取值27~97 m^1/3·s^-1，Samagorinsky模型中常数c_s取0.28，以此来计算涡黏性系数.ω_s参照泥沙粒径分布取平均值0.000 5 m·s^-1；临界侵蚀切应力τ_ce根据唐存本公式^[20]估算后经率定调整取0.51~1.00 N·m^-2，τ_cd取0.22~0.44 N·m^-2；E取1.5×10^-5 kg·m^-2·s^-1；n取1；床面糙率k_n取0.001 m；床层泥沙干密度ρ′_s取具代表性的恒定值811 kg·m^-3.对河床冲淤模型进行技术约减，模型采用取值为12的地形演变加速系数，外海开边界条件为典型中潮差潮位过程，经计算1个月后即得全年的冲淤演变.

2.4 模型验证

采用2012年12月18至19日S₁和S₂点(图 1)流速和流向过程以及S₂点实测潮位对潮流模型进行验证，从图 3的验证结果可见，S₂点的模拟潮位过程与实测潮位过程一致.同时，涨潮过程显示双龙河口海域属于不正规半日潮，涨潮历时约5 h，落潮历时约7 h.S₁和S₂点的模拟流速、流向过程与实测值在相位和数值上都基本吻合.2个点的潮流过程显示，潮流由外海传入闸下河道发生较大变形，具有涨潮历时缩短、落潮历时增长的涨落潮不对称特点.

泥沙输运模型验证采用2012年10月15至16日曹妃甸海域S₃和S₄点垂向平均悬沙浓度.对2010年9月及2011年9月实测地形数据进行差值计算，得到闸下3.5 km河道深泓线一年的沿程冲淤变化，以此用作河床冲淤模型验证.图 4、图 5分别为S₃和S₄点悬沙浓度验证过程和河床冲淤验证，结果显示计算悬沙浓度与实测过程基本一致，河床纵剖面年厚度变化与实测值吻合较好，表明数学模型选取的泥沙计算参数较合理，能反映双龙河口泥沙输运过程和河床演变规律.

进一步使用Willmott^[21]所提出的统计学方法对模型进行评价，计算所得的潮流模型、泥沙输运及河床冲淤模型效率系数分别为0.70~0.96(极好)、0.50~0.57(非常好)、0.90(极好)，整体而言模型计算误差在可接受范围内，可用于建闸河口整治工程与河床演变特征研究.

3 结果分析

双龙河口挡潮闸闸下引河从长度上属短引河，潮流为其控制性动力.在整治工程前，河道滩槽分明，平面形态不但受主槽落潮流控制，也受滩面归槽水影响.整治工程后，工程河段水流不受滩面的束缚，河道水流和冲淤发生改变.为探究整治工程对闸下水动力调整和河床演变的影响，运用已经验证的合理的数学模型计算分析整治工程前后闸下河道的潮流变化和河床演变特征.工程后岸线采用一期、二期工程后形成的新岸线，修改工程河段地形以反映河道清淤效果，工程后河道横截面方向网格数略有调整，约为15个，其他计算参数与模型验证时相同.

3.1 潮流响应分析

地形岸线变化直接改变河道水流流态，图 6为3种工况下双龙河口典型时刻流场.2011年工程前水流主要流经河道主槽，总体表现为顺岸往复流，涨潮流向西北，落潮流向东南，主槽涨急流速大于落急流速.一期工程后，未整治河段流态并未发生较大改变，主流向仍保持与主槽走向一致，而工程河段由于断面单一、规整，涨急和落急时刻水流变得比较平顺.涨急时刻，潮流从外海进入河道，由于整治工程河道纳潮量增加，水流以较大流速从下段河道主槽流向上段宽阔工程河道，在地形突然陡降展宽的区域，水流呈喷射状涌入整治河段并在河道右岸产生逆时针回流.之后主流继续向闸下推进，虽然流态平顺，但流速因边滩大量清除、过水断面增加，流速发生大幅度降低.落急时刻，工程河段较弱的水流经河槽汇水进入下游河段，汇水段主槽形成高流速区.河道弯段凸岸有浅滩分布，主流在凹岸一侧更为集中，水流经弯段转向后，受河槽束水作用主槽流速增大，断面流态紊乱.二期工程后，虽涨急流速仍大于落潮流速，但河道流态十分规整，从河口至闸下自然过渡.

潮流涨落潮过程中的时均流速是衡量潮流强度的重要特征，也是影响泥沙输运的主要因素，它是潮差、潮周期和下垫面条件的综合反映.表 1和表 2分别显示了工程前后河道各断面的特征流速及其变化率，表中V和V_max分别为平均流速和最大流速，V₁和V₂分别为涨潮和落潮平均流速，α为涨潮与落潮平均流速的比值，P₁和P₂分别为一期工程和二期工程后(与工程前相比)各特征流速的变化率.结合图 6可以看出：一期工程后，受河道浚深拓宽影响，工程河段内的D₁断面平均流速减小60.0%，但因工程河段地形平整、水深增加，潮波传播过程中耗能减小，造成断面涨落潮平均流速比由3.90增至5.00.工程下游河段断面平均流速增长46.7%~76.1%，断面涨落潮平均流速比由一期工程前的1.26~5.70减至0.87~3.01.这是受纳潮量增加和边滩综合影响导致的非工程河段涨落潮流速不对称程度减缓.二期工程后，整个河口地形按固定坡降均匀变化，河口纳潮量得到显著增加，在地形变化和纳潮量增加的双重作用下，河道除拓宽段(D₃断面)外平均流速均有较大的提升，特别落潮流速增长幅度较大，使得涨落潮流速差距明显减小，涨落潮平均流速比值维持在1.13~1.27.其中位于河道弯段的D₂断面落潮流速较工程前可增加862.5%，涨落潮平均流速比值由工程前的5.70降至1.33，流速增加及涨落潮流速不对称性锐减可有效降低该段淤积强度，有利于维持通航水深.综合以上分析，一期工程造成工程段涨落潮流速不对称性加剧，易导致该段输沙不平衡增强，而二期工程后河道纳潮量显著增加，涨落潮流速不对称性减小，输沙不平衡程度减小.

3.2 河床冲淤演变趋势

3种工况下河道各段河床比降和断面变化较大，潮流受地形控导作用，挟沙能力在不同工况及同工况不同河段间存在差异，冲淤结果不一(图 7).工程前，河口段受外海潮流动力影响和主槽束水作用，较强的潮流冲刷河槽.而在局部河道断面增大的河道弯段，由于断面涨落潮流速小于其上、下游河段，潮流挟沙能力显著降低，泥沙纷纷落淤，呈现出大面积淤积态势.海相泥沙可被涨潮流带至闸下，但落潮时流速极小且历时较长，致使淤积主要发生在主槽内，闸趾处最大淤积厚度为0.65 m.一期工程后，由于上段河道浚深拓宽，纳潮量增大，受主槽束缚的下游水流冲刷能力和挟沙能力增强，继续冲刷下游滩槽型河段，而上段河道断面浚深拓宽后断面流速减小，挟沙能力骤减，涨潮时大量泥沙淤积在此且微弱的落潮流不能将泥沙带出工程段，造成工程段回淤强烈，年淤积厚度为0.10~0.65 m.整个河口工程河段与非工程河段出现相反的冲淤态势，全河冲淤相抵后净淤积量为3.4×10³m³.二期工程后，口门段仍为强动力区，潮流冲刷河床(最大年冲刷厚度0.6 m)，泥沙随涨潮流进入拓宽段，但该段断面流速减小，挟沙能力降低，部分泥沙在拓宽段落淤.潮流绕过弯段时，受凸岸挑流作用，局部流速增大，冲刷凸岸，冲刷的泥沙被潮流带至闸下，且相对较小的落潮流不能将落淤的泥沙再次扬起带出河道，导致闸趾处淤积厚度达0.6 m.虽然二期工程后，下段河道疏浚，输水输沙量增加，年河道淤积总量为1.2×10⁴ m³，但由于河床坡降均匀一致，水深增加，冲淤分布较均匀.

3.3 开闸放水冲淤效果

对于建闸河口，河道淤积是其必然结果，为了维护整治河道的通航能力，通常采用开闸放水的方式进行主动冲刷，因此开闸放水的冲刷效果也是评价工程优劣的标准.考虑到双龙河口的泄流能力，分别以40%，20%和10% 3种不同频率丰水年设计流量100 m³·s^-1，173 m³·s^-1，271 m³·s^-1为特征流量，对一期工程、二期工程后每年放水7d的年河床冲淤演变进行模拟.从图 8、图 9、图 10可以看出，随着上游开闸流量的增大，河道年淤积面积和淤积厚度逐渐减小，但沿程冲刷趋势2种工况存在明显差异.一期工程后，由于下游河道并未整治，过流断面相对较小，使得整治河段被冲起的泥沙并不能被带出河道，工程段放水冲淤效果不明显.在271 m³·s^-1的流量下闸趾处的冲刷厚度仅为0.30 m，而对于0.40~1.45 km的严重淤积河段，最大冲刷厚度不超过0.12 m，因此工程段仍需采用疏浚措施对其清淤.在未整治的河口段，水流主要集中在深泓处，流速较大且集中，开闸放水后，深泓处流速增大较快，造成该河段冲刷严重.二期工程后，开闸放水冲刷效果在近闸段(距闸趾0~700 m)冲刷程度最大，而在拓宽段(距闸趾2 500~3 300 m)冲刷效果最差.在100 m³·s^-1流量下闸趾处和拓宽段仍存在淤积，但随流量增至173 m³·s^-1和增至271 m³·s^-1，除拓宽段仍存在少量淤积外，其余河段均呈冲刷态势.

从模拟结果可见，闸下河道整治后的冲淤分布与纳潮量增加、局部浚深拓宽密切联系.河道整治导致的纳潮量增加在工程初期会引起中下段的冲刷，而局部河段的拓宽和浚深会使该段产生淤积.一期工程段和二期工程的拓宽段都相对宽阔，产生了淤积且开闸放水减淤效果不明显.水流进入拓宽河道时，过水断面扩大，导致河道流速减小，挟沙能力降低而使泥沙纷纷落淤.若遇坡降不均或水深突然增加(如一期工程)，在河道突然展宽区域很可能形成由水流横向流速梯度引起的回流(如图 6b)，并造成主流流速不对称性加剧，泥沙在回流区与主流交界带严重淤积.虽然在无径流的建闸河口河道淤积不可避免，但为延长工程寿命，在闸下整治时应尽量对整条河道进行整治并使河床坡降均匀一致，避免河道突然展宽.不宜采用河道局部浚深拓宽的整治方案，避免加重水流流速不对称性和回流发生.同时可辅以定期开闸放水和局部清淤等措施维护河道生态环境.

4 结论

基于MIKE 21建立和验证的双龙河口潮流、泥沙输运及河床演变模型，通过对双龙河口整治一期工程和二期工程后的水动力和河床演变进行模拟分析，得到以下结论：

(1) 建闸河口水动力主要受地形的控导作用和纳潮量改变的影响.河道局部整治(一期工程)增加纳潮量并改善非工程段涨落潮不平衡状况；而工程段由于浚深拓宽，涨落潮不对称加剧，局部产生回流区.全河口整治(二期工程)后，河道坡降均匀，纳潮量显著增加，涨落潮流速增加且不对称性减小.

(2) 工程后，水动力改变使河道内冲淤趋势发生调整.上段河道局部整治时，工程段因浚深拓宽流速大大降低，发生强烈淤积，非工程河段因水流冲刷能力增强而发生冲刷.下段河道疏浚后，口门段增大的流速冲刷河床并将泥沙送至上段河道，使拓宽段和闸趾处发生淤积.

(3) 开闸放水冲刷效果随工程方案和流量不同存在差异.上段河道局部整治时，由于其下游过流断面相对较小，使得被冲起的大部分泥沙不能被带出河道，整治段放水冲刷效果并不明显.全河整治后，冲刷范围随开闸放水流量增大而不断下移，但高效冲刷范围主要集中在距闸趾0~700 m处.

(4) 整治工程中相对宽阔的河段易产生淤积且开闸放水冲刷效果较弱.因此，在闸下整治时应尽量对整条河道进行整治并使河床坡降均匀一致，避免河道突然展宽，不宜采用河道局部挖深拓宽，以免加重水流流速不对称性和造成回流，同时可辅以定期开闸放水和局部清淤等措施维护河道水深.