三峡工程建成后在长江上游形成了河道型水库，将改变库区和支流的水动力条件.随着区域社会经济的发展，大量的点源和非点源造成的污染物排入河流水体中，在多种因素作用下引起了库区水质的下降.近年来，三峡库区支流(如香溪)多次爆发了水华事件(李媛等，2012)，给流域的生态环境和当地居民的生活造成了不利影响.

河流在复杂边界驱动力的作用下(如水库调度和海洋潮流)，河流水体中的污染物输移也变得较为复杂.基于拉格朗日法的粒子跟踪模型可以模拟单个粒子或粒子群的运动轨迹，因此，被广泛地应用于研究和模拟自然界中颗粒物体在水体中的输运轨迹、掺混和交换等问题.例如，Andre(2008)采用粒子跟踪模型模拟了水体中浮游生物的浮游运动，发现粒子模拟不仅可以了解浮游生物个体的微观尺度的运动，也可以了解宏观尺度的种群分布状态，模型的关键是对浮游生物行为特征的物理描述.Chen等(2010)采用三维粒子跟踪模型研究了淡水河口内污染物的运动情况，发现潮流驱动和风生流动是影响河口污染物运动的主要原因.Liu等(2011)应用粒子跟踪模型研究发现，淡水河口的水体密度差异引起的立面环流使粒子不能沉降至河床底部.Ulf(2011)总结比较了若干种粒子跟踪算法的模拟精度，指出粒子的对流为流场驱动的确定性运动，粒子的扩散由水流的随机紊动扩散引起，并基于此对粒子跟踪算法做了改进，取得了精度较好的模拟结果.可见，粒子跟踪模型有着广阔的应用前景.而传统的基于欧拉场的模拟水藻中叶绿素浓度的生态模型不能模拟单个细胞颗粒的运动轨迹，也难于了解某一区域水藻细胞颗粒的生长和输移运动情况.

因此，本文将应用海洋动力学模型ELcirc(Zhang et al., 2004)对2007年秋季三峡水库蓄水期内香溪河的水动力情况进行模拟，在此基础上应用考虑营养物质浓度、水下光照、水温和流速影响水藻细胞生长的非守恒三维粒子跟踪模型对库湾内的水流及水藻颗粒输移进行模拟，以期对三峡库区水质的治理提供科学依据.

ELcirc模型基于在欧拉场的非结构网格离散控制方程，可求解自由水面波动方程、流场.采用基于Boussinesq假定的静水压力方程时，ELcirc模型为准三维模型，即不求解垂向运动方程，通过控制体内的水量平衡及求解得到的水平向流速分量，即可得到垂向流速，这可明显减少计算量，并可应用于复杂长河段(如香溪河)进行三维场的模拟(李健，2012).控制方程包括质量守恒方程、动量守恒方程和物质输移方程：

式中，(x，y)为平面笛卡尔坐标(m)；z为垂向坐标(m)，向上为正；t为时间(s)；Z_ini为计算初始水位(m)；η为自由水位波动(m)；h为水深(m)；u、v、w分别为法向、切向和垂向流速(m · s^-1)；g为重力加速度(m · s^-2)；f为柯氏力系数；K_mh和K_mv分别为动量方程中的水平和垂向紊动粘性系数(m² · s^-1)；K_h和K_v分别为物质浓度方程中的水平和垂向扩散系数(m² · s^-1)；C_i为营养物质(包括总磷和总氮)的浓度(mg · L^-1)，本文模型不考虑短时间内营养物质的沉降及生化反应的源项.

粒子运动轨迹的模拟是在基于ELcirc的水动力场计算结果下进行的(Zhang et al., 2004)，模拟粒子群在非结构网格内的运动轨迹，粒子在每一时间步内的运动由确定性的对流和随机性的紊动扩散两部分组成，具体由下式计算：

式中，(X，Y，Z)是粒子在t时刻和t+Δt时刻的三维空间位置(m)，(u，v，w)是水动力模型计算得到的三维流速(m · s^-1)，K_h、K_v分别为物质浓度的水平和垂向扩散系数(m² · s^-1)，(R_x，R_y，R_z)为均匀分布的随机数，范围在-1~1之间，粒子从扩散度高的位置向扩散度低的位置移动.

为了模拟水深方向上受到营养物质浓度、水下光照、水温和水动力限制作用影响下的水藻生长过程，本文在守恒性三维粒子模型的基础上，考虑了以上4种因素影响下的粒子增加过程.水藻细胞颗粒的增殖过程的计算公式为：

式中，p_n是水藻细胞颗粒在t时刻和t+Δt时刻的数目(个)，α为水藻细胞颗粒的增值率(个 · m^-3 · h^-1)，f_N、f_I、f_T和f_U分别为营养物质浓度、水下光照、水温和流速对水藻细胞颗粒增殖的限制因子，数值均在0~1之间，计算公式分别如下(Stefan et al., 1983; 李锦秀等，2005; 李健，2012)：

式中，C_TP和C_TN分别为总磷和总氮的浓度(mg · L^-1)，K_mP和K_mN分别为水藻细胞吸收总磷和总氮的半饱和系数，取值分别为0.01和0.003(李健，2012)，I₀为每日水面的日照强度(lx · d^-1)，I_m为浮游植物的饱和光强(lx · d^-1)，Z_d为水面至目标水层的垂直距离(m)，Δz为计算单元的厚度(m)，T_m为蓝藻生长最适宜温度，取值20℃，K_Tg为温度对蓝藻生长的限制系数(当T>T_m时K_Tg=0.006，当T<T_m时K_Tg=0.008)，K_e为水下消光系数，U为网格节点计算流速(m · s^-1).

在平面网格单元中进行粒子轨迹跟踪计算，算法示意图如图 1所示.粒子的新时刻的位置坐标(x_new，y_new)采用下列公式计算(Martin et al., 2005)：

图 1(Fig. 1)

图 1 粒子轨迹跟踪算法示意图 Fig. 1 Schematic map of particles tracking algorithm

式中，(x_new，y_new)、(x_old，y_old)分别为新、旧时刻的粒子位置坐标(m)，Δx和Δy分别为控制单元x方向和y方向的边长(m)，U_x1、U_x2和V_y1、V_y2分别为进出控制单元的x和y方向流速(m · s^-1)，x₁、x₂和y₁、y₂分别是控制单元节点的x和y坐标(m)，τ_e为每一时间步内粒子轨迹跟踪计算的子时间步长(s).

以上计算式中粒子运动速度由周围网格节点上的计算流速(水动力学模型计算结果)得到.由于水动力模型中采用的时间步长较大，粒子跟踪模型需将其分解为多个子步(τ_e)，以满足粒子下一时刻的位置处于计算范围的条件.另外，t+Δt时刻新增加的粒子空间位置，采用线性插值的方法添加到在新旧时刻的粒子的中点处，具体如图 1中的P_add；然后添加新粒子后的粒子群作为下一时刻粒子轨迹计算的初始条件.可见，本文的粒子轨迹跟踪模型是基于欧拉背景网格的跟踪计算，计算效率可以满足天然河流中的物质输移模拟研究.

采用弯道水槽试验资料对水动力模型进行验证，结果如图 2所示.弯道横断面为矩形，弯道半径为8.53 m，两弯道连接过渡段为4.27 m的直水槽，弯道的进出口由长为2.13 m的直段过渡，水槽宽2.34 m.实验进口流量Q=0.09486 m³ · s^-1，进口平均流速为0.366 m · s^-1，平均水深为0.115 m，水槽比降为J=0.35‰，糙率为0.015.

图 2(Fig. 2)

图 2 弯道示意图(单位：m) Fig. 2 Sketch map of the flume

如图 2所示，实验中在水槽的第2个弯道处进行了若干断面的纵向和横向流速测量，测量断面号分别为：S0、S2、1、PI/8、2、3、3PI/8和4，每个断面布置2个测点，限于篇幅，本文仅给出PI/8断面的流速验证结果(图 3).可知，在水槽中部纵向和横向的流速模拟计算值与实测值符合良好.

图 3(Fig. 3)

图 3 PI/8断面测点流速验证(a.纵向流速U，b.横向流速V) Fig. 3 Velocity verification at PI/8 section(a.longitudinal velocity U，b.transverse velocity V)

在水槽进口处释放5颗粒子，粒子在剖面上的初始位置如图 4左侧所示，粒子将随水流运动，水槽粒子运动模拟不考虑粒子的增加过程.如图 4a所示，当不考虑水流紊动的随机扩散作用时，粒子的运动完全取决于欧拉流场的驱动，运动轨迹是规则光滑的曲线.由图 4b可知，在考虑紊动扩散作用后，粒子的运动轨迹不再光滑，可以反映出弯道水流输移的特点，在弯道水体表面附近横向上粒子有从凸岸游移至凹岸的趋势，并且只要粒子处于水体当中，粒子即可从进口运动至出口(满足数量守恒).可见，考虑对流扩散作用的拉格朗日粒子跟踪模型可以模拟出弯道水槽污染物输移的特性，粒子跟踪模型可以模拟自然河流中的物质输移.

图 4(Fig. 4)

图 4 粒子运动轨迹(a.不考虑随机扩散，b.考虑随机扩散) Fig. 4 Particles tracking modeling(a. without r and om dispersion，b. with r and om dispersion)

香溪河位于湖北宜昌市境内，流经兴山县和秭归县，是长江三峡水库湖北省库区段内的第一大支流.流域面积2939 km²，控制水文站为兴山水文站，年平均流量47.4 m³ · s^-1，香溪支流高岚河的控制水文站为建阳坪水文站.本文将应用上文经过验证的ELcirc模型和粒子跟踪模型，研究2007年9—10月三峡水库蓄水期内香溪河的水动力和水藻颗粒的运动情况.

计算区域包括香溪河主河道和支流高岚河，各断面位置见图 5.从白沙河和南阳河交汇处至香溪河口，水平向划分19512个四边形网格单元，垂向上分为30层.进口边界条件采用对应时段的兴山水文站和建阳坪水文站的流量，出口边界采用香溪河口的日平均水位(1985黄海高程系).如图 6所示，三峡水库蓄水时长江干流水体将大量倒灌入香溪河， 而三峡水库泄水时香溪河口表层水流向下流动，但底部仍有水体缓慢倒灌入香溪河，水流条件较为复杂.表层水流流速明显大于底部，但均小于0.01 m · s^-1，处于准静止状态，水体的复杂流动将影响污染物及水藻颗粒的输移.

图 5(Fig. 5)

图 5 香溪河示意图 Fig. 5 Sketch map of Xiangxi River

图 6(Fig. 6)

图 6 香溪河口局部计算流场矢量图 Fig. 6 The flow field at the river mouth

本文中水藻细胞颗粒的轨迹跟踪计算是在欧拉框架下进行的，模型模拟了水动力场和营养物质(总磷和总氮)浓度场.由图 7可知，在峡口镇和刘草坡2个测点的总磷(TP)的计算值与实测值的变化过程吻合良好，而总氮(TN)的模拟效果欠佳.峡口镇的总磷浓度在9月30日达到峰值0.14 mg · L^-1，刘草坡的总磷浓度有2个峰值，分别是9月28日的0.16 mg · L^-1和10月3日的0.12 mg · L^-1，上游总磷浓度大于下游的总磷浓度.这是因为香溪河上游有大量磷矿开采企业和磷化肥厂，实测值具有明显的规律;而总氮受到水藻吸收空气中的氮元素及面源污染入汇的影响，实测过程无明显规律，而数值模拟过程为连续光滑过程，因此，模拟效果欠佳，但能反映出总氮的大致变化过程.

图 7(Fig. 7)

图 7 营养物质浓度计算值与实测值对比 Fig. 7 Measured and simulated nutrients′ concentrations

野外测量中一般采用叶绿素浓度来表征水体中藻类的生物量状态，2007年2—5月峡口镇的实测叶绿素a浓度(杨正健等，2008)结果表明：香溪库湾水体中的浮游藻类受光照、水温等因素的影响，叶绿素浓度在水面处较高，沿垂向逐渐衰减，至10 m水深处已很低，且表层叶绿素浓度增殖较快，具体如图 8所示.

图 8(Fig. 8)

图 8 峡口镇实测叶绿素浓度垂向分布 Fig. 8 Measured vertical distribution of chlorophyll at Xiakou Station

下面将应用考虑粒子数目增加的三维粒子轨迹跟踪模型对香溪河在2007年9—10月发生的水华爆发过程进行模拟.从2007年9月25日开始，在秭归至贾家店河段释放粒子(图 9a)，共分为10层粒子进行计算.为模拟浓度的垂线分布，按图 8垂向上各测点实测叶绿素浓度占沿水深积分的叶绿素总量的百分比布置各层的粒子数目，初始时刻(2007-9-25)在自由水面处释放656200颗粒子，沿水深方向依次递减，至9 m水深处的粒子为9580颗.水藻细胞颗粒的增值率α根据室内蓝藻培养皿实验研究得到，取值2.0×10³ 个 · m^-3 · h^-1.

图 9(Fig. 9)

图 9 香溪库湾内的粒子运动(a.粒子群运动(表层)，b.单个粒子运动轨迹) Fig. 9 Particles′ movement in Xiangxi River(a.particles swarming at water surface，b. single cell′s moving path)

三峡水库蓄水时长江干流水体将大量倒灌入香溪河，水流条件较为复杂，库湾下游水体流速普遍小于 0.01 m · s^-1，水体的复杂流动将影响污染物及水藻颗粒的输移.由图 9a可知，在模拟时段内粒子群向上游运动，2007-9-26至2007-10-8在最上游的粒子仅从贾家店运动到峡口镇附近，运动距离约600 m，边岸处由于流速较河道中间位置小，因此，边岸处的粒子有滞留现象，如2007-10-8的粒子位置图中秭归至三闾的河段同时有部分颗粒进入支流高岚河道内.单个粒子的运动轨迹(图 9b)表明：细胞颗粒受复杂水动力场的搬运作用，形成随机游走的运动轨迹，且运动距离不远(约50 m)，说明粒子模型的计算可以反映出粒子的物理运动图景.

河流中浮游藻类的浓度分布受水流条件、营养物质浓度、水温和水下光照等因素的综合影响，尤其在水深方向上光照强度由于水体背景反射及水体中的水藻和泥沙颗粒的遮蔽，光照很快衰减，这对水藻垂向上的分布及生长影响明显.对考虑粒子数目增加的三维粒子轨迹跟踪模型各层粒子数目经过若干天的计算值统计(图 10)发现，粒子模型可以模拟出水藻颗粒沿水深方向的分布趋势及表层水藻颗粒的急剧增殖过程，非守恒三维粒子模型模拟的物理现象与实测到的叶绿素浓度垂向变化(图 8)相似.

图 10(Fig. 10)

图 10 各层水藻颗粒增殖过程 Fig. 10 The vertical algal cells′ reproduction process

1)本文首先采用室内弯道水槽试验资料对海洋模型ELcirc的水动力模块和粒子轨迹跟踪模型进行了初步验证，发现计算值与实测值符合良好.然后应用验证后的数学模型对2007年9—10月三峡水库蓄水期内香溪河的水动力和水藻颗粒的运动进行了数值模拟，结果发现，ELcirc模型可以精确地模拟弯道水流特性和地形及边界条件变化剧烈的自然河流的水动力场，可较好地模拟三峡水库蓄水期间长江干流水体倒灌入香溪库湾而产生回流的过程.

2)三维粒子轨迹跟踪模型可以模拟出弯道水流中物质的输移特性，物质输移由确定性的对流运动和随机性的紊动扩散两部分构成，粒子在弯道处由凸岸向凹岸偏转.

3)三峡水库蓄水期内香溪河道中的浮游藻类颗粒在回流作用下，有向上游输移的现象，并且考虑营养物质浓度、水温、流速和水下光照影响下的水藻颗粒增殖的非守恒粒子跟踪模型可以较好地模拟藻类颗粒在三维空间上的输移和紊动扩散，以及水藻急剧生长下的水华过程.