基于水流擾動參數的復雜河湖系統藻類數值模擬

閔志華+王曉琴+李傲

摘要： 為研究水流擾動對藻類生長的影響，基于室內藻類生長動力學試驗，提出藻類生長過程中的水流擾動參數α，并對藻類生長公式進行修正；運用一維河網與二維湖泊耦合的藻類生長動力學模型對太湖西北部重污染區藻類生長的時空分布特征進行動態模擬。結果表明：（1） 水動力條件對藻類生長影響明顯，研究區內東西向河流流速較大，水流擾動強，藻類濃度低于南北向河流；（2）各研究斷面藻類濃度均呈現時間分布不均性，枯季（11月至次年4月）水流擾動較弱，為藻類生長創造較好的環境，藻類平均濃度較洪季（5月至10月）增加了9.22%；（3）河湖系統內湖口區水流條件更適宜藻類生長，位于湖口區的S2與S4點位濃度平均較河網區增加了18.5%。

關鍵詞： 太湖；河湖系統；富營養化；藻類；擾動；模型

中圖分類號： X171 文獻標志碼： A

文章編號：1002-1302（2017）22-0264-05

水體富營養化導致的水華暴發是全球范圍內面臨的一項重要環境問題[1]，藻類等浮游植物的大量生長是其發生的重要原因[2]。藻類生長數值模擬是藍藻水華預警的重要基礎性工作，對于維持健康穩定的水生態系統具有重要意義。諸多學者開展了藻類數值模擬研究。逄勇等進行了太湖藻類的數值模擬，探討了太湖藻類的動態變化機制[3-4]；Cole等構建了包括溫度、營養鹽及溶解氧等狀態變量在內的富營養化模型，開展了湖灣及淺水湖泊的藻類生長數值模擬[5-7]。這些成果提高了藻類數值模擬的水平，但基本沒有考慮水流擾動對藻類生長數值模擬的影響。同時，關于藻類數值模擬的研究大多圍繞單一湖泊或河道開展，對于河網縱橫交錯、湖泊較多的復雜河湖系統所進行的藻類數值模擬研究未見報道，這類地區湖泊與河網之間存在著頻繁的物質與能量交換，藻類數值模擬更為復雜[8-9]。本研究提出了影響藻類生長的水流擾動參數α，運用一維河網與二維湖泊耦合的藻類生長動力學模型對太湖西北部復雜河湖系統藻類生長的時空分布特征進行了動態模擬，旨在為藻類生長數值模擬和水華預警提供依據。

1 研究區域

研究區為太湖西北部重污染區，該區域河流眾多并有多個湖泊，屬于典型的河湖系統。研究范圍北起京杭大運河，南至江蘇省宜興市東、西氿，西至滆湖，東至太湖梅梁灣和竺山灣，涉及江蘇省無錫市惠山、北塘、南長、新區、崇安、濱湖等6個區和宜興市，以及江蘇省常州市鐘樓、天寧、戚墅堰、新北、武進等5個區，共41個鎮，56個街道。研究區域總面積為 5 272 km2，占太湖流域的5.4%。該河網區水質達標率為17%，占太湖入湖污染物通量的 80% 左右[10]，是影響太湖湖體，特別是梅梁灣、竺山灣水質的主要區域（圖1）。該河湖系統內主要有12條河流，其中入太湖河流7條，入太湖竺山灣的河流分別為太滆運河、太滆南運河、漕橋河、官瀆港、陳東港，入太湖梅梁灣的河流分別為直湖港、武進港；出滆湖河流4條，分別為太滆運河、漕橋河、太滆南運河（殷村港）、燒香港。研究區內河道見表1。河湖系統內主要湖泊有滆湖、東氿、西氿、太湖竺山灣、梅梁灣。滆湖是該區域重要的行蓄洪湖泊，位于太湖西北部，進出水量基本平衡，換水周期較短，為每年7.03次。東氿、西氿位于太湖以西，滆湖以南，是流域眾多河流入湖前所必經的重要湖泊。竺山灣和梅梁灣為太湖北部的2個湖灣，竺山灣為太湖西北部的半封閉性富營養湖灣，北起百瀆口，南至馬山咀一線；梅梁灣是太湖北部的湖灣，占太湖總面積的5.2%。2個湖灣是承接經入湖河流攜帶的大量陸域污染物的主要水域。研究區各湖泊基本信息見表2。

2 藻類生長動力學試驗

2.1 試驗條件與方法

為了研究水流擾動對藻類生長的影響，試驗過程中維持溫度、光照、水質等藻類生長影響因子的一致性，通過改變恒溫振蕩箱轉速來代替水流擾動。試驗過程中選擇引起藍藻暴發的優勢藻種銅綠微囊藻作為試驗藻種[11-12]。藻種來自中國科學院水生生物研究所藻種庫，保存在20 ℃的光照培養箱中。運用pH值為8.6的MA培養基對銅綠微囊藻進行單獨培養[13]。在具硅膠塞的250 mL錐形瓶中加入100 mL培養基，高壓滅菌，將預先培養至對數生長期的藻種離心、洗滌后接種，使得試驗培養液中藻類的初始密度為 5.0×104 cell/mL，接種后將錐形瓶放入恒溫振蕩箱培養；布置6個培養樣品，每個樣品分別設置轉速0、60、90、120、150、200 r/min，各轉速重復試驗2次；試驗過程所有操作均在無菌條件下進行，不同擾動強度下，均維持溫度為25 ℃，光照度為2 000 lx，光—暗周期為12 h—12 h。自接種當日起，每天在同一時間取樣， 每種樣品取0.1 mL藻液置于0.1 mL微囊計數

板中用雙目顯微鏡進行多次計數，取平均值，直到藻衰敗為止。以藻類生長期內平均比增長率作為藻類在某條件下增長率的估計值[14]，公式為

式中：μ為藻細胞比增長率；n為藻類生長期時間；xi為當日藻細胞數；xi-1為前1 d的藻細胞數。

2.2 試驗結果分析

根據試驗結果，繪制不同擾動強度下藻類生長過程線及平均增長率（圖2）。結果表明：（1）在不同轉速下，銅綠微囊藻生長情況隨時間分布規律基本一致，可分為生長期（第1～5天）、穩定期（第5～8天）、衰亡期（第8天后）。（2）在靜止條件下，藻類生長過程中沒有明顯的穩定期，僅有生長期與衰亡期，且生長期較擾動條件下延長了1倍。（3）在靜止條件下，藻細胞數為0.7×106個/mL；當轉速為60、90 r/min時，藻細胞數明顯增加，分別達1.50×106、0.84×106個/mL；當轉速增加至120 r/min以上時，藻細胞數明顯減少。（4）在靜止條件下，銅綠微囊藻增長率較低，為0.29；當轉速為60、90 r/min 時，增長率分別達0.57、0.71；當轉速增加至 120 r/min 以上時，藻類增長率呈現降低趨勢。綜上所述，擾動轉速對銅綠微囊藻的生長有影響，適當的擾動轉速可以促進銅綠微囊藻的生長，而靜止狀態和過高的擾動轉速對藻類生長有抑制作用。endprint

2.3 水流擾動參數確定

由于室內試驗條件與野外環境存在差異，如何將室內擾動轉速轉換成野外環境中的水體流速尤為關鍵。本研究選擇起懸泥沙濃度作為中間因子，通過野外同步水流泥沙監測及同步泥沙起懸試驗，建立了擾動轉速與水體流速的等量轉化關系。以上轉化建立在“引起相同起懸泥沙濃度的室內擾動轉速與水體流速其動力效應是相當的”。野外同步監測分別在2012年1月（枯季）與7月（洪季）2個典型季節開展，共設置了24個同步監測點位，其中湖體6個，河網區18個（圖3）。基于各監測垂線0.8H處的水流泥沙監測結果，建立水體流速與起懸泥沙濃度關系。室內泥沙起懸試驗同樣設置0、60、90、120、150、200 r/min共6個轉速處理，試驗用水和泥樣完全取自野外監測點位；每組擾動設置2個重復開展起懸試驗，根據試驗結果，建立室內擾動強度與起懸泥沙濃度關系曲線。基于野外實際流速、室內擾動強度與起懸泥沙濃度的相關性建立起兩者間的轉換方程；結合藻類生長動力學試驗結果，繪制藻類生長水流擾動參數曲線（圖4），具體數學方程如下：

式中：α為水流擾動參數；f（u）與g（u）分別為不同流速區間的相關函數；u為水體流速，m/s；uopt為最適宜水體流速，根據試驗結果，取0.15 m/s。

3 河湖藻類數學模型

3.1 基本方程

太湖西北部重污染區河湖系統水流復雜，藻類數值模擬采用一維河網與二維湖泊耦合模型。一維模型實現河網區水流過程與藻類生長過程的模擬計算，同時為二維湖泊模型提供相關的計算邊界；二維模型計算范圍主要為系統內河網受納湖泊，包括東氿、西氿、太湖竺山灣和梅梁灣，二維模型用于河湖交換條件下湖泊藻類生長的數值模擬。

3.1.1 一維河網藻類生長模型

一維模型選用TN、TP、Chl-a作為主要計算指標，基于生態動力學方程建立綜合考慮光照、溫度、營養鹽及水動力條件等多因子耦合的一維河網區藻類生長模型，其中光照、溫度等參數在藻類生長參數模塊中統一考慮，模型基本控制方程如下：

式中：t為時間坐標，s；x為順河向長度坐標，m；Q為流量，m3/s；Z為水位，m；U為斷面平均流速，m/s；n為糙率；A為過流斷面面積，m2；B為主流斷面寬度，m；BW為水面寬度，m；R為水力半徑，m；g為重力加速度，m2/s；q為旁側入流流量，m3/s。N、P、Chl-a為斷面總氮、總磷、葉綠素a的垂線平均濃度，mg/L；Ex為縱向分散系數；Kn、KP為N、P污染物降解系數；SChl-a、SN、SP為葉綠素a和氮、磷的源匯項。

3.1.2 二維湖泊藻類生長模型

二維湖泊藻類生長模型在有限體積法框架下構建，模型同樣耦合了光照、溫度、營養鹽及水動力條件等多項藻類生長影響因子，具體控制方程如下：

式中：h為水深；u、v為x、y方向垂線平均水平流速分量；g為重力加速度；s0x為x向的河底底坡；s0y為y向的河底底坡；sfx為x向摩阻底坡；sfy為y向摩阻底坡；ΔT為溫差；tx、ty分別為水溫在x、y向擴散系數；ρ和Cp分別為水的密度和比熱；Ks為水面綜合散熱系數，模型采用Gunnerberg經驗公式計算水面綜合散熱系數。N、P、Chl-a分別為總氮、總磷、葉綠素a的垂線平均濃度；SN、SP、SChl-a分別為氮、磷和葉綠素a的源匯項。

方程（3）、（4）中，SN、SP、SChl-a可具體表達為以下形式：

式中：RN、RP分別為底泥中氮、磷的釋放項；SEDN、SEDP分別為氮、磷的沉降項；CN、CP分別為藻類死亡分解對氮、磷的貢獻項；UN、UP分別為藻類對氮、磷的吸收項；μ為藻類的生長率；sed為藻類的沉降率；d為藻類的損失率。水流擾動參數α對藻類生長率μ進行修正如下：

式中：μmax為藻類最大生長率；f（T）為溫度函數；f（I）為光照函數；f（N）與f（P）為氮、磷營養鹽函數；α為水流擾動參數，按照公式（2）計算取值。

3.2 數值解法

河湖系統耦合模型采用分段求解。一維河網水動力微分方程組采用四點隱式差分格式數值求解，對水動力方程組以Preissmann四點線性隱式差分格式將其離散，輔以連接條件，形成河道方程，以微段、河段、汊點三級聯解的方法[15]求解。TN、TP、Chl-a求解采用隱式迎風格式，時間項采用前差分對每個單一河道進行離散[16]，將整個河網TN、TP、Chl-a濃度離散成節點和斷面水質濃度進行數值求解。二維湖泊模型將水流水質方程與藻類生長方程聯合求解，公式（4）可統一寫為以下形式[17-19]：

式中：q為守恒物理量；f（q）、g（q）分別為x、y方向通量；b（q）為源匯項。

在任意形狀的單元Ω上采用有限體積法對式（4）進行積分離散，運用FVS格式求解法向數值通量。具體解算過程及模型率定驗證過程詳見參考文獻[20-21]。

3.3 率定驗證

基于2012年1月（枯季）、7月（洪季）的水文水質同步監測結果，對所建耦合數學模型進行率定驗證。計算區域為研究區域內的河湖系統，運用Gambit軟件對湖區進行無結構網格剖分。以西部滆湖、北部京杭運河實測水量水質作為計算邊界。考慮計算穩定性及精度，取時間步長為1 s。結果表明，計算結果與實測值擬合效果較好，COD、氨氮的平均相對誤差分別為22.0%、17.9%。所建模型能夠較準確地反映河湖系統水量水質動態變化特征。

4 河湖系統藻類時空分布規律

基于所建的考慮水流擾動參數的一維、二維耦合藻類模型，對太湖西北部重污染區河湖系統藻類時空分布規律進行數值模擬。根據太湖流域茅東、丹陽、溧陽、趙村、宜興、常州、陳墅等40多個雨量站近10年時間序列內降水量資料頻率分析結果，確定2011平水年為計算典型年，以西部滆湖、北部京杭運河實測TN、TP、Chl-a濃度作為計算初始邊界，對系統內藻類生長及遷移過程進行數值模擬。選擇6個典型斷面對藻類時空分布特征進行分析，分別為太滆南運河上游S1 （31°30′37.18″N，119°48′45.92″E）、竺山灣湖口S2 （31°27′12.31″N，120°0′25.39″E）、直湖港河流上游S3 （31°40′20.14″N，120°8′12.13″E）、梅粱灣湖口S4 （31°30′31.44″N，120°7′2167″E）、永安河下游S5 （31°42′8.90″N，119°59′14.91″E）、太滆運河S6 （31°34′20.31″N，119°57′14.94″E），其中S2與S4為湖泊點位，其余為河網點位（圖5）。endprint

計算結果表明：（1） 平水年水文條件下，S1～S6斷面全年平均Chl-a濃度分別為8.50、10.32、10.02、10.92、11.49、9.15 mg/m3；S1與S6斷面藻類水平相對較低，而S4與S5斷面Chl-a濃度明顯增加，S2與S3斷面濃度處于中間水平（圖6）。根據多年實測水文數據，研究區域內東西向河流流速大于南北向河流，太滆運河、太滆南運河等東西流向河流全年平均流速0.23 m/s，南北流向的直湖港、永安河等河流全年平均流速為0.16 m/s；S1與S6均位于東西向河流上游，水流擾動較強，故藻類濃度最低；而S5斷面位于南北向永安河流下游，水動力條件較適合藻類生長，故Chl-a濃度最高。（2） 6個斷面Chl-a濃度過程都呈現時間分布不均性，洪季5—10月水動力條件較好，Chl-a濃度相對較低，S1～S6斷面洪季平均Chl-a濃度分別為8.15、10.07、9.68、10.57、1103、8.80 mg/m3；而枯季11月至次年4月，水流擾動減弱為藻類生長創造較好環境，Chl-a濃度有所增加，6個斷面枯季平均水平分別較洪季增加了9.26%、5.58%、7.42%、6.58%、9.18%、8.09%。S1與S5斷面洪枯2季水流差異較大，Chl-a濃度變化最為明顯，而S2與S4位于湖口區，Chl-a濃度隨季節變化幅度相對較低。（3）6個斷面中，S2與S4斷面分別位于太滆南運河、直湖港下游入湖口，其Chl-a濃度較其上游S1與S3斷面分別增加了28.0%、8.9%；湖泊區水流條件更適宜藻類生長，Chl-a水平整體高于河網區。根據孔繁翔等對大型淺水湖泊中藍藻水華形成機理的研究結果[21]，以Chl-a濃度10 mg/m3作為水華暴發臨界值。S1與S6斷面Chl-a濃度，除枯季局部時段，基本未達水華暴發水平；S4與S5斷面Chl-a濃度則超過了臨界值；S2與S3斷面Chl-a濃度在臨界值附近波動，枯季Chl-a濃度較高時會超過臨界水平，而洪季又有所降低。

5 結論

以太湖西北部重污染區的復雜河湖系統為例，針對該區域水系連通發達，水動力條件復雜多變的特征，對水動力影響下的藻類生長過程進行了數值模擬，得到以下結論：（1） 水動力條件的變化會影響藻類生長，適當的動力擾動能促進藻類生長。（2） 提出了水動力影響參數α，對藻類生長公式進行了修正，完善了藻類生長數學模型，從而更合理地對水體富營養化進行模擬和預測。（3）河湖系統內河網區藻類生長與水動力條件關系密切，區內Chl-a濃度呈現顯著時空分布不均性。洪季水流擾動較強，Chl-a濃度低于枯季。區內南北向河流的水體流速低于東西向河流，其水動力條件更適宜藻類生長，Chl-a濃度相對較高。研究成果對于區域水華預測預警具有重要指導意義，也為太湖流域解決水體富營養化問題提供了參考。

參考文獻：

[1] Lürling M，van Oosterhout F. Controlling eutrophication by combined bloom precipitation and sediment phosphorus inactivation[J]. Water Research，2013，47（17）：6527-6537.

[2]Guo L. Doing battle with the green monster of Lake Taihu[J]. Science，2007，317（5842）：1166-1166.

[3]逄 勇，丁 玲，高 光. 基于生態槽實驗的藻類生長參數確定[J]. 環境科學，2005，26（3）：78-82.

[4]張巖松. 北京水系多藻類生態動力學模型[J]. 清華大學學報（自然科學版），2009（12）：1992-1996.

[5]Cole F C. Three-dimensional eutrophication model of Chesapeak Bay[J]. Journal of Environmental Engineering，1993，119：1006-1025.

[6]Cerco C F，Noel M R. Twenty-one-year simulation of Chesapeake Bay water quality using the CE-QUAL-ICM eutrophication model[J]. Journal of the American Water Resources Association，2013，49（5）：1119-1133.

[7]Nyholm N. A simulation model for phytoplankton growth and nutrient cycling in eutrophic，shallow lakes[J]. Ecological Modelling，1978，4（3）：279-310.

[8]郝秀平，王中根. 河湖水系連通分類體系探討[J]. 自然資源學報，2011（11）：1975-1982.

[9]Zhang Y，Zhang Y，Gao Y，et al. Water pollution control technology and strategy for river-lake systems：a case study in Gehu Lake and Taige Canal[J]. Ecotoxicology，2011，20（5）：1154-1159.

[10] 馬 倩，劉俊杰，高明遠. 江蘇省入太湖污染量分析（1998—2007年）[J]. 湖泊科學，2010，22（1）：29-34.

[11]史小麗，王鳳平，蔣麗娟，等. 擾動對外源性磷在模擬水生態系統中遷移的影響[J]. 中國環境科學，2002，22（6）：537-541.

[12]章宗涉，劉永定. 微囊藻柵藻資源競爭的動力學過程 Ⅰ. 光能和磷營養的半飽和參數及其生長率動態[J]. 環境科學學報，2000，20（3）：349-354.

[13]朱廣偉，秦伯強，高 光. 強弱風浪擾動下太湖的營養鹽垂向分布特征[J]. 水科學進展，2004，15（6）：775-780.

[14]張東生. 河網非恒定流的三級聯合解法[J]. 華東水利學院學報，1982，10（1）：1-13.

[15]Dewey R J. Application of stochastic dissolved oxygen model[J]. Journal of Environmental Engineering，1984，110（2）：412-429.

[16]Zhao D H. Finite-volume two-dimensional unsteady-flow model for river basins[J]. Journal of Hydraulic Engineering，1994，120（7）：863-883.

[17]趙棣華，戚 晨，庾維德，等. 平面二維水流-水質有限體積法及黎曼近似解模型[J]. 水科學進展，2000，11（4）：368-374.

[18]Zhao D H. Approximate riemann solvers in FVM for 2D hydraulic shock wave modeling[J]. Journal of Hydraulic Engineering，1996，122（12）：692-702.

[19]譚維炎. 無結構網格上二維淺水流動的數值模擬[J]. 水科學進展，1995，6（1）：1-9.

[20]丁 玲. 水體透明度模型及其在沉水植物恢復中的應用研究[D]. 南京：河海大學，2006.

[21]孔繁翔，高 光. 大型淺水富營養化湖泊中藍藻水華形成機理的思考[J]. 生態學報，2005，25（3）：589-595.endprint