附着性藻类生长-迁移模型及其在南水北调中线工程中的应用

胡豫英， 胡 鹏, 朱 杰， 廖卫红

(1.浙江大学 海洋学院， 浙江 舟山 316021； 2.珠江水利委员会 珠江水利科学研究院， 广东广州 510611；3.北京工业大学， 北京 100124； 4.中国水利水电科学研究院， 北京 100038)

1 研究背景

附着性藻类是指附着于水体基质表面生长的藻类，由于其环境适应能力强，可通过孢子/配子等多种途径繁殖，广泛分布于全球淡水及浅海区域，常见的有刚毛藻、浒苔等。一般而言，附着性藻类的适度生长能够吸收水体中的氮、磷等营养盐，净化水质，可以被用于水体富营养化的治理[1]，某些种类甚至可以食用(如浒苔)。但附着性藻类异常增殖后其冲刷脱落物往往会堵塞取水管道和水利设施，影响其正常功能；其次，脱落物分解腐烂后往往漂浮于水面，产生难闻气味并污染水质[2]。南水北调中线工程运行中的实际监测资料显示，附着性藻类往往会在春、秋季大量增殖，附着面积甚至能达到干渠面积的43%，且在冲刷脱落后，脱落的残体会随流输移进而堵塞分水口门，影响干渠的正常输水[3]。虽然目前总干渠沿线分水闸均安装了拦藻除藻装置[3-5]，但当附着性藻类爆发期到来时，单个除藻装置的清洗频率达到5次/d，清运的藻类多达3 t，消耗了大量的人力物力[3]，因而对附着性藻类增殖迁移过程的进一步研究迫在眉睫[6]。

随着计算机及计算技术的不断发展，水动力-水质模型广泛应用于各种水环境问题的分析研究中，成为解释自然水体中物质迁移变化规律的重要科学工具。通用的水动力-水质模型(如the environmental fluid dynamics code, EFDC和water quality analysis simulation program, WASP等)将藻类视为水中溶质的一部分，通过考虑营养盐、光照等因素对悬浮态藻类生长的影响来模拟生态系统中物质能量循环和传递过程[7]。然而附着性藻类的生长除了受上述因素影响外，水动力条件可能是一个更重要的影响因素[8]。水动力的影响主要体现在两个方面：首先水动力过程可通过影响藻光合作用效率、形态发展等影响附着性藻类的生长密度，南水北调中线工程中的实际监测也显示，附着性藻类更容易聚集于分水口和弯道等水动力条件较弱的部位[3]；其次，水体也是附着性藻类完成繁殖过程的重要媒介，成熟藻体释放的孢子会随着水体输移至合适的部位并附着，以完成藻类附着区域的扩展[9-10]，相关研究表明，藻类孢子的大量释放，是藻类爆发的主要原因和物质基础[11]。通用水质模型往往以悬浮态藻类为研究对象且没有考虑水动力作用对藻类生长的影响，所以无法用于附着性藻类的模拟；现有针对附着性藻类的模型虽然考虑了藻类附着的特性且量化了水流对藻体的冲刷作用，但忽略了附着性藻类的繁殖迁移过程，往往导致模型无法有效描述附着性藻类的迁移扩散现象[12]。

本文通过引入对附着性藻类的孢子释放-迁移-附着过程的数学描述，使模型能够正确体现出附着性藻类的生长迁移过程，弥补现有模型体系的不足，并将该模型用于南水北调中线工程总干渠的藻类增殖模拟中，以探究附着性藻类的增殖规律，为实际工程中的藻类防治提供科学参考。

2 附着藻类生长-迁移模型

2.1 水动力方程

本文的水动力控制方程组由水流的质量方程和动量方程组成，表达式如下：

(1)

(2)

(3)

式中：x、y为笛卡尔坐标系下的空间坐标，m；u、v分别为x、y方向上各自沿深度方向上的积分平均流速，m/s；t为时间，s；h为水深，m；g为重力加速度，m/s2；Sfx=n2u(u2+v2)1/2/h4/3、Sfy=n2v(u2+v2)1/2/h4/3分别为水流在x、y方向的水力坡度；n为底床的曼宁糙率系数；Sbx、Sby分别为x、y方向的底坡源项。

2.2 基于孢子/配子的藻类生长-迁移模型

现有针对附着性藻类的生长模型普遍将区域内藻类密度的变化分为藻类自身的生长与受水流的冲刷作用两部分，经过参数修正，模型能够很好地描述藻类密度在光照、温度和营养盐等条件下的动态变化过程[12]。但对于附着性藻类的迁移扩散过程，现有模型无法很好地表现。本小节基于Fovet等[12]提出的充分考虑了水流冲刷作用的附着藻类生长模型，通过引入对藻类孢子/配子释放-迁移-附着过程的描述及概化，为现有附着性藻类生长模型引入迁移模块，补全模型体系。

附着性藻类作为低等生物，可以依靠孢子、配子以及体细胞营养再生等多变的繁殖方式和大量的后代延续种类进行繁殖，其中孢子/配子繁殖是主要的繁殖方式。当藻体成熟时，位于顶端的藻体细胞会逐渐转化为生殖细胞，之后生殖细胞的细胞囊壁破裂，孢子/配子被释放至环境水体。进入水体后，孢子/配子将主要随着水流运动，寻找合适的基质表面并附着、萌发，形成新的藻类附着区。基于以上对附着性藻类增殖过程的理解，本文将附着性藻类概化为两种不同的形态：附着在基质表面的附着态藻类B和藻类释放的悬浮态的孢子/配子S，这两种形态藻类的相关关系如图1所示。

图1 附着性藻类的两种概化形态及相关关系图解

综上所述，总干渠中藻类的两种形态之间存在着相互耦合的关系，基于对它们迁移变化过程的理解，附着态藻类的密度变化方程和悬浮态孢子/配子的浓度变化方程可分别表示为方程(4)和方程(5)：

(4)

(5)

式中：B为附着性藻类生长密度，cells/m2；S为藻类孢子/配子浓度，cells/m3；G为附着性藻类的生长量， cells/(m2·d)；E为附着性藻类受水流作用的冲刷脱落量，cells/(m2·d)；D为悬浮态孢子沉降量，cells/(m2·d)；rsurvival为孢子/配子沉降后的萌发率，在实验室条件下，孢子/配子的萌发率可以达到96%[11]，但在自然水体中，由于水蚤等生物作用，萌发率会锐减至4 %～6 %[13]；R为孢子/配子的释放量，cells/(m2·d)。

指数形式的藻类生长过程被众多水质模型广泛采用[12]，藻类生长量G可以表示为藻类实际生长速率μ与当前藻类生物量的乘积，计算式如下：

(6)

式中：μ为藻类的实际生长速率，d-1。人工渠道(如南水北调中线工程)中由于水环境简单，当气温适宜时，光照、营养盐等因素往往也能满足藻类生长需要，这时藻类自身密度将成为限制藻类生长的重要因素[12]，故本文将藻类实际生长速率表示为生长速率μ0与藻类密度限制因子的乘积，Bmax为藻类的最大生长密度，cells/m2。

水流对附着性藻类的冲刷作用E是影响藻类密度的另一重要机制，不同模型对该项的考虑有较大的区别。Fovet等[12]通过总结对比分析现有模型对冲刷作用的表达，指出使用水体的切应力是表达冲刷作用的合适参数，其表达式为：

(7)

式中：cd为冲刷效应系数，d-1；τ=ρgu2n2/h2/3为水流切应力，N/m2;τcr为临界切应力，本文取值为0.087 N/m2，只有当水流的切应力大于临界切应力时，附着的藻类才会受到冲刷作用；Br为由于底床粗糙而难以被冲走的藻类生物量，β为无量纲的经验系数，一般取值1～1.5。

沉降作用D作为一个从水体中去除物质的重要物理机制，早已被众多水质模型、水沙模型纳入至模型体系中。本模型中，孢子沉降表现为水体中孢子/配子转移到底床，成为新的附着性藻类，沉降通量可以表示为：

D=ωS

(8)

式中：ω为孢子/配子的沉降速度，取值范围为0.01～0.1 mm/s[14]。

孢子/配子释放作用R表征藻类细胞向水体中释放繁殖体的速率，Lotze等[13]发现附着性藻类可以源源不断地放散孢子/配子，且释放速率可以高达1.2×109cells/(m2·h)；在自然环境下30%的藻体细胞可以转化为生殖细胞并释放孢子[11]，而在实验室条件下87%甚至全部藻体细胞均可以转化为藻体细胞[15]，每个生殖细胞可产生4～32个孢子/配子，从开始转化至孢子/配子释放完毕大约需要3～4 d，故孢子/配子释放作用R可以表示为：

R=αB·rrelease

(9)

式中：α为藻体中生殖细胞占比，本文取值为0.3；rrelease为生殖细胞孢子/配子平均释放速率，cells/d，根据上文数据推算，约为1～8 cells/d。需要指出的是，虽然孢子/配子的释放会损失生殖细胞，进而减小附着藻类的生长密度，但该质量损失效应可直接体现在藻类的综合生长速率μ中，故本文没有在方程式(4)中增加由于孢子/配子释放造成的密度损失表达式。

附着性藻类增殖模型的相关参数总结于表1。

表1 附着性藻类增殖模型参数及取值一览表

2.3 数值算法

利用有限体积法将控制方程组离散，为便于书写将控制方程(1)～(5)写成如下向量形式：

(10)

其中：

(11)

(12)

(13)

(14)

将控制方程组对任意网格单元i进行积分，得：

(15)

式中：E=(F,G)；Ai为计算单元面积，下标i为单元编号，i=1,2,...,Nc，Nc为单元总数。首先，对方程(15)左侧第2项及右侧第1项运用高斯散度定理，将面积分转化为线积分来对流通量进行计算；其次，假设源项在计算单元内均匀分布，以此简化对方程(15)右侧第2项的处理。则方程(15)变为：

(16)

式中：上标n和n+1为时间层级的编号，时间层级间相差一个时间步长Δt；Ei,j为第i个计算网格下第j条边上的对流通量，j=1,2,3。ni,j=(nx,ny)为单元i第j条边上的法向量；Li,j为对应边的边长。本文使用HLLC(Harten-Lay-van Leer-contact)算子[16]计算对流通量Ei,j; 底坡源项(ghSbx)i、(ghSby)i采用Slope flux method求得[17]，如方程(17)所示:

(17)

本文采用非结构三角形网格离散计算区域，如图2所示。图2(a)表示单个网格单元对应了3个单元角点、3条单元边界和3个相邻单元，且这些元素均呈逆时针排序；图2(b)表示共有单元角点的单元拓扑方式，再结合反向距离加权等插值方式，可将存储在网格中心的变量值插值至网格角点；图2(c)表明了单元边界与其两侧网格单元的拓扑关系，模型网格间通量交换利用边界两侧的黎曼状态量，通过HLLC近似黎曼算子计算。

图2 非结构三角形网格的拓扑结构

2.4 模型验证

为探究水流的冲刷作用对附着藻类生长密度的影响，Fovet等[12]在室外搭建了3条8 m×0.22 m (长×宽)、粗糙程度相同的灌溉水渠，分别通入18.9、6.7和5.1 L/s 3种流量，通过调节尾堰使水渠中的水位保持在0.24 m，形成相应流速分别为0.358、0.127和0.097 m/s(水流切应力分别为1.19、0.15和0.087 N/m2)的近似均匀流，并设置了初始密度为1×1010cells/m2的附着性藻类(羽纹脆杆藻，fragilariapinnata)，培养30 d，观察藻类密度的变化情况。

本文将采用表1中所列参数及取值，模拟上述3条水渠中的藻类生长状况，并与Fovet等[12]采用的模型进行结果对比，如图3所示。由图3可看出，当水渠流量为18.9 L/s时(图3(a))，水体流速最大，切应力也最大，藻体受到强烈的冲刷作用而无法维持原有密度，故藻类密度逐渐减小，但由于渠壁的粗糙性，使部分藻类难以被水流冲走，因而藻类密度最终保持为残留密度Br；流量减小至6.7 L/s时(图3(b))，附着藻类虽然也受水流的冲刷作用，但藻类自身的生长量大于被水流的冲刷量，藻类密度逐渐增大；而当流量为5.1 L/s时，水流切应力等于藻类冲刷的临界切应力τcr，藻体将完全不受水流的冲刷作用而始终增长。图3中的比较结果表明，本文模型较好地模拟出了3种流量情况下附着藻类生长密度的变化趋势。

图3 3种不同流量下水渠中附着藻类生长密度变化曲线对比

本文模型由于增加了基于孢子/配子藻类的迁移机制，孢子/配子的沿程释放将使水中孢子/配子的浓度沿程增大，进而使孢子/配子的沉降量也沿程增加，最后渠尾的附着藻类密度大于渠首，如图4所示。虽然这一机制在当前算例中并没有提升模型精度，但对解释附着藻类的迁移扩散现象至关重要。

图4 模拟2 d后水渠中藻类密度、孢子/配子浓度沿程分布(流量为6.7 L/s)

3 附着性藻类增殖模型在南水北调中线工程中的应用

3.1 水动力验证

算例选取南水北调中线工程中的陶岔至刁河站节制闸段进行水动力验证，陶岔位于丹江口水库东岸，也是南水北调中线工程的渠首所在；刁河节制闸位于陶岔下游约15 km处，其间在陶岔下游约4 km处还设置一个肖楼分水口，固定分水流量为12～15 m3/s。计算区域用10 365个三角形计算单元、6 159个单元角点进行离散，网格平均边长约15 m，计算区域位置及网格示意图如图5所示。在模型建立中，干渠地形数据采用实测断面数据；糙率根据前人研究设置为0.014[18]；入口处采用陶岔站2018年4月实测流量过程；出口处设置为刁河站实测水位过程，同时在肖楼分水口处固定设置分水流量为15 m3/s，进、出口边界条件即流量、水位过程线如图6所示。

图5 模型干渠水动力验证区域地形及网格示意图

图6 模型干渠入口、出口边界设定的流量、水位过程线(2018年4月实测)

将刁河站节制闸的出口流量模拟值与实测值进行比较，如图7所示。由图7可以看出，本文模型较好地模拟出了干渠内水体的流动过程，由于渠道输水存在一定量的水流损失，一般情况下占总水量的5%[19]，故后期流量的模拟值较实测值偏高约5%，属于正常误差，因此模型可以进一步用于该区段的水质分析。

图7 刁河站节制闸出口流量模拟值与实测值对比图(2018年4月)

3.2 南水北调中线工程干渠附着藻类增殖的初步模拟

本小节将在表1的参数设定和3.1小节的计算网格条件下，探究干渠边壁附着藻类的生长扩散情况。入口边界给定3种特征流量180、250和350 m3/s，出口水位设置为固定值146.75 m，初始设定水渠中没有附着藻类及孢子/配子，在计算稳定后设置入口孢子/配子浓度为4×107cells/m3，模拟30 d，观察附着藻类在3种特征流量下的生长迁移情况。

图8为3种特征流量下干渠内藻类平均密度随模拟时长的变化曲线。如图8所示，含孢子/配子的水流进入干渠后，孢子/配子会逐渐沉降并附着萌发，成为边壁上的附着藻类，同时由于此时藻类密度较小，渠壁粗糙的表面给附着下来的藻类提供了足够的空间并免受水流的冲刷作用，藻类密度迅速增大，干渠流量对初期藻类生长的影响很小；约7 d后，附着藻类密度达到可冲刷的临界值Br，此时受不同流量水流的冲刷作用，藻类的密度变化曲线开始分化，大流量下附着藻类密度明显低于小流量的藻类密度。模拟30 d后，350 m3/s流量下干渠藻类平均密度为9.57×109cells/m2，仅为180 m3/s流量下干渠藻类平均密度的20%，这一结果与南水北调中线工程中的实际情况相符。

图8 3种特征流量下干渠内藻类平均密度随模拟时长的变化曲线

3种特征流量下模拟30 d后干渠中附着藻类密度及水流流速空间分布分别见图9～11。如图9～11所示，由于渠道中局部水动力条件的差异，附着藻类生长密度的空间分布也十分不均。定性上，弯道顶部由于水流的顶冲作用，藻类密度要明显小于弯道内侧。而当水流流动稳定时，渠道中心水流流速要高于渠道两侧流速，导致附着藻类主要分布于渠道两岸；定量上，当水流流速小于0.4 m/s时，附着藻类便会呈明显淤积态势[3]。

图9 180 m3/s流量下模拟30 d后干渠中附着藻类密度及水流流速空间分布

图10 250 m3/s流量下模拟30 d后干渠中附着藻类密度及水流流速空间分布

图11 350 m3/s流量下模拟30 d后干渠中附着藻类密度和水流流速空间分布

4 讨 论

近年来，附着性藻类的增殖迁移问题随着堵塞输水管线[3]、绿潮[20-21]等现象的频发而被关注，也有众多研究者尝试通过数学模型研究附着藻类受营养盐、光照、水动力条件等影响的生长过程[22]，并且已经能够准确模拟出区域内附着藻类的生长密度等关键指标[22]。但对绿潮等藻类爆发的溯源研究表明，藻类爆发时漂浮藻体以藻类孢子/配子等微观繁殖体为主[21]，是藻类繁殖迁移过程的体现。现有附着藻类模型仅针对固定区域的藻类成熟体的生物量进行模拟，而忽略了附着藻类随水流的繁殖迁移过程，这一方面导致模型无法解释附着藻类生长范围的扩散性，这一点与藻类孢子/配子等微观繁殖体随水流的输移沉降密切相关[10]；另一方面也导致这些模型无法适用于与附着藻类繁殖迁移密切相关的藻类爆发等现象的进一步研究。

将藻类孢子/配子的释放迁移过程引入现有附着藻类生长模型能够补全对附着藻类全生命周期的数学描述，拓展模型的应用场景。同时，得益于附着藻类生长发育过程相关实验的逐渐完善[11,13-15,23,24]，本文模型中与孢子/配子释放沉降的相关经验参数才能够合理的取值，使本文针对南水北调中线干渠中附着藻类的模拟结果符合实际分布规律[3]。但同时也必须指出，目前模型考虑的还不够全面，如附着藻类受营养盐、温度、光照等影响而导致其生长速率必然呈现出时空差异性，但本文模型中由于缺少具体的实测数据，只能将这些影响综合考虑进生长速率μ0中，后续仍需进一步补充完善。

5 结 论

(1)本文通过将藻类孢子/配子的释放迁移过程引入现有附着藻类生长模型，构建了一个能正确描述附着性藻类生长迁移过程的水动力-水质模型，完善了对附着藻类全生命周期的数学描述，并根据实验数据进行了率定验证。

(2)本文采用附着藻类生长迁移模型对南水北调中线干渠工程中藻类增殖问题进行了定性研究，模拟结果显示，在本文参数条件下，干渠内水流流速小于0.4 m/s时，附着性藻类容易大量增殖，故渠道两侧边壁处的附着藻类密度普遍高于渠道中心；当增大干渠输水流量时，渠内水流平均流速会相应增大，从而抑制附着藻类的整体生长水平。

(3)通过参考相关文献和调研，确定了本文附着藻类繁殖迁移过程模型中一些关键参数的取值，但在模拟过程中仍需根据具体监测数据对现有模型参数的取值进行进一步修正；同时，模型中藻类生长受营养盐、光照等因素影响的定量描述也需进一步补充，以增强模型的适用性。