微囊藻属一日内垂向分布的数值模拟

俞 茜，刘昭伟，陈永灿，朱德军（清华大学，水沙科学与水利水电工程国家重点实验室，北京 100084）

微囊藻属一日内垂向分布的数值模拟

俞 茜，刘昭伟*，陈永灿，朱德军（清华大学，水沙科学与水利水电工程国家重点实验室，北京 100084）

将水质模型与微囊藻属浮力调节模型相结合，建立微囊藻属垂向浓度分布模型，该模型不仅保留了以往模型的优点且改进了以往模型的不足，如考虑了微囊藻群体的直径大小及微囊藻属的浮力调节能力.模型结果与荷兰的Vinkeveen湖的实际测量数据吻合良好.此外，模拟结果也表明：微囊藻属在垂向的分布受水体紊动程度影响显著.微囊藻群体的最大浓度并不在水体表面，而是在水下一定深度处.小尺寸群体微囊藻群体（20～50μm）在水深方向的分布比大尺寸群体（100～200μm）更加均匀.日落后，微囊藻属会因为密度减小上浮至水面，反之，白日里，微囊藻属密度增加向水底方向下沉，因此日出前水面的微囊藻属浓度高于日落前的水面浓度.

微囊藻属；垂向分布；浮力调节；移流传输方程；数值模拟；日变化

微囊藻属是地球上适应性最强的藻类之一［1］，在长期的进化过程中，微囊藻属演进出的伪空胞帮助其可以在垂向自如地调节浮力［1-2］，较其他藻属能适应更多极端环境［3］.此外，微囊藻属产生的微囊藻毒素会造成人体的肝脏、神经系统和皮肤疾病［4］，因此，微囊藻属是典型的有害藻属之一［5］.微囊藻属因为适应生长的温度较高［4］，所以通常在夏季爆发［6］，而随着全球气温的改变，微囊藻属水华爆发的时间可能会延长，从而带来更大的危害［4］.

除了季节性爆发这个特点以外，微囊藻属另一个较少被关注的特点为日垂向迁移能力.具有浮力调节能力的微囊藻属在一日内随着光强的改变会在水体垂向迁移，选择最合适的生长位置［1］，这种调节能力是它们与其他藻属竞争光源时的巨大优势［4］.当水下光强不能满足微囊藻属生长时，它们会通过调节浮力上浮到水体表面获取足够的光强［4］.因为一日内水面的光强改变，所以在不同时刻微囊藻属在水中所处的位置有所不同，因此在不同水深处采集的水样中所包含的微囊藻属的群体大小和数量都不尽相同.考虑到目前在较浅水体特别是河流中多为一个采样点采集一个水深处的水样［7］，因此有些样品可能不能完全反应水体中微囊藻属的真实含量，甚至造成误判.所以了解微囊藻属在一日内的垂向分布非常重要，且对于采集水样进行针对性的研究有指导意义.

综观目前的生态模型，多数成熟模型并没有考虑微囊藻属在垂向的迁移［8-9］，而专门研究微囊藻属在垂向迁移的模型也存在诸多问题［1，10-13］，如Huisman等［10-11］的模型中没有考虑微囊藻群体密度随水下光强变化而发生的改变，只给定了一个固定上浮速度，但是Ibelings等［14］在野外测量中发现湖泊中的微囊藻属在不同的水深处有不同比例的下沉.此外Huisman等［10-11］的模型也没有考虑水体表面的光强在一日内的变化，这在实际问题中是不可能存在的.而Visser等［15］和Chen等［16］的模型只研究了微囊藻群体的密度变化和轨迹，却没有研究其在垂向的浓度分布，从而限制了实际应用.野外试验发现在营养盐较为充足的水体中，水动力条件改变和悬浮颗粒物导致的水下光强衰减等是限制藻类生长的主要因素［10，17-18］，因此本研究将微囊藻属的浮力调节与水质模型相结合，建立具有浮力调节能力的微囊藻属在一日内随光强和紊动扩散系数改变的浓度垂向分布模型.该模型保留了已有模型的优点，并加入了以往微囊藻属垂向模型中没有考虑却非常重要的微囊藻属的特点，如微囊藻群体的直径大小及其垂向浮力调节能力等.本文模拟了一日内湖泊水面的微囊藻属浓度变化及其在一日内不同时刻、不同水体紊动条件下的垂向分布.

1 数学模型和数值方法

本文提出的微囊藻属垂向浓度分布模型主要由3部分内容组成：移流扩散方程［10，19］；浮力调节模型［15］；静止水体中群体上浮/下沉模型.

1.1 移流扩散方程

假设研究水体营养盐充足、温度适宜，因此模型中微囊藻属的生长速率在垂向只受光强限制，遵循Monod方程规律.此外，假设研究的深水湖泊垂向水体密度不变，即暂不考虑环境水体密度分层.微囊藻属的损失包括被细胞溶解和被动物捕捉等［11］.微囊藻属在垂向的迁移为自身上浮/下沉以及被紊动水体夹带两个作用综合的结果.式（1）～式（3）即为微囊藻属在水深方向的移流传输控制方程：

式中：C为微囊藻属的垂向浓度，cells/L；V为微囊藻群体自身在水体中的速度（V＞0，微囊藻群体下沉；V＜0，微囊藻群体上浮），m/s；w为水体垂向紊动流速，m/s；D为水体的紊动扩散系数/垂向涡旋系数，m2/s；p（I）是微囊藻属在光强为I时的生长速率，1/h；l为微囊藻属的损失速率，1/h；pmax为微囊藻属的最大生长速率，1/h；H为半饱和生长速率对应的光强，μmol photons/（m2·s）.

1.2 随光强变化浮力调节模型

Visser等［15］基于室内试验发现微囊藻属细胞密度随光强变化具有一定的规律性，并根据实验数据拟合得到细胞密度随光强变化的式（4）.白天，细胞密度改变率在光强为277.5μmol photons/（m2·s）时达到最大值，当光强未达到此数值之前，细胞密度增加率是一个增长函数，而当光强超过277.5μmol photons/（m2·s）时，细胞密度增加率为一个下降函数.在夜晚，细胞密度减少速率取决于当时的细胞密度.

式中：Im为一日内水体表面接受的最大光强，μmol photons/（m2·s）；DL为光周期，s；K为水体的消光系数，1/m；z为群体所在的水深位置，m.

1.3 静止水体中群体上浮/沉降模型

微囊藻群体在静止水体中的上浮/下沉速度采用在该领域被广泛应用的Stokes方程：

式中：g为重力加速度，一般取为9.8m/s2；ρw为水体密度，kg/m3；ρcol为微囊藻群体密度，kg/m3；d为球形微囊藻群体直径，m；v为水体的黏滞系数，kg/（m·s）；Ф为群体形状系数.

由于每个微囊藻群体都由黏液包围大量的内含有伪空胞的单细胞组成，因此微囊藻群体密度可由下式计算得到［12］：

式中：ncell是微囊藻群体中细胞所占据的比例，％；ngas为单细胞中伪空胞所占据的比例，％；ρmuc为微囊藻属体内黏液的密度，kg/m3.

1.4 数值离散

采用全隐差分格式对上述控制方程进行离散.同时，假设微囊藻属在水面和水底（z=0和z=h）处无通量交换，给定如下水面和水底的边界条件：

2 模型验证

为验证模型，将其应用在位于荷兰的Vinkeveen湖［14］，与Ibelings等［14］的实测数据进行比对.Ibelings等［14］于1989年8月23～24日期间在Vinkeveen湖（湖泊水体表面积为0.6km2，最大水深14m）监测了微囊藻群体的上浮/下沉状态.气象数据参考Medrano等［12］文中提供的有关Vinkeveen湖的数据.且根据Medrano等［12］文中的设置，假设湖泊水体密度为998kg/m3.所有的参数设置见表1.由于本文的模型每次只能设定一个群体尺寸，而Ibelings等［14］在实地测量中发现Vinkeveen湖水面的微囊藻属由多种尺寸的群体组成（图1），因此分别模拟直径为20，50，70，100，150，200μm的微囊藻群体在0m和1m的浓度，然后根据图1中的比例加权平均得到两个水深处的微囊藻属浓度［12］，再将0m和1m的浓度取平均［1］，得到图2中的模拟结果.

表1 文中所有参数取值Table 1 Parameters in the mathematical model

图1 Lake Vinkeveen湖泊表层微囊藻群体组成（图中数字表示群体直径）［14］Fig.1 Size distribution of Microcystis colonies（numbers indicate the diameters of colonies）measured in the surface layer of Lake Vinkeveen［14］

图2 Lake Vinkeveen 1989年8月23～24日水面微囊藻属生物量模拟结果与实测结果比较Fig.2 Biomass comparisons of Microcystis between the present modelling results and field data in Lake Vinkeveen on August 23-24，1989

图2表明：本文模型的模拟结果与实测数据吻合良好，且与Serizawa［1］的模拟结果相似.此外，模拟结果能够反映水体表面的微囊藻属浓度随着光强的改变进行自身垂向位置调整的变化趋势［1，4］：当水面接触的光强足够时（如8月23日的12：00～19：00和8月24日的07：00～11：00），处于光强丰富的水体表面的微囊藻属发生光合作用，体内碳水化合物含量不断积累，密度增加，微囊藻属向水底方向下沉；而当进入夜晚（如8月23日的19：00至8月24日03：00），发生呼吸作用的微囊藻群体体内的碳水化合物含量不断减少，群体密度减小，微囊藻属会上浮至水面，导致水面的微囊藻属浓度增加.由图2可见，和本文的模拟结果一致的是，Serizawa等模拟的白日水面微囊藻属生物量（8月23日的12：00～19：00和8月24日的07：00～11：00）也在不断降低，但是其减少变幅较缓（07：00～11：00），与实测值不相符，相较之下，本文的模拟结果与实测变化率更加贴合.此外，在本文的模拟结果中，微囊藻属浓度最大值出现在8月24日清晨7h左右，此时，水体表面已经有较为充足的光强，导致这个结果产生的原因可能是Wallace等［21］提出的碳水化合物的积累需要一定的反应时间.根据模拟结果可以确定：一日内水面处微囊藻属浓度峰值只有一个，出现在清晨，这个结果也与Serizawa等［1］和Medrano等［12］的模拟结果一致.

3 微囊藻属的垂向分布特征

3.1 水体紊动对于微囊藻群体垂向分布的影响

图3 直径为20μm微囊藻群体在不同紊动扩散条件的垂向分布Fig.3 Vertical concentration profiles for colonies of 20μm in diameter under different turbulent diffusion coefficients图中D的单位为cm2/s

表1中的所有参数不变，改变紊动扩散系数的数值，模拟直径为20μm的微囊藻群体在水深方向的浓度分布（图3）.由图3可见，当紊动扩散系数较低时（D=0.01cm2/s），水体紊动携带作用较弱，因此微囊藻属自身的上浮/沉降对于垂向分布起主要作用，较深水体处密度较小的微囊藻群体上浮，而上层水体中密度较大的微囊藻群体下沉，共同集中于群体密度与水体密度相同处，导致该水深处的微囊藻属浓度最大，本案例中的最大浓度约在水面以下2.4m处.而当紊动扩散系数增加（D=0.1cm2/s和D=0.5cm2/s），虽然2.4m水深处仍然是微囊藻属最大浓度所在处，但是由于水体紊动携带能力的增强，导致微囊藻属的垂向分布较为均匀.该模拟结果也解释了深水湖泊较之河流更容易发生微囊藻属水华的原因［17］.湖泊的水体紊动主要由水体表面风生成，深水湖泊在微囊藻属水华高发的夏季垂向有稳定分层，因此水体整体紊动扩散系数较小，微囊藻属易聚集在最大浓度处；相较之下，浅水河流的水体紊动主要由水底剪切力产生，加之纵向流动，导致河流整体的紊动扩散系数较高，微囊藻属不易在河流垂向某一位置集中，因此也就不易在河流中看到明显的蓝绿色水华［22］.其他尺寸的微囊藻属浓度分布与紊动扩散系数的关系与直径20μm的微囊藻群体的结果相似，本文中不再赘述.

3.2 不同尺寸微囊藻群体在不同时刻的垂向分布

表1中的所有数据不变，设置紊动扩散系数D=0.5cm2/s，并假设初始微囊藻属浓度沿水深方向均匀分布，根据第2节的分析，发现太阳升起后和太阳落山后，微囊藻属分别因为光合作用和呼吸作用导致密度发生改变，浮力会自动调节，改变其在水中的垂向位置.因此，本文模拟浓度分布稳定后一日内太阳升起时（04：00）和太阳落山时（20：00）Vinkeveen湖中直径为20，50，100，200μm的微囊藻群体在垂向的浓度分布（图4），04：00和20：00分别为当地太阳升起和落下的时刻（根据Medrano等［12］提供的当地气象数据确定）.

图4 不同直径大小的微囊藻属一日内不同时刻的垂向浓度分布Fig.4 Vertical concentration profiles for colonies of different size at different times in one day（a）4h（b）20h

图5 直径为20μm和100μm的微囊藻属在太阳升起时和太阳将落后的垂向密度和浓度分布Fig.5 Vertical density and concentration profiles for colonies of 20μm and 100μm in diameter at sunrise and at sunset

由图4可见，较小尺寸的微囊藻群体（20～50μm）在水深方向分布较为均匀，而较大尺寸的微囊藻群体（100～200μm）在水深方向分布较为集中.如直径为200μm的微囊藻属主要都集中在水下2.4m左右，水面和水深4m以下的浓度几乎为零，这与Medrano等［12］在其他案例中的模拟结果相似.因此，在较深的水深处，如6m和8m处，小尺寸微囊藻群体占据绝对优势，而随着尺寸的增加，对应尺寸群体占据的比例越小，甚至完全消失，这也与Ibelings等的观测结果一致［14］.相较于Huisman等［10-11］模型模拟得到的微囊藻群体都聚集在水体表面，本研究的结果更贴近实际.而相较于Serizawa等［1］模型中没有考虑微囊藻群体直径大小，本研究的模拟结果说明不同直径大小的微囊藻群体在垂向分布不同，证明了模型中考虑直径的重要性.

鉴于小尺寸群体和大尺寸群体的巨大差别，本文选取直径20μm和100μm分别作为代表分析微囊藻群体在太阳刚升起时和太阳落山后的密度和浓度在垂向分布情况.由图5可见，小尺寸微囊藻群体虽然沿垂向分布相对均匀，但是也和大尺寸群体一样有浓度最大处，且从图4的分布中可以大致看出浓度最大处的位置基本一致，即上文总结的群体密度与水体密度相同的位置（图5），在该水深处群体的上浮/下沉速度为零，而最大浓度处位置在一日内不同时刻也会发生变化.由图5可见，经过一夜的呼吸作用之后（20：00～04：00），微囊藻属垂向最大浓度位置向水面方向上移了一定距离，同时水面的微囊藻属浓度相应地有所增加［14］；相反，当经历了白天主要的光合作用后（04：00～20：00），微囊藻属密度明显增加（图5），群体会加快向下沉降的速度.不同尺寸的微囊藻群体在水体垂向的差异性分布再次说明群体尺寸的重要性，而这一点是Serizawa等［1］微囊藻属垂向迁移模型所没有考虑的.

4 结论

4.1 本文建立的微囊藻属浮力调节模型与移流扩散模型相结合的复合模型能够较为准确地模拟出深水湖泊水面微囊藻属昼沉夜浮的现象.夜间，水体中的微囊藻属发生呼吸作用，群体密度减少，向水面方向上浮，因此水面微囊藻属浓度增加；白日，水面的微囊藻属光合作用充分，群体密度增加，大于水体密度，下沉，因此水面微囊藻属浓度减少.此外，由于模型中考虑了Serizawa等模型中没有考虑的群体尺寸问题，使模型与实际更加贴合.

4.2 不同尺寸的微囊藻群体浓度在深水湖泊的垂向分布不尽相同，尺寸越大的微囊藻群体浓度在垂向分布越不均匀，而尺寸越小的群体浓度则容易在垂向分布相对均匀.微囊藻群体的浓度最大集中区不是在水面，而是在水下一定深度处，不同尺寸群体的最大浓度位置一致，即集中在群体密度与水体密度相同的水深处.

4.3 一日内相同尺寸的微囊藻群体在不同时刻的深水湖泊中的浓度分布也不尽相同，夜间水面微囊藻属浓度比白日浓度高，且微囊藻属最大浓度集中处会向水面方向上移一定距离.

［1］Serizawa H，Amemiya T，Rossberg A G，et al.Computer simulations of seasonal outbreak and diurnal vertical migration of cyanobacteria［J］.Limnology，2008，9（3）：185-194.

［2］Reynolds C S，Oliver R L，Walsby A E.Cyanobacterial dominance： the role of buoyancy regulation in dynamic lake environments［J］.New Zealand Journal of Marine and Freshwater Research，1987，21（3）：379-390.

［3］Ibelings B W.Changes in photosynthesis in response to combined irradiance and temperature stress in cyanobacterial surface water blooms［J］.Journal of Phycology，1996，32（4）：549-557.

［4］Paerl H W，Huisman J.Blooms like it hot［J］.Science，2008，320（5872）：57.

［5］陈永灿，俞 茜，朱德军，等.河流中浮游藻类生长的可能影响因素研究进展与展望［J］.水力发电学报，2014，33（4）：186-195.

［6］王靖国，邹 华，张 强，等.太湖微囊藻毒素的时空分布特征［J］.环境科学研究，2014，27（7）：696-703.

［7］姚绪姣，刘德富，杨正健，等.三峡水库香溪河库湾水华高发期浮游植物群落结构分布特征［J］.四川大学学报（工程科学版），2012，4（2）：211-220.

［8］Elliott J A，Irish A E，Reynolds C S.Predicting the spatial dominance of phytoplankton in a light limited and incompletely mixed eutrophic water column using the PROTECH model［J］.Freshwater Biology，2002，47（3）：433-440.

［9］Chen Q，Mynett A E.Modelling algal blooms in the Dutch coastal waters by integrated numerical and fuzzy cellular automata approaches［J］.Ecological Modelling，2006，199（1）：73-81.

［10］Huisman J，Van Oostveen P，Weissing F J.Species dynamics in phytoplankton blooms： incomplete mixing and competition for light［J］.The American Naturalist，1999，154（1）：46-68.

［11］Huisman J，Sharples J，Stroom J M，et al.Changes in turbulent mixing shift competition for light between phytoplankton species［J］.Ecology，2004，85（11）：2960-2970.

［12］Medrano E A，Uittenbogaard R E，Dionisio Pires L M，et al.Coupling hydrodynamics and buoyancy regulation in Microcystis aeruginosa for its vertical distribution in lakes［J］.Ecological Modelling，2013，248：41-56.

［13］Guven B，Howard A.Modelling the growth and movement of cyanobacteria in river systems［J］.Science of the Total Environment，2006，368（2）：898-908.

［14］Ibelings B W，Mur L R，Walsby A E.Diurnal changes in buoyancy and vertical distribution in populations of Microcystisin two shallow lakes［J］.Journal of Plankton Research，1991，13（2）：419-436.

［15］Visser P M，Passarge J，Mur L R.Modelling vertical migration of the cyanobacterium Microcystis［J］.Hydrobiologia，1997，349（1-3）：99-109.

［16］Chen Y，Qian X，Zhang Y.Modelling turbulent dispersion of buoyancy regulating cyanobacteria in wind-driven currents［C］//Bioinformatics and Biomedical Engineering，ICBBE 2009.3rd International Conference on.IEEE，2009：1-4.

［17］Yu Q，Chen Y，Liu Z，et al.The influence of a eutrophic lake to the river downstream： Spatioternporal algal compostion changes and the driving factors［J］.Water，2015，7（5）：2184-2201.

［18］朱 伟，姜谋余，赵联芳，等.悬浮泥沙对藻类生长影响的实测与分析［J］.水科学进展，2010，21（2）：241-247.

［19］Huisman J，Weissing F J.Light-limited growth and competition for light in well-mixed aquatic environments： an elementary model［J］.Ecology，1994：507-520.

［20］Reynolds C S，Jaworski G H M，Cmiech H A，et al.On the annual cycle of the blue-green alga Microcystis aeruginosa Kütz.emend.Elenkin［J］.Philosophical Transactions of the Royal Society of London.Series B，Biological Sciences，1981：419-477.

［21］Wallace B B，Bailey M，Hamilton D.Simulation of vertical position of buoyancy regulation Microcystis aeruginosa in a shallow eutrophic lake［J］.Aquatic Science，2000，62（4）：320-333.

［22］颜润润，逄 勇，赵 伟，等.环流型水域水动力对藻类生长的影响［J］.中国环境科学，2008，28（9）：813-817.

Modelling daily variation in the vertical distribution of Microcystis.

YU Qian，LIU Zhao-wei*，CHEN Yong-can，ZHU De-jun（State Key Laboratory of Hydroscience and Engineering，Tsinghua University，Beijing 100084，China）.China Environmental Science，2015，35（6）：1840～1846

By coupling advection-diffusion equation and buoyancy regulation，a coupled model for simulating diurnal changes of vertical concentration distribution of Microcystis was developed.Different from existing models，the present model incorporated the buoyancy regulation of the Microcystis colonies with different diameters.The numerical simulations showed a good agreement with the field data acquired in Lake Vinkeveen，the Netherlands.It was found that the highest concentration showed up below rather than on the water surface and turbulent diffusion had a significant impact on the vertical distribution.The colonies with smaller sizes（20～50μm in diameter）distributed more uniformly than those with larger sizes（100～200μm in diameter）.The densities of colonies decreased in the evening，and thus Microcystis float upwards.However，the densities increased and the Microcystis sank in the daytime.Hence，the surface concentration of the Microcystis reached the highest value at the sunrise and reached the lowest value at the sunset.

Microcystis；vertical distribution；buoyancy regulation；advection-diffusion equation；numerical simulation；daily variations

X703.5

A

1000-6923（2015）06-1840-07

俞 茜（1987-），女，江苏扬州人，清华大学博士研究生，主要从事水力学与水环境研究.发表论文9篇.

2014-10-27

国家自然科学基金资助项目（51039002）；教育部新世纪优秀人才支持计划（NCET-12-0309）；教育部创新团队发展计划（IRT13025）

* 责任作者，副教授，liuzhw@mail.tsinghua.edu.cn