1. 引言
三峡工程在发挥防洪、发电、航运等多方面巨大效益的同时，也对长江流域的生态环境产生了一定的影响，受到社会的广泛关注 [1] [2] 。自2003年，三峡工程开始试运行，采取蓄放交替的方式运行。长江干流库区段的水位在145 m至175 m之间波动，导致了库区一级支流河口的水动力学特性也呈周期性变化 [3] ，在天然河流与河流型水库之间切换，主要表现在：蓄水过程中，由于长江干流水位高于支流，形成回灌，改变了支流的水动力特性，进而影响了水体中沉积物中迁移转化过程 [4] ；在放水过程中，虽然表层流速较大，但底层流速与蓄水期相差不大，有可能导致营养盐、重金属在底层蓄积。这些周期性变化造成了整个支流中物质与能量的迁移规律发生变化。
由于天然河流尺度较大，采用实测研究具有相当难度，近年来，模拟结合实测研究成为了研究者应用最多的手段。目前，众多针对河流、湖泊以及海洋的水动力学模拟软件已越来越多的被应用，如MAPGIS [5] 、MIKE [6] 以及EFDC [7] 等，这些软件都有各自的着重点和应用条件。深入研究一级支流的水动力学特征及演变规律，对于库区支流水环境保护，改善三峡库区水环境状况，保证供水安全具有十分重要作用。本文以三峡库区库尾的较大一级支流御临河为研究对象，结合多年的实测结果，应用EFDC [8] 模型，研究御临河河口段在三峡水库蓄放过程中的水动力学特征，为御临河以及库区其它支流的水环境保护与治理提供参考。
2. 研究对象
御临河发源于四川省大竹县铜锣山，在洛碛太洪岗注入长江，研究区域如图1所示。御临河全长约218.2 km，流域面积3861 km2，重庆市境内流域面积772.8 km2，水域面积52.97 hm2，河口多年平均流量50.72 m3/s。御临河是长江上游北岸的一级支流，也是三峡水库库尾的重要一级支流水系。三峡大坝蓄水后，在库区回水顶托的作用下，御临河形成了约20 km的回水区域。
3. 水动力模型
3.1. 模型概述
EFDC (Environmental Fluid Dynamics Code)于1988年由美国弗吉尼亚海洋科学研究所(IMS)和威廉玛丽学院的海洋科学学院共同开发。该模型是一个开源的地表水建模系统，将水动力、泥沙、污染物和水质模块完全集成在一个单一的源代码内。可以很好的用于包括河流、湖泊、海湾、湿地的一维、二维、三维水动力模拟 [9] [10] 。
3.2. 模型水动力模块的主控方程
为适应实际边界，把正交曲线坐标转换成
坐标；在EFDC模型垂向上采用
坐标变换，能较好的拟合近岸复杂的岸线和地形 [11] ，从而实现在垂直方向上结合重力矢量和水下地形边界，以及允许长波运动的自由液面：
(1)
式中：Z取值范围[0,1]；
代表
转换前的实际垂向物理坐标；总深度
，
为底床高程，
为自由水面高程。
水动力主控方程包括动量方程、连续性方程，其中动量方程采用了湍流运动方程的垂直流体静力学边界层形式的转换以及Boussinesq近似作为变密度结果：
动量方程：
(2)
(3)
连续性方程：
(4)
式中：
和
为曲线正交坐标系下
和
方向的流速分量；
为
坐标下垂向流速；
、
和
为Jacobian曲线正交坐标转换系数，
；
为相对静水压力；
为柯氏力参量；
为垂向紊流黏滞系数；
为垂向紊动扩散系数；
和
为动量源汇项。
为提供垂直湍流粘度和扩散系数，EFDC采用了
阶Mellor-Yamada紊流模型：
(5)
(6)
(7)
式中：
为紊动混合长尺度；
和
为稳定性函数。
4. 网格生成及边界条件
4.1. 网格划分
网格是模型求解的最小单元，为了使自然河道的轮廓显现在软件上，就必须利用矩形网格或正交网格来构造河流的轮廓，使模型的边界能与实际边界吻合。网格划分的质量直接关系到模拟的精度。但是自然水体水岸线的几何形状往往是不规则的，在进行数值模拟时，如果采用矩形网格会受到各种限制，模型的边界不能与实际边界很好的吻合，边界上曲折岸线造成的锯齿效应明显，采用贴体坐标系建立正交曲线网格可以省去大量的交叉导数项，模拟计算的速度也能有一定的提高。Delft3D所建立的网格具有正交性、光滑性(连续性)及疏密度可控性等优点 [12] ，所以本次模拟采用由荷兰Delft Hydraulic公司开发的Delft3D作为贴体网格的生成工具，研究区域水平网格数为1820个，垂向共分5层。如图2所示。建立并完善后的网格文件导入EFDC软件。
4.2. 河底地形的获取与概化
御临河是典型的山地河流，河底起伏比较大，上游和下游落差很大，部分河底为的碎石，在放水期，流速
Figure 2. Boundary and curvilinear orthogonal grid based on Delft3D
图2. 基于Delft3D的河流边界与曲线正交网格
很急。对于这样的河流，河底粗糙率是一个模型构建中的关键因素之一，与河床地形、水深、雷诺数、弗劳德数有关 [13] ，河底粗糙率确定的准确与否直接关系到最终模拟的准确性。但是，由于水下地形的获取具有较大的难度，很多研究在建立模型时常常忽略水下地形的起伏，将河底视为平底或是经验取值 [14] [15] 。本次模拟利用声学多普勒流速测量仪(ADCP)对研究区域设置的10个横向断面及中泓线进行了反复测量，获得河床地形数据，再结合网格文件构造出了御临河河口段的三维地形，最终确定河底粗糙率为0.035。御临河的三维地形概化图，如图3所示。
4.3. 边界条件设置
边界条件是模型模拟的基础设置之一，本次模拟设置的边界条件为流量边界和水位界，流量的入口边界位于
Figure 3. 3-dimentional map of river bed terrain
图3. 三维河底地形概化图
河流上游，出口边界位于河流的下游，南部的水位边界数据从长江海事局获得。模型的最小水深限制为0.1 m，利用干湿网格法模拟河口的流场，消除负水深的出现。气压、气温、风速以及风向等气象参数均来自于当地气象站。
5. 结果分析
5.1. 模型验证
5.1.1. 流速验证
本次模拟，利用2012年~2014年两个周期12次多普勒流速仪的实测流速数据对水动力模型进行了验证，验证结果如图4所示，误差如表1所示。从验证结果可知，模拟值与实测值吻合度较高，平均相对误差分别为9.20%和11.50%。由于精确了河底粗糙率，因此，本次验证误差较文献中的验证误差小 [16] [17] 。
5.1.2. 水位验证
在很多河流、河口、海湾及湿地模拟研究中，水位往往是进行水动力模型验证的首选指标，采用水位进行模型验证的案例远多于流速验证 [18] 。本次模拟的水位验证同样利用了12次实测的水位数据进行，验证结果如图5所示，误差如表1所示。从验证结果可以看出，2013~2014年度，两次的模拟误差比较大，第一次可能是由于降雨的原因，正好在2013年8月中旬有一段时间降雨比较多，河口水位比较高，导致误差比较大；本模型设置河口水位与长江水位等同，但是实际上，在放水末期，河口水位和长江水位存在一定得差值，这可能是导致第二次误差比较大的原因。实测值与模拟值的平均相对误差分别为8.9%和5.0%，本模型具有较高的吻合度。
5.2. 流场分析
5.2.1. 平面二维流场分析
模拟了2012~2014年三峡蓄、放过程御临河河口的流场情况，结果如图6所示。由于御临河处于三峡库区库尾，因此，在蓄水过程中的响应慢于库首河流。在蓄水初期，如图6(a)所示，由于长江干流水位并不高，河口流场主要受御临河上游来水作用，正向流为主导流向，回水现象并不明显，只是在北岸一侧出现部分回流。蓄水中期(图6(b))，此时长江干流水位已超过御临河河口水位，回水已成为河口主要流向。根据流场图显示，此时在河口形成了一个大的环流，整个回流区可以划分为三个区域，其中北岸为主流区，南岸为负回流区，中心为正回流区 [19] 。在蓄水末期(图6(c))，随着长江干流水位持续增至最高水位，由于受到干流水体顶托，流速趋近于0，且方向由中期的负方向变为正方向，河口的流场变得紊乱，这可能由于以下原因导致：一是主流区与
Figure 4. Water velocity verification of estuary
图4. 河口流速验证
回流区水流的剪切作用起到了重要的作用, 它还引起两个区域水体、动量和能量的交换 [20] ，使得回流区向上移动；二是在回水末期，河口水位已经与长江干流水位相同，长江与河口交互影响的区域往上游移动。
在放水过程中，御临河响应快于库首支流，根据模拟结果显示，放水过程，御临河河口流速呈逐渐增大趋势，到放水末期流速到达最大，呈现出自然河流的流态特征。
5.2.2. 纵向流场分析
图7是2013~2014年的纵向流场图，模拟时刻与图6相对应。在蓄水初期，当长江水位达到168~169 m左右，此时，河口表层速率为0.083 m/s，方向为正方向，底部速率为0.07 m/s，方向为负方向；到蓄水中期，河口水位变为173~174米，回水区域变大，表层速率为0.128 m/s，方向为正，底层速率为0.122 m/s，方向为负，这说明回水以下切的方式流向御临河。蓄水中期的速率比蓄水初期的速率大，这是因为在蓄水中期长江水位与御临河水位差变大，从而静水压差变大，导致河口的速率变大。在蓄水末期河口水位变为175 m，流速方向开始变的紊乱，表层流速和底层流速的方向开始逐渐由负方向变为正方向，表层速率变为0.076 m/s，底层速率变为0.064 m/s。
Table 1. Analysis of relative error about velocity and water level
表1. 流速及水位相对误差分析
(a) (b)
Figure 8. Comparison of surface and bottom velocity in the estuary
图8. 河口表底层流速对比分析
并且从图7(a)、图7(b)可以看出，在蓄水期，河口处在纵向上同样存在漩涡，平面上的竖轴环流与立面的径向环流的复杂叠加，主、回流交界面存在较大的流速梯度，而压强梯度是产生流速梯度的主要原因之一，同时造成水体之间产生紊动掺混层，紊动掺混层的存在加剧了主、回流之间的动量、能量和质量交换 [18] 。
在放水初期，河口水位降至174 m，表层速率为0.028 m/s，底层速率为0.033 m/s，方向都变为正方向，速率变到最小。放水中期，水位继续降至168 m，此时，御临河上游来水成为河口流速的主要驱动力，表层速率和底层速率开始变大，分别为0.062 m/s和0.051 m/s。放水末期，表层流速和底层流速的方向，都是变为正向，这时河口流速继续增大，已经完全呈现出自然河流的流态特征，表层最大速率为0.1 m/s，底层最大流速为0.08 m/s。
5.2.3. 表底层流场分析
图8是御临河河口段的表层和底层流速随时间的变化。在图8(a)中，2012~2013年度的开始模拟时间是11月10号，是在蓄水的中后期，这其中的部分时间，表层流速为正，底层流速为负。2013~2014年度，从2013年8月10号开始模拟，从图8(b)中看出，御临河在9~10月左右刚刚开始形成回水，在这短时间内，部分时间的表层流速是正的，底层流速是负的，正好印证了前面所述；到蓄水的后期，底层流速和表层流速变为正的，蓄水完成。
在2012~2013年度的蓄水期，河口表层和底层的平均速率各为0.026 m/s和0.027 m/s；放水期，表层和底层的平均速率各为0.076 m/s和0.061 m/s。在蓄水期表层速率与底层速率差值与表层速率比值为3.85%，在放水期，这个比值为19.74%；这说明在放水期，表层速率和底层速率的差值更大。同样，从图8中看出，从每年的1~2月左右开始，表层和底层速率差值开始变大，正好是三峡大坝开始放水的时间段。
6. 结论
御临河代表了位于三峡库区库尾的一级支流。根据流场实测及模拟结果显示，受到三峡工程蓄放过程的影响，御临河出现回水与自然流交替运行的现象，其中，蓄水过程受长江干流水位的影响呈现较为复杂的流场变化规律；而放水过程，则是一个随着干流水位降低，流速由小逐渐增大最后变为自然河流的过程。并且，在蓄放过程中，受到干流下切式回流的影响，在河口会出现横向和纵向环流的现象，改变了河口物质与能量的迁移途径与转化规律。
基金项目
国家自然科学基金面上项目(51478061)。