微囊藻屬一日內垂向分布的數值模擬

俞 茜，劉昭偉，陳永燦，朱德軍（清華大學，水沙科學與水利水電工程國家重點實驗室，北京 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