王海燕， 陳凱麒， 祁昌軍，張貝貝， 黃 茹， 溫靜雅

(1:環(huán)境保護(hù)部環(huán)境工程評(píng)估中心，北京 100012;2．國(guó)家環(huán)境保護(hù)環(huán)境影響評(píng)價(jià)數(shù)值模擬重點(diǎn)實(shí)驗(yàn)室，北京 100012;3.水電環(huán)境研究院，北京 100012; 4.河海大學(xué)水文水資源學(xué)院,南京 210098)

1 引 言

當(dāng)前，數(shù)值模擬已經(jīng)成為研究水環(huán)境問(wèn)題的一項(xiàng)重要手段，以往針對(duì)淺水湖泊的研究，通常采用二維數(shù)值模擬以減少計(jì)算量，同時(shí)因?qū)崪y(cè)數(shù)據(jù)難以獲取而無(wú)法進(jìn)行更深層次的機(jī)理性研究。隨著計(jì)算機(jī)水平發(fā)展和數(shù)據(jù)獲取手段的豐富，水環(huán)境數(shù)值模擬的精度亟需進(jìn)一步提升。對(duì)于大型淺水湖泊來(lái)說(shuō)，研究域水平空間尺度遠(yuǎn)遠(yuǎn)大于垂向尺度，如太湖南北長(zhǎng)68.5km，東西平均寬34km，而湖泊平均水深只有1.9m，平面和垂向尺度相差萬(wàn)倍。根據(jù)2003年7月秦伯強(qiáng)[1]等在太湖烏龜山附近開展的湖流垂直結(jié)構(gòu)觀察發(fā)現(xiàn)，即使是太湖這樣最大水深不足3m的淺水湖泊，也存在垂向湖流分層的現(xiàn)象，且溶解氧、葉綠素及營(yíng)養(yǎng)鹽濃度等要素也與水深有著密切的響應(yīng)關(guān)系[2]，為了進(jìn)一步研究淺水湖泊流場(chǎng)內(nèi)部結(jié)構(gòu)，如補(bǔ)償流和上、下流矢的切變等，更精確的描述富營(yíng)養(yǎng)化的形成及轉(zhuǎn)化過(guò)程，勢(shì)必要更多的采用三維數(shù)值模擬，建立淺水湖泊生態(tài)水動(dòng)力及水質(zhì)模型，模擬水動(dòng)力及水質(zhì)要素遷移轉(zhuǎn)化過(guò)程。

本文從水動(dòng)力以及與垂向流場(chǎng)分布相關(guān)的水環(huán)境要素角度出發(fā)，總結(jié)、概括并分析了不同水質(zhì)要素垂向分布規(guī)律及與水深響應(yīng)關(guān)系，歸納了建模過(guò)程中關(guān)鍵參數(shù)的取值，并提出下一步淺水湖泊水環(huán)境模擬值得進(jìn)一步深入研究的方向。

2 淺水湖泊三維水環(huán)境數(shù)值模擬研究

水環(huán)境數(shù)值模擬研究始于1925年，至今已有92年發(fā)展歷程，美國(guó)、英國(guó)、韓國(guó)、澳大利亞、丹麥等國(guó)家曾先后發(fā)布了相應(yīng)的環(huán)境模型開發(fā)及應(yīng)用指南，并推薦了多種不同類型的數(shù)值模擬軟件。我國(guó)大陸地區(qū)的數(shù)值模擬研究始于20世紀(jì)70年代，在水環(huán)境數(shù)值模擬研究方面也已開展了大量工作[3]，其中針對(duì)湖泊的數(shù)值模擬已取得很多成果[4-5]。早期二維水環(huán)境數(shù)值模擬被廣泛用于淺水湖泊研究，但研究與水深有相關(guān)關(guān)系的要素時(shí)，則更多地采用三維數(shù)值模擬。

2.1 淺水湖泊三維風(fēng)生流結(jié)構(gòu)

湖泊流場(chǎng)主要有三種：風(fēng)生流、密度流和吞吐流，其中淺水湖泊中風(fēng)生流占主導(dǎo)地位。在定常風(fēng)作用下，淺水湖泊表層為風(fēng)生流，流向與風(fēng)向平行，中間為過(guò)渡層，底層流向相反，為補(bǔ)償流[6]。在大小相同，方向相反的定場(chǎng)風(fēng)作用下，湖流在水平方向上呈現(xiàn)流型相似，但環(huán)流方向相反的特征[7]，且風(fēng)速主要影響環(huán)流流速，風(fēng)向影響環(huán)流結(jié)構(gòu)[6, 8～10]，在垂向上，增水區(qū)域下沉運(yùn)動(dòng)強(qiáng)烈，減水區(qū)域存在較強(qiáng)的上升補(bǔ)償流，形成垂向封閉環(huán)流[11]。此外，當(dāng)風(fēng)生流主導(dǎo)湖泊流場(chǎng)時(shí)，存在臨界風(fēng)速值，超過(guò)臨界值的風(fēng)速會(huì)使流場(chǎng)結(jié)構(gòu)發(fā)生顯著改變[12]。關(guān)于鹽度對(duì)艾比湖[13]、美國(guó)切薩皮克灣[14]三維流場(chǎng)影響的研究結(jié)果表明，咸水湖與淡水湖水動(dòng)力空間分布相似，水體鹽度增加對(duì)湖泊環(huán)流結(jié)構(gòu)基本無(wú)影響。黃平[15]建立了基于交錯(cuò)網(wǎng)格差分方法的湖泊風(fēng)生流三維數(shù)值模型，成功運(yùn)用于武漢市墨水湖；曹慧江等[9]用三維數(shù)值模型計(jì)算了夏季滴水湖吞吐流、密度流和風(fēng)生流三者共同作用下的流場(chǎng)；湯露露等[16-17]建立了基于ECOMSED模式的太湖風(fēng)生流三維數(shù)值模型，其計(jì)算的太湖三維水流結(jié)構(gòu)與現(xiàn)有成果基本吻合[8, 18]，田勇[19-20]提出了能應(yīng)對(duì)實(shí)際工程計(jì)算需要的湖泊三維水動(dòng)力計(jì)算模型，該模型能夠有效模擬湖泊三維流場(chǎng)結(jié)構(gòu)以及風(fēng)生環(huán)流、補(bǔ)償流、吞吐流等復(fù)雜流速，揭示不同分層流場(chǎng)的流速差異；張利民、濮培民等建立了日本琵琶湖[21～23]水動(dòng)力模型，基本掌握了其水動(dòng)力過(guò)程及機(jī)制；Schmalz[24]建立了Okeechobee湖水位模型。

2.2 沉積物及水體污染物和藍(lán)藻空間分布

風(fēng)場(chǎng)可以影響水動(dòng)力過(guò)程，使湖泊局部水位產(chǎn)生變化，繼而改變水體運(yùn)動(dòng)速度和方向，影響各種物質(zhì)在湖泊內(nèi)的輸移擴(kuò)散，對(duì)湖泊的水質(zhì)及生態(tài)系統(tǒng)的結(jié)構(gòu)和功能產(chǎn)生重要影響[12]。水動(dòng)力引起的沉積物再懸浮是一次污染物大量釋放的過(guò)程，其污染物釋放量遠(yuǎn)遠(yuǎn)大于靜態(tài)條件下由孔隙水和上覆水體的濃度差而引起的釋放量[25]。研究表明太湖底泥懸浮的臨界風(fēng)速在5～6.5m/s之間，一旦超過(guò)臨界風(fēng)速，底部的泥沙會(huì)發(fā)生懸浮，引起內(nèi)源釋放[26-27]。三維數(shù)值模擬可以更精確的模擬水土界面的物質(zhì)交換過(guò)程及污染物到上覆水中的時(shí)空分布規(guī)律，Kopasakis K等[28]建立了淺水湖泊生態(tài)水動(dòng)力模擬系統(tǒng)，包括水動(dòng)力、沉積物和生態(tài)動(dòng)力學(xué)模塊；陳彧等[29](2010)運(yùn)用生態(tài)動(dòng)力學(xué)模型CAEDYM與三維水動(dòng)力模型ELCOM的耦合模型模擬了太湖三維流場(chǎng)下的水質(zhì)分布，并獲得一套適用于太湖的水質(zhì)模型參數(shù)；馮峰等[30]以淺水湖泊武漢東湖為研究對(duì)象，分析并比較了沉積物中碳、氮、磷等生源要素主要形態(tài)的垂向分布，結(jié)果表明pH值和氨氮都隨深度增加而有所增加最終趨于穩(wěn)定。陸健剛[31]等研究表明中水動(dòng)力和強(qiáng)水動(dòng)力下，水體中溶解態(tài)Cu、Zn、Pb含量從底層水體至表層水體呈對(duì)數(shù)增長(zhǎng)。

底泥懸浮引起的泥孔隙水中的營(yíng)養(yǎng)物質(zhì)釋放到上覆水中，水體營(yíng)養(yǎng)物質(zhì)濃度大大增加，為藻類提供營(yíng)養(yǎng)元素，進(jìn)而與湖泊富營(yíng)養(yǎng)化產(chǎn)生密切聯(lián)系。吳挺峰等[32]2008年9月在梅梁灣及貢湖灣水文水質(zhì)及氣象的觀測(cè)結(jié)果表明，觀測(cè)期間水動(dòng)力強(qiáng)度對(duì)太湖北部湖區(qū)葉綠素a濃度垂向分層及藍(lán)藻水華水平漂移均具有重要影響。王超等[33]建立了基于水動(dòng)力空間分布的淺水湖泊模型，研究了大型淺水湖泊中微囊藻素的運(yùn)輸過(guò)程，結(jié)果表明，微囊藻素垂向分布與風(fēng)力條件和對(duì)應(yīng)的水動(dòng)力條件有明顯的響應(yīng)關(guān)系，微風(fēng)條件(<3m/s)下，藍(lán)藻產(chǎn)生的區(qū)域從表層變?yōu)樗?.3m[34]。在水動(dòng)力滯緩水域，藍(lán)藻水華易在水表發(fā)生漂移堆積，在水動(dòng)力強(qiáng)度較大水域，強(qiáng)烈的垂向混合作用能使藍(lán)藻沿水深方向混合均勻，降低水華暴發(fā)風(fēng)險(xiǎn)。

風(fēng)的持續(xù)時(shí)間及其動(dòng)力性也是影響太湖藍(lán)藻水華時(shí)空間變化的關(guān)鍵因子之一[35]。Wu等[36]利用衛(wèi)星圖像結(jié)合實(shí)地抽樣調(diào)查研究發(fā)現(xiàn)，藍(lán)藻水華面積與風(fēng)速成負(fù)相關(guān)，較小風(fēng)速有利于藍(lán)藻水華形成，但短時(shí)間強(qiáng)風(fēng)(30分鐘平均風(fēng)速超過(guò)6m/s)引起的強(qiáng)烈混合作用則會(huì)增加大面積藍(lán)藻水華形成的機(jī)會(huì)，并且太湖地區(qū)4月和10月大風(fēng)頻繁可能正是藍(lán)藻水華大面積發(fā)生的重要原因。李未等[37](2016)建立了太湖三維水動(dòng)力—水質(zhì)耦合模型，對(duì)短期內(nèi)太湖湖泛易發(fā)水域及發(fā)生面積進(jìn)行預(yù)測(cè)，準(zhǔn)確率在80%以上。歐陽(yáng)瀟然等[38]利用FVCOM對(duì)太湖梅梁灣三維水溫及水體中溶解氧的時(shí)空分布進(jìn)行了模擬，相關(guān)研究標(biāo)明，水下光場(chǎng)對(duì)于混合藻類初級(jí)生產(chǎn)力具有明顯的影響，由于水體對(duì)于光照的衰減作用，太湖梅梁灣口出現(xiàn)光抑制的深度主要分布在水下0.1～0.3m區(qū)域[39]。

3 淺水湖泊數(shù)模關(guān)鍵參數(shù)研究

數(shù)值模擬過(guò)程中參數(shù)的率定至關(guān)重要，直接影響到模型的穩(wěn)定性及計(jì)算結(jié)果的精度。目前對(duì)數(shù)值模擬過(guò)程中的參數(shù)研究很多，但在規(guī)范化和標(biāo)準(zhǔn)化上仍有欠缺，亟需形成一套完整的適用于淺水湖泊水動(dòng)力模型的參數(shù)方案，本節(jié)內(nèi)容對(duì)前人研究過(guò)程中的相關(guān)參數(shù)取值進(jìn)行梳理和對(duì)比，以期為今后研究中的建模過(guò)程提供一定參考。

3.1 湖泊水動(dòng)力模型關(guān)鍵參數(shù)

水動(dòng)力過(guò)程主要涉及連續(xù)方程和動(dòng)量方程，連續(xù)方程如下：

?t(mζ)+?x(myHu)+?y(mxHv)+?z(mw)=0

(1)

(2)

其中H是水深，u、v分別是流速分量；w是垂向速度。mx和my是拉梅系數(shù)。動(dòng)量方程如下：

?t(mHu)+?x(myHuu)+?y(mxHvu)+?z(mwu)-(mf+v?xmy-u?ymx)Hv=-myH?x(gζ+p)-my(?xh-z?xH)?xp+?z(mH-1Av?zu)+Qu

(3)

?t(mHv)+?x(myHuv)+?y(mxHvv)+?z(mwv)+(mf+v?xmy-u?ymx)Hu=-mxH?y(gζ+p)-mx(?yh-z?yH)?zp+?z(mH-1Av?zv)+Qv

(4)

其中f是柯氏力系數(shù)，Av是垂直渦動(dòng)粘滯系數(shù)，Qu和Qv是動(dòng)量源匯項(xiàng)，g為重力加速度。

有關(guān)切應(yīng)力項(xiàng)的處理如下，水面風(fēng)切應(yīng)力：

(5)

(6)

式中wx、wy為10m高空處風(fēng)矢量在x、y 方向上的分量；Cf為風(fēng)的拖曳系數(shù)，由下式計(jì)算：

(7)

其中，ρa(bǔ)與ρw分別是空氣的密度與湖水的密度。

床面切應(yīng)力：

(8)

(9)

式中，cb為湖底床面的拖曳系數(shù)，由下式計(jì)算：

(10)

其中，Δbl是無(wú)量綱底部層的厚度，σ0是無(wú)量綱的底部粗糙度。水動(dòng)力數(shù)值模擬計(jì)算過(guò)程涉及多個(gè)參數(shù)輸入，選取其中對(duì)數(shù)值模擬計(jì)算結(jié)果影響較大的5個(gè)參數(shù)進(jìn)行分析。

(1)風(fēng)應(yīng)力拖曳系數(shù)

風(fēng)是大型淺水湖泊水流運(yùn)動(dòng)最主要的驅(qū)動(dòng)力之一，在風(fēng)生流的形成中風(fēng)應(yīng)力起決定性作用。風(fēng)應(yīng)力拖曳系數(shù)決定了大氣與湖泊間的動(dòng)量傳輸率。在選擇大型淺水湖泊水動(dòng)力模型參數(shù)時(shí)，要著重率定風(fēng)場(chǎng)參數(shù)，其中風(fēng)應(yīng)力拖曳系數(shù)對(duì)計(jì)算結(jié)果不確定性影響最高[40]。

在淺水湖泊水動(dòng)力模擬中，風(fēng)拖曳系數(shù)取值從1×10-6～3×10-6[41～43]各不相同，也有把拖曳系數(shù)作為隨風(fēng)速變化的變量，周婕等[44]指出，風(fēng)應(yīng)力拖曳系數(shù)為常數(shù)時(shí)，風(fēng)應(yīng)力和風(fēng)速平方呈線性關(guān)系，風(fēng)應(yīng)力拖曳系數(shù)取為和風(fēng)速有關(guān)的表達(dá)式時(shí)，兩者為非線性關(guān)系，模擬結(jié)果顯示，非線性關(guān)系的考慮相對(duì)更為合理。

(2)底部拖曳系數(shù)

底部拖曳系數(shù)Cb是空間變量，反映了水底糙度對(duì)能量損失的影響，文獻(xiàn)中的底部拖曳系數(shù)通常采用公式C=h1/6/n計(jì)算[45]，h可近似取計(jì)算區(qū)域內(nèi)的平均水深，n為曼寧系數(shù)(即底部糙率)，其取值直接影響底部拖曳系數(shù)(式10)，底部糙率是影響模型穩(wěn)定的重要參數(shù)，淺水湖泊底部通常比較平坦，故文獻(xiàn)中多采用經(jīng)驗(yàn)常數(shù)[46]，也有根據(jù)不同水生植物分布及下墊面不同采取分區(qū)取值方案[47]。

(3)渦粘系數(shù)

渦粘系數(shù)分為水平渦粘系數(shù)和垂向渦粘系數(shù)，部分研究者將渦粘系數(shù)取為常數(shù)，另有部分研究者將垂向渦粘系數(shù)取為與風(fēng)速有關(guān)的值，Koutitas[48]采用如下公式計(jì)算垂向渦粘系數(shù)

(11)

王惠中等[49]采用垂向變化的渦黏系數(shù)取值對(duì)太湖進(jìn)行數(shù)值模擬，結(jié)果表明，水平流速的梯度變化在變渦黏系數(shù)時(shí)要小于常渦黏系數(shù)，底部切應(yīng)力在常渦黏系數(shù)時(shí)要大于變渦黏系數(shù)，說(shuō)明取常渦黏系數(shù)實(shí)際上夸大了湖水流層間的切變作用和作用于湖底泥層的作用力，認(rèn)為在進(jìn)行太湖三維風(fēng)生流場(chǎng)模擬時(shí)不宜簡(jiǎn)單將垂向渦黏系數(shù)取常數(shù)值；Yongqi Wang[50]建立了Constance湖三維模型，采用Smagorinsky和2階Mellor-Yamada模型計(jì)算水平和垂向渦粘系數(shù)，計(jì)算結(jié)果表明，水平和垂向湍流渦動(dòng)粘滯系數(shù)對(duì)環(huán)流和溫度的垂向分布起著重要的作用，而不能取為常數(shù)。

(4)擴(kuò)散系數(shù)

擴(kuò)散系數(shù)主要是對(duì)水質(zhì)擴(kuò)散而言，其對(duì)水動(dòng)力模型計(jì)算結(jié)果影響較小，當(dāng)河流/水庫(kù)斷面形態(tài)單一，岸坡傾角大，水深增加快，擴(kuò)散系數(shù)呈現(xiàn)常數(shù)擴(kuò)散系數(shù)特征，可簡(jiǎn)化取為常數(shù)[7, 17, 40, 51]，對(duì)灘槽形復(fù)式斷面，深槽遠(yuǎn)離岸邊，岸坡傾角小，灘地水深淺，水深增加緩慢，擴(kuò)散系數(shù)呈現(xiàn)變擴(kuò)散系數(shù)特征[52]。

3.2 參數(shù)經(jīng)驗(yàn)值及經(jīng)驗(yàn)公式

本文查閱了前人在淺水湖泊數(shù)值模擬研究中的參數(shù)取值及經(jīng)驗(yàn)公式，現(xiàn)匯總供參考(見下表)。

表 淺水湖泊水動(dòng)力模型參數(shù)取值表Tab. Values of hydrodynamic model parameters in shallow lakes

4 展 望

隨著現(xiàn)場(chǎng)觀測(cè)手段的發(fā)展和計(jì)算機(jī)水平的提高，淺水湖泊數(shù)值模擬研究受到的限制越來(lái)越小，富營(yíng)養(yǎng)化是湖泊特別是淺水湖泊面臨的主要環(huán)境問(wèn)題，治理難度很大，其相關(guān)要素又與水深及水力條件空間分布有密切的響應(yīng)關(guān)系，對(duì)淺水湖泊進(jìn)行精細(xì)化的三維數(shù)值模擬變得越來(lái)越有必要。針對(duì)現(xiàn)有研究已取得的成果和存在的局限性，提出以下展望。

4.1 復(fù)雜邊界條件下的生態(tài)動(dòng)力學(xué)模型構(gòu)建

由于湖泊現(xiàn)場(chǎng)觀測(cè)難度大，大型淺水湖泊全湖區(qū)氣象及水文觀測(cè)數(shù)據(jù)極難獲得，通常采用假設(shè)的定常風(fēng)作為氣象邊界條件，未來(lái)可結(jié)合水利、氣象等管理部門，開發(fā)數(shù)據(jù)共享平臺(tái)，利用好現(xiàn)有的淺水湖泊現(xiàn)場(chǎng)觀測(cè)數(shù)據(jù)，開發(fā)基于大氣-水-底泥多邊界全過(guò)程生態(tài)模型，研究水-氣交界面動(dòng)力條件及碳循環(huán)過(guò)程，探索變化邊界條件下的氣候與富營(yíng)養(yǎng)化及沉積物再懸浮的循環(huán)和響應(yīng)機(jī)制，建立更符合天然情況的全過(guò)程生態(tài)動(dòng)力學(xué)數(shù)值模型。

4.2 多條件耦合的富營(yíng)養(yǎng)化微觀響應(yīng)機(jī)制

富營(yíng)養(yǎng)化產(chǎn)生需要在合適的光照、溫度、營(yíng)養(yǎng)鹽濃度及微生物生長(zhǎng)所需水力條件，水體中的氮磷營(yíng)養(yǎng)鹽濃度對(duì)微生物的降解功能有顯著影響，營(yíng)養(yǎng)鹽分布也與微生物分布有密切聯(lián)系，水中植物會(huì)影響水流結(jié)構(gòu)分布，需要進(jìn)一步開展上述要素之間的協(xié)同作用和遷移轉(zhuǎn)化規(guī)律，研究多條件耦合下的富營(yíng)養(yǎng)化產(chǎn)生微觀機(jī)制。

4.3 數(shù)值模擬精細(xì)化及參數(shù)規(guī)范化研究

現(xiàn)有數(shù)學(xué)模型均基于一定假設(shè)條件下，且模擬過(guò)程中各項(xiàng)參數(shù)取值差異較大，尚未針對(duì)淺水湖泊形成一套規(guī)范化的參數(shù)取值方案，可結(jié)合前人研究，通過(guò)物理水槽試驗(yàn)，得出參數(shù)初始值或公式，并通過(guò)湖泊原型觀測(cè)數(shù)據(jù)進(jìn)行驗(yàn)證，最終得到適用于淺水湖泊的規(guī)范化參數(shù)取值方案。