输水渠系着生藻类水力生消动态建模及仿真研究

2021-01-27杨梦薇管光华刘大志

中国农村水利水电 2021年1期

杨梦薇，管光华，刘大志

(1.武汉大学水资源与水电工程科学国家重点实验室，武汉 430072；2. 湖北金浪勘察设计有限公司，武汉 430070 )

开敞明渠是我国广大灌区及调水工程的主要输水方式，经过一段时间运行的渠道系统取水口拦污栅处以及衬砌表面往往固着大量藻类生物，严重降低了渠道的输水能力以及输水效率。将这类生长在浸没于水中基质表面上的低等藻类植物称为着生藻类，主要包括蓝藻、绿藻、硅藻等多种藻类[1]。

受着生藻类生长以及混凝土表面碳化等因素影响，渠道糙率增加[2]，造成更多的沿程能量损失，从而降低渠道输水能力，且藻类的附着作用及酸性物质也会破坏渠道衬砌从而增加渗漏损失[3]。常见的除藻方法有物理、化学以及生物方法。其中物理方法为人工采用机械作用刮除水藻[4]，该方法能耗很高、效率低下，在灌区实际应用中推广价值低；化学方法则为在水源中添加化学制剂[5]例如铜制剂、O3、高锰酸钾等或者絮凝剂[6]，以抑制藻类生长，该方法需要的制剂量非常大，且对渠内水质造成二次污染，加大了水处理的成本；生物方法则通过放养喜草鱼类、微生物或高等植物，通过食物链平衡以及化感作用实现抑制藻类生长的调控措施[7, 8]，然而我国大部分灌区渠道系统为开放式，鱼类、微生物和植物的投入会给灌区管理引入新的难题，从而无法达到治理效果。

水力除藻是从水动力学角度出发，通过加大输水流量或者提高渠道平均流速以增大剪切层内水流对边界的剪应力，达到剥离水藻的目的。部分学者研究了大坝汛期洪水调度对下游藻类生物量及群落结构的影响[9]。而渠道中的水力调控方法目前在法国的普罗旺斯灌区得到应用[10]。在水力除藻的研究中，水藻并不能视为固定、静态的边界。在两次冲洗之间，藻类仍以一定的速率生长，因此需要建立藻类的生态系统模型。影响藻类初级生产力的主要因素包括：日照[11]、温度[12]、养分[13]以及水力条件[14]。贾海峰[15]等采用WASP水质模型模拟藻类生长；林春焕[16]等基于多智体建模理论，构建了藻类生长动态模拟模型；刘海东[17]等通过半隐式有限差分方法建立了藻类生长扩散的三维模型，可为Milford 湖的水质管理和藻华预警提供帮助。目前国内外的藻类生长模型主要由天然河道[18]、湖泊[17, 19]和室内水槽[16]实验率定，然而在应用于人工渠道时，由于尺度、环境条件以及藻类密度、组成的差异，实验结果适用性有待商榷，因此需要结合人工明渠的特点研究藻类生长模型并采用实测数据加以率定。

本文建立了渠道中着生藻类的生长消散模型，引入南方某灌区部分干渠作为仿真渠段用以检验模型的有效性。其中生长模型考虑了包括流速、温度、营养物和光照在内的四种影响着生藻类生长的外部因素，而消散模型则包含藻类的自然消散以及水力冲刷作用下由于剥离作用而产生的消散。论文的研究对于我国广大灌区的运行管理具有重要意义，为水力除藻的实现提供了研究基础。

1 着生藻类生消模型计算原理

1.1 渠道着生藻类生长模型

影响藻类生长的环境因素主要包括光照、温度、营养盐、流速以及水体的pH值。由于自然条件下，水中藻类可以自动调节水体pH环境以使其满足自身生长要求[20]，本文将不讨论pH值这一影响因素。论文采用单位水面面积上藻类干重度来表示着生藻类生物量的衡量指标。藻类生长用比生长速率表示，假定诸多因子对藻类比生长速率的综合影响等于各因子单个影响作用的乘积，即：

P=PmaxGuGTGNGI

(1)

式中：Pmax为最大比生长速率，其变化范围在1.00～3.00 d-1之间，模型取为1.27 d-1[21]；Gu为流速对着生藻类生长的影响因子；GT为温度对着生藻类生长的影响因子；GN为营养物质对着生藻类生长的影响因子；GI为光照对着生藻类生长的影响因子。

1.1.1 流速

一般情况下，渠道中过大的流速会冲刷渠道，使着生藻类从渠道表面剥离。然而当渠道中流速较小时，流速的适度增加对藻类生长有促进作用，其作用机理为：渠道水流流速较小时，渠道水流处于层流状态，该水流状态中营养物质的稀释扩散速度相对于紊流更慢，使得藻类细胞可以吸收利用的营养物质在层流中比在紊流中多，因而可以在一定程度上促进藻类生长[22]。对于常年水动力条件基本不变且流速很小的水域，如湖泊、水库等，流速对生长于其中藻类的影响可以忽略不计，则Gu=1[21]。而由于本研究在灌溉渠道中进行，灌溉期与非灌溉期的渠道水动力条件差异较大，则无法忽略流速的影响。

由上述分析可知，在一定流速区间内增大流速会促进藻类生长，即藻类生长存在一个临界流速，该流速较小但不等于零。论文假定一个临界流速a，流速影响因子计算如下[21]：

(2)

式中：u为渠道内流速，m/s;a为临界流速，其变化范围为0.01～0.03 m/s，论文取a=0.018 m/s；b为流速影响参数，变化范围为0.1～0.3 m2/s2，论文取b=0.238 m2/s2。

1.1.2 温度

温度通过控制藻类光合作用的酶促反应、营养物质的吸收和呼吸作用之间的平衡关系，来直接影响藻类的生理活动以及生长繁殖，并通过控制水体中各类营养物质的溶解度、离解度或分解率等物理化学过程来间接地影响藻类的生长。

论文由于缺少仿真渠道实测水温资料，且着生藻类多半生长于靠近水面处，可假定生长藻类的水域水温和当地气温一致，即利用气温资料来代替水温资料。论文所选取的仿真渠道处于亚热带季风气候区，多年平均气温为16.6 ℃，取渠道中藻类生长的最适温度为20 ℃。温度对着生藻类生长的影响[23]如下：

GT=θT-Topt

(3)

式中：T为渠道水温，℃；θ为温度修正系数，取值在1.0～1.5之间，论文取θ=1.12[24]；Topt为藻类生长的最适温度，取Topt=20 ℃。

1.1.3 营养物

水体的营养状况是调节藻类生长、繁殖和内部生物化学反应的基本因素之一。藻类为了满足正常的生理活动所需要的元素有14～21种，而在这些元素之中，对藻类生长限制最为显著的是氮和磷[25]。渠道水流可分为两层，上层流速较快即为流动水域，下层流速较慢即为静水区域。藻类均生长于静水区域中，因而静水区域中的营养物质浓度直接影响和限制着生藻类的生长。静水区域中营养物质的浓度是小于上层流动水域的，因此营养物质就会在这两层水域之中进行交换，即渠道中藻类营养物质的外部来源[26]。着生藻群内部营养物质的来源涉及许多复杂的生物变化，将该生物变化简化如下：首先假定着生藻类内部营养物质的变化主要取决于藻类的繁殖；然后假定，在藻类细胞的繁殖过程中，细胞中的营养物质全部由母细胞进入到子细胞中。基于上述假定，着生藻类细胞对营养物质的消耗则主要来源于呼吸作用[27]。

论文依照水中营养物质过少会影响藻类生长的限制原则，运用Mickaelis Menten模型[28]计算营养物质对着生藻类生长的影响因子如下：

(4)

式中：KP为磷的半饱和常数，取KP=0.33 mg/L；KN为氮的半饱和常数，取KN=3.8 mg/L；TP为渠道中磷元素含量, mg/L；TN为渠道中氮元素含量, mg/L。

1.1.4 光照

藻类利用太阳能将自然界中丰富的二氧化碳和水合成有机物，同时提供动物、植物以及多种微生物所必需的氧气，该过程称为光合作用。光照通过影响藻类的光合作用强弱进而影响藻类生长。计算光照对着生藻类的生长影响系数主要包括以下几个方面：

首先计算太阳总辐射。太阳辐射总量是指某一接受面在一日、一季度、一年或一时段内所接受的太阳直接辐射和散射辐射之和。本论文为简化计算，假定大气对太阳辐射不发生削弱作用，则一天之中太阳直接辐射总量按照下式计算[29]：

(5)

式中：So为太阳常数；T为单位时间，1 d；wo为太阳时角，可由cosw0=-tgφtgδ求出后[29]，化为弧度；δ为赤纬，为太阳的平行光线垂直照射地球表面所在的纬度值，数值上等于阳光直射点所在地的地理纬度，北半球取正值，南半球取负值。赤纬随季节变化，其值由近似公式δ=23.5 sinN求得[30]，其中N以度为单位，其值等于距春分日或秋分日最近的总天数，春分日至秋分日取正值，反之取负值；φ为当地的地理纬度[31]。

由式(5)可以得到理论太阳辐射的总量。然而实际的辐射总量除了受到太阳高度角、大气的吸收和散射的影响，还与云量和海拔高度有关，其计算采用经验公式如下[32]：

(6)

式中：QO为日总辐射量；n为每日的日照时数；No为日可照时数；a和b均为经验系数，a=0.248，b=0.752。

对于散射辐射，由于受到太阳高度角、大气厚度、纬度、海拔、云量等众多因素的影响，很难用一个公式去计算。在粗略计算中可以将散射辐射日总量的计算看成是直接辐射日总量的1/2[29]:

D=0.5Q

(7)

植物的光合作用光谱比短波辐射的整个波长区(300～3 000 nm)狭窄得多，光合有效辐射的波段大体在可见光(400～760 nm)范围以内，即400～700 nm或380～710 nm的范围。论文通过气象学计算的方法求解光合有效辐射，计算光合有效辐射的简便公式为[29]：

Qφ=0.43Q+0.57D

(8)

式中：Qφ为光合有效辐射；Q为水平面上的直接辐射；D为散射辐射。

将辐照度与光照度进行转化[33]：

ΦV(λ)=KmV(λ)Qφ

(9)

式中：ΦV(λ)为波长为λ的光的辐射通量；Km为最大光视效能；V(λ)为光视效率。

然后考虑光在水中的衰减。由于水体中含有悬浮的有机物质和无机物质的阻挡、水体的颜色、被藻类吸收以及水对光的散射作用等情况的出现，光在穿透水体时会呈指数形式的衰减。论文将水体从顶部到底部共分为M层，第z层可用于光合作用的光照幅度Iz就可通过Beer-Lambert方程求解[34]：

(10)

式中：kbg为在没有水生生物时的浊度(即水的遮光度)，论文研究的是渠道中的水流，主要用于城市用水以及灌溉用水，则取kbg=1.2 m-1；k为上层浮游植物吸收光照的系数；Io为入射到水面的光照强度，Io=ΦV(λ)；P为在s层，时间t时的水中生物总量，则从0到z的积分其意义就是在t时刻，0到z层水中生物的总量，记为Qz。

在本论文中，对Qz的计算进行简化。假设水体中的生物量是均匀分布，则0到z层水中生物的总量可以表示为：

(11)

式中：QM为0到M层水中生物的总量，即水体由顶部到底部的生物总量。

不考虑藻类对光的吸收，只考虑可见光在水中的衰减，则式(10)简化为如下形式：

(12)

求得每一层的光照强度后，采用下式求出平均光照强度[34]：

(13)

式中：h为渠道中的水深。

最终得到光照对着生藻类生长的影响计算如下[35]：

(14)

式中：Iopt为最适光照强度。

1.2 渠道着生藻类消散模型

渠道中藻类生长的同时也伴随着藻类消散，藻类消散主要分为两个部分：其一为藻类的自然消亡，包括藻类的自然死亡以及生物的捕食，本论文仿真环境为灌溉渠道，忽略生物捕食作用，则该部分以自然死亡为主；其二为渠道中的水动力条件改变导致的着生藻类剥离。

论文藻类的消亡模型部分计算死亡速率公式如下[36]：

(15)

式中：K1R是内源呼吸速率常数，取K1R=0.125 d-1，即30 ℃时浮游植物呼吸速率常数；θ1R为温度修正系数，取θ1R=1.045；K1D为死亡速率常数，主要是指由于细菌感染、物理性破坏和藻类生长的自然死亡等过程引起的死亡，取K1D=0.02 d-1[35]。

藻类分离率模型如下[37, 38]：

(16)

式中：D为藻类分离生物量；Cd为藻类分离系数，取Cd=0.035；τo为渠道剪切力， N；τ0,cr为使藻类分离的临界剪应力，取τ0,cr=0.087 N；β为指数系数，取β=1.5；B为渠道生物量；Br为粗糙的渠道表面隐藏的无法冲刷的藻类，即残留生物量。

渠道剪切力计算公式如下[39]：

(17)

式中：γ为水容重， N/m3；n为渠道糙率；h为渠道水深，m；v为渠道流速， m/s。

2 目标渠段建模及仿真分析

2.1 环境参数

论文选用南方某灌区部分干渠为仿真渠道，其环境参数选取及仿真结果如下。

2.1.1 流速

为便于仿真计算，论文在进行藻类生长计算时，假定渠道各处流速均相同，且根据仿真渠道流速数据取v= 0.15 m/s。将其代入式(2)中，计算得到流速对着生藻类生长的影响因子Gu=0.929 4。

2.1.2 温度

为了更好地反映出藻类的生长期、平衡期以及下降期的特点，论文选用仿真渠段2014年3月至11月的观测温度资料代入式(3)中计算。绘制温度T和影响因子GT随时间变化过程图如图1所示。

由图 1可知，影响因子GT的取值随着温度T的变化而改变。在温度较低的初春和深秋，GT的数值较低，即低温对着生藻类生长的抑制作用较为明显；而在温度较高且气温比较接近藻类最佳生长温度的夏季，GT的数值则相对较高，即该时段的温度有助于藻类生长。

2.1.3 营养物

水流中促进藻类生长的营养物质主要是含有氮、磷两种元素的营养盐类。考虑到仿真渠段没有实时监测的水质资料，论文选用其取水水库坝首的监测数据作为水质资料进行计算。将数据代入式(4)中，得到总氮TN、总磷TP以及影响因子GN随时间变化过程如图 2所示。

由图中3条曲线的变化规律可知，GN的波动趋势与TP大致相同，则渠道中着生藻类的生长主要受水流中磷元素含量的限制，因为水流中的磷元素相较氮元素其含量是不足的。从整体上看，GN的数值较小，则仿真渠段中营养物质的含量较低，其对着生藻类生长的促进作用较小，水库的水质良好。

2.1.4 光照

由于无法实时监测渠道水面的光照强度，论文通过数值计算的方法将日照时数转化为光照强度。其中日照时数资料来自2014年《地面气象记录月报表》。仿真从2014年3月1日开始，至2014年11月30日结束。2014年的春分日是3月21日，秋分日是9月23日，代入赤尾δ的近似计算公式中，进而通过式(5)计算得到太阳直接辐射总量Q0。通过式(6)至式(14)的一系列计算，最终得到日照时数与影响因子GI随时间的变化过程图如图3所示。其中式(14)中的最适光照强度理论上应通过室内的着生藻类培养实验进行确定，而论文由于时间和试验条件的限制无法进行实验，则选用香溪河中的藻类生长最适光照强度4 700 lx[35]代入计算。

由图3可知，GI随着日照时数的变化而变化。日照时数较小时，GI的数值也相应较小；日照时数较大时，GI也相应较大。这是由于阴雨天的太阳辐射较小，不利于藻类的光合作用，进而抑制了藻类的生长。此外，在仿真中期156～165 d时段，日照时数都较大，但GI的数值却不升反降，这是因为夏天光照强度太大，以至于光照强度远大于最适光照强度，抑制了藻类生长。

2.2 仿真结果

2.2.1 死亡率

仿真得到的死亡率由自然死亡率和水力剥离两部分组成，可由式(15)计算得到死亡率与温度随时间变化过程如图4所示。

由图4可得，死亡率与温度的变化趋势完全一致。温度较高时，死亡率也相应较大；温度较低时，死亡率也相应较小。这是因为温度的高低影响着藻类的呼吸作用，且呼吸作用对高温不敏感，只有当温度超过40 ℃时才会抑制藻类的呼吸作用，低于40 ℃时呼吸强度与温度的变化规律一致，因此图中呈现出死亡率与温度变化一致的规律。

2.2.2 生物量

在确定生长模型和自然消亡模型后就可以仿真在自然环境下渠道中着生藻类生物量的变化。由式(18)可计算得到仿真渠段内的生物量：

B(t+1)=(1.27×Gu×GT×GN×GI-

0.125×1.045(T-30)-0.02)B(t)

(18)

式中：B(t)和B(t+1)分别为t时刻和t+1时刻对应的仿真渠段中的生物量。

将计算得到的生物量随时间的变化过程作图如图5所示。图5定量展现了着生藻类的生长和消散过程。其中藻类的大量繁殖发生于仿真开始后149 d左右，此时正值7月末，为夏季，该现象与所预测的夏天藻类会大量繁殖的结论一致。而在仿真开始后244 d左右，生物量恢复至较低水平，此时正值10月末，进入秋季，温度降低，光照强度和日照时数均下降，导致藻类的光合作用强度下降，藻类生物量减少。分析生物量的数值变化，观察到仿真初期生物量是有所下降的，这是因为仿真初期温度较低，光照强度和日照时数较短，导致藻类光合作用小于其呼吸作用，故而生物量下降。而随着仿真时间的增加，温度逐渐上升，光照因素也逐渐有利于光合作用，藻类的数量才重新开始增加，其最大值为4.629 8 kg/m2。该数值为初始值的3 307 倍，则在没有水流冲刷时，渠道中的藻类会疯狂生长。此外，在严寒的冬季渠道中的藻类数量并不完全为零，有极少数的藻类会保存下来，使得在仿真初期渠道中藻类数量并不为零而是一个较小的数值。

2.3 参数敏感性分析

渠道藻类生长与环境因子密切相关。定义O为模型输出结果，在论文中输出结果为仿真时段内藻类生物量最大值。选取用于计算环境影响系数的物理量作为O的影响因子，包括流速u、温度T、总氮TN、总磷TP以及日照时数与可照时数的比值RN。

设数学模型基本表达式为：

O=f(F1，F2，…，Fn)

(19)

式中：O为模型模拟输出结果；Fi为影响O的因子。

本文采用参数敏感性局部分析法[40]，参数敏感性指标计算使用下式：

(20)

式中：等式右边表示模型输出结果O对Fi因子变化ΔFi的响应S。

为了对参数之间的敏感性进行对比，常采用的方法是使用相对敏感性值，进行参数敏感性的归一化，计算敏感性指标I。则式(20)变换为：

(21)

敏感性指标I不受O和Fi的单位尺度的影响，使模型参数之间具有可比性。

对原始数据的流速u、温度T、总氮TN、总磷TP以及日照时数与可照时数的比值RN的逐日数据的平均值按照+10%、+5%、+2.5%、+1%、-1%、-2.5%的比例求得变化值ΔFi，并叠加至逐日数据中，通过仿真求得生物量的变化过程，并得到生物量B在变化过程中的最大值作为输出量。由式(20)计算得到各影响因子对应的敏感性指标|I|见表1，敏感性分类见表2。