李一平,施媛媛,姜 龍,朱向宇,龔 然
(1.河海大學環(huán)境學院,江蘇 南京 210098; 2.清水源(上海)環(huán)保科技有限公司,上海 201107;3.南京工程學院環(huán)境工程學院,江蘇 南京 211167)
水環(huán)境保護是國家水安全保障的重中之重,掌握新形勢下水環(huán)境演變過程是開展水環(huán)境保護工作的關(guān)鍵[1]。研究水環(huán)境問題的手段主要包括野外觀測、室內(nèi)試驗和數(shù)值模擬,其中水環(huán)境數(shù)學模型扮演著越來越重要的角色。地表水環(huán)境數(shù)學模型(surface water environment numerical models,SWENM)主要分為水動力學模型、水質(zhì)模型和水生態(tài)模型,基本原理是將氣象條件、水動力條件、邊界條件等因素進行定量化約束,通過求解方程組,獲得污染物的時空分布特征及遷移轉(zhuǎn)化規(guī)律,分析和判別各環(huán)境因子間的相互關(guān)系,實現(xiàn)模擬與預測等功能。
目前,已研發(fā)的地表水環(huán)境數(shù)學模型有數(shù)百種,算法及應用的差異給水環(huán)境決策帶來了一定的困難?,F(xiàn)有綜述論文對3類模型的介紹不夠全面,水文領(lǐng)域更關(guān)注水動力模型,環(huán)境領(lǐng)域又更注重水質(zhì)模型,水生態(tài)模型因其復雜性及涵蓋面廣難以進行整合歸納。此外,關(guān)于模型模擬精度影響因素的綜合探討比較欠缺。因此,歸納比較各種模型的適用范圍、強化模型敏感性和不確定分析、辨析模型精度的影響因素等均有利于地表水環(huán)境數(shù)學模型的精準優(yōu)化發(fā)展。鑒于傳統(tǒng)模型應用的時空局限性,“3S”、云計算、人工智能等新技術(shù)的引入以及模型綜合化和法規(guī)化,將成為水環(huán)境模擬與預測領(lǐng)域未來的研究熱點。
近代水動力學起步于Navier-Stokes方程(簡稱N-S方程),隨后,Saint-Venant于1871年首次提出了計算一維河道及河網(wǎng)水流的一維水流運動基本方程——Saint-Venant方程。由此,水動力學模型逐步應用于河流、河網(wǎng);20世紀60—70年代,隨著Saint-Venant方程的廣泛應用,各種求解方法的數(shù)值穩(wěn)定性和精度問題得到深入研究。平原地區(qū)河道交錯,流向流速易受潮汐、水利調(diào)度等影響而催生了涉及面更廣的河網(wǎng)模型[2]。在湖庫方面,Hesen于1956年最早提出淺水平面二維水動力學模型。之后的模型側(cè)重于水流運動,包括風生環(huán)流[3-4]和吞吐流,以及在深水湖庫還因溫度分層而存在垂向密度流[5-6]。
隨著計算機技術(shù)引入水動力學,計算水動力學在計算方法、網(wǎng)絡技術(shù)、紊流模型、大渦模擬等方面的新理論和新方法也為水動力學模型領(lǐng)域提供了有效借鑒[7]。直接數(shù)值模擬(direct numerical simulation,DNS)、Reynolds平均模型(Reynolds average Navier-Stokes,RANS)和大渦模擬(large eddy simulation,LES)等方法得以應用,進一步修正了N-S方程,典型的紊流模型有k-ε模型和v2f模型等[8-9]。此外,水動力學模型在洋流領(lǐng)域的應用研究推動了紊流機理的深入探索,洋流模型由最初的普林斯頓海洋模型(Princeton ocean model,POM)發(fā)展至淺海三維水動力模型(3D estuarine coastaland ocean model,ECOM),直至如今更為完善的有限體積海岸海洋模型(finite volume coastal ocean model,F(xiàn)VCOM)。在河口海岸區(qū)域,水動力學模型解決了潮汐作用下的紊流混合[10]、因鹽度差而形成的密度流以及因徑流和鹽度密度流產(chǎn)生的入??跍魡栴}[11]。
當前,水動力學模型與其他模型(流域水文模型、波浪模型和泥沙模型等)的耦合應用實現(xiàn)了空間多維度的模擬。例如,外夾河流域的耦合模型以水文模型計算的流量作為水動力學模型的邊界條件,成功實現(xiàn)了丘陵、平原混合地區(qū)的洪水預報[12];太湖透明度模型[13]基于淺水波浪數(shù)值模型(simulating waves nearshore,SWAN),耦合了波浪模塊和湖流三維模型,動態(tài)模擬了太湖波浪和湖流的生消過程;三維水動力模型也常常應用于探究泥沙輸移規(guī)律[14]及魚道模擬[15]。
自1925年Streeter和Phelps建立了BOD-DO耦合模型(S-P模型)以來,水質(zhì)模型已經(jīng)發(fā)展了90多年。水質(zhì)模型基于水動力學模型,根據(jù)物質(zhì)守恒原理概化污染物在水中發(fā)生的物理、化學、生物化學變化過程。主要分為4個發(fā)展階段[16]:
第一階段(1925—1965年):簡單的氧平衡模型階段。集中研究氧平衡,部分涉及非耗氧物質(zhì),不斷修正S-P模型。
第二階段(1966—1985年):水質(zhì)模型迅速發(fā)展階段??紤]污染物不同形態(tài)的影響機制,20世紀70年代初,美國的研發(fā)機構(gòu)開始推出QUAL-Ⅰ、WASP(water quality analysis simulation program modeling system)等綜合水質(zhì)模型軟件,后續(xù)列入其他污染源、底泥、邊界的作用,研發(fā)的QUAL系列模型應用于河流綜合水質(zhì)規(guī)劃管理[17]。
第三階段(1986—1998年):水質(zhì)模型深入研究、廣泛應用階段。更多復雜的水質(zhì)過程被納入模型系統(tǒng),例如,考慮沉積物“土-水”界面動態(tài)過程的沉積成巖模型[18];模擬對象涵蓋河流、湖泊(水庫)、河口、海岸帶等。
第四階段(1999年至今):水質(zhì)模型集成化、設計人性化階段。MIKE、EFDC(environmental fluid dynamics code)、DELFT-3D等模型集成水動力學、水質(zhì)、泥沙、生態(tài)等模塊,并提供網(wǎng)格生成和前后處理工具。新興技術(shù)也逐步被引入水質(zhì)模型,人工智能提高了水質(zhì)模型的預測水平[19],遺傳算法、模擬退火算法強化了參數(shù)識別[20],神經(jīng)網(wǎng)絡明晰了河網(wǎng)物理結(jié)構(gòu)[21]。
水生態(tài)模型是描述水生生態(tài)系統(tǒng)中生物個體或種群間的內(nèi)在變化機制,及構(gòu)建水文、水質(zhì)、氣象等因素連接的復雜模型,主要用于研究水體富營養(yǎng)化、生物富集及水域系統(tǒng)食物網(wǎng)[25]。
20世紀70年代初期誕生的簡單總磷模型成為水生態(tài)模型的基石;20世紀80年代,美國和日本開發(fā)了第一批三維生態(tài)數(shù)學模型[25]?,F(xiàn)代水生態(tài)模型考慮了自然界中多因素相互作用及時空變化。例如,Delft3D BLOOM/GEM模型能成功模擬4種不同海域的水生系統(tǒng)[26];三維ERSEM(european regional seas ecosystem model)能夠用于研究紅海的水動力和生化動力[27]。
近年來,水生植物、魚類遷移及生物棲息地等模塊得到不斷開發(fā)。江志超等[28]建立了硅藻中肋骨條藻與氮、磷關(guān)系的非線性動力學模型,發(fā)現(xiàn)光衰變率和營養(yǎng)鹽是影響中肋骨條藻赤潮生消過程的關(guān)鍵因子。RIVER2D模型已被用于識別在研究范圍內(nèi)最小水深和寬度處成年大鱗大馬哈魚[29]?;贖ABITAT模型建立的生物棲息地評價模型,能夠反映瀘沽湖水質(zhì)變化對寧蒗裂腹魚的影響[30]。此外“水生態(tài)足跡”概念的提出衍生出了水生態(tài)足跡計算模型、水生態(tài)承載力計算模型和水生態(tài)赤字或盈余計算模型[31]。
當前國內(nèi)水環(huán)境模型研究主要關(guān)注水動力-水質(zhì)-水生態(tài)耦合模型的應用。例如,太湖富營養(yǎng)化機理模型耦合了太湖三維風生湖流模型、垂向平均的二維水質(zhì)模型和富營養(yǎng)化模型,考慮了水溫、總氮、總磷和太陽輻射等因子對藻類生長的影響,模擬了藻類生消過程以及其隨風生流遷移的規(guī)律[32-33];王生愿等[34]以湖泊水動力水質(zhì)響應機制為研究背景,構(gòu)建了水-生態(tài)-底泥耦合的湖泊水動力水生態(tài)模型。
近30年來,先后涌現(xiàn)出許多高品質(zhì)的地表水環(huán)境數(shù)學模型,例如WASP模型素有“萬能水質(zhì)模型”之稱;還有近年應用廣泛的EFDC(environmental fluid dynamics code)模型[35],集成了水動力、水質(zhì)、風浪、泥沙、重金屬及有毒物質(zhì)、沉積成巖和水生植物等模塊。近些年國內(nèi)也開發(fā)出了諸如CJK3D、IWIND等商用軟件。表1列出了國內(nèi)外常用的18個地表水環(huán)境數(shù)學模型。
模型不確定性和敏感性分析是數(shù)值建模研究的熱點。不確定性分析是將模型輸出中的不確定性進行量化評定[36],而敏感性分析是研究模型輸入因素和輸出變化的響應關(guān)系。
水質(zhì)模型不確定性的研究始于20世紀70年代,O’Neill等[37]指出單純尋找模型最優(yōu)化參數(shù)沒有意義,需要確定參數(shù)的分布;之后,水質(zhì)模型的不確定性得到了系統(tǒng)性歸納[38],采用不同的判定標準能夠確定水質(zhì)模型不確定性來源及分類[39];生態(tài)模型和統(tǒng)計模型開啟了水環(huán)境數(shù)學模型不確定性分析的時代[40];各種不確定性分析方法基本原理的提出,為模型參數(shù)排序及決策評估提供了有效依據(jù)[36]。
模擬結(jié)果的不確定性主要來自參數(shù)、輸入數(shù)據(jù)和模型結(jié)構(gòu)不確定性3個方面,其中參數(shù)不確定性指參數(shù)估計存在誤差,輸入數(shù)據(jù)不確定性指模型邊界條件和初始條件的不確定性,而模型結(jié)構(gòu)不確定性是由于人類對復雜環(huán)境系統(tǒng)認識的局限性,在系統(tǒng)建模過程中常常對一些現(xiàn)象和變化過程進行抽象和概化。
常用的不確定性分析的數(shù)學表達方法主要有區(qū)間數(shù)學法、模糊理論法以及概率分析法。區(qū)間數(shù)學法用于計算測量和參數(shù)估值誤差引起的不確定性;模糊理論法解決具有模糊性的系統(tǒng)不確定問題,但難以實現(xiàn)定量評估;概率分析法常用于描述物理系統(tǒng)的不確定性,根據(jù)模型輸入的概率分布來確定模型輸出的概率分布,最終以概率分布的形式來表達不確定性,如蒙特卡羅(Monte Carlo)法[41]、拉丁超立方抽樣法(Latin hypercube sampling,LHS)[42]、普適似然不確定估計方法(generalized likelihood uncertainty estimation,GLUE)[43]以及單純多邊形進化算法(shuffled complex evolution algorithm,SCE-UA)[44]等。
表1 常用地表水環(huán)境數(shù)學模型
敏感性分析方法主要分為局部分析方法和全局分析方法。常采用的局部分析方法是檢驗單個參數(shù)的變化對模型結(jié)果的影響程度(one factor at a time,OAT)。局部分析方法簡單、易于實施,但響應結(jié)果較片面,無法解決“異參同效”問題。全局分析方法克服了局部分析方法的缺點,能夠反映整體參數(shù)組合對結(jié)果輸出不確定性的影響。全局敏感性分析方法有很多[45]:①基于回歸或相關(guān)分析技術(shù)的方法,如多元回歸法、響應曲面方法(response surface methodology,RSM);②全局篩選法,如LH-OAT方法、Morris方法;③基于方差理論的方法,如傅里葉振幅敏感性檢驗法(Fourier amplitude sensitivity test,F(xiàn)AST)、Sobol方法和擴展傅里葉振幅敏感性檢驗法(extend FAST);④因子設計實驗[46];⑤摩爾斯分析法[47];⑥取樣分析法[48]等。
敏感性分析常用的蒙特卡羅法屬于取樣分析法,能簡單有效地評價多個參數(shù)對模型輸出結(jié)果不確定性的貢獻。LHS克服了蒙特卡羅法計算成本高的缺點,抽取的樣本能更精確地反映輸入概率函數(shù)的分布,不僅高度控制抽樣值,又為它們留有變化的余地?;贚HS的思想,運行模擬的次數(shù)由輸入變量數(shù)決定,最少可為隨機變量數(shù)的1.5倍,一般為幾百次[49]。
當前國內(nèi)外已有對模型敏感性和不確定性研究的成功案例。李一平等[50-51]利用EFDC模型和LHS抽樣方法分析了水動力模塊中的5個重要參數(shù)以及4個重要的外部輸入條件對太湖水位和流場分布的敏感性和不確定性,并借助原位觀測簡化了模型參數(shù)的率定驗證;Gong等[52]基于SWAT(soil and water assessment tool)模型采用GLUE不確定性分析方法分析模型參數(shù)對模型不確定性的影響,其中針對不同的研究目標需要選取不同的似然函數(shù)[53];Reder等[54]采用LHS和全局敏感性分析法(global sensitivity analysis,GSA)在42個參數(shù)中篩選出4個最敏感的參數(shù)以提高模型性能。
除了傳統(tǒng)分析方法以外,基于熵的不確定性敏感性分析方法也被提出[55],并用以識別不確定變量對隨機變量和模糊變量組合系統(tǒng)的影響[56]。全局靈敏度測算方法(sobol)是一種兩階段不確定性量化方法,簡化了不確定度量化計算過程[57]。風場、邊界及底部地形的空間變化導致研究區(qū)域內(nèi)的參數(shù)分布也有所差異,因此結(jié)合GIS-Lab衍生的參數(shù)空間不確定分析將成為新的研究方向。
參數(shù)不確定性和敏感性分析已成為模型構(gòu)建的必要工作,通過分析可表征輸出結(jié)果的不確定性對輸入?yún)?shù)不確定性的依賴程度。對于特定的模擬目標,篩選出敏感參數(shù)、參數(shù)合理分布范圍以及各個參數(shù)對模型結(jié)果不確定性的貢獻率,可大幅度減少模型后續(xù)率定驗證工作,同時可指導關(guān)鍵參數(shù)的監(jiān)測工作。
除由參數(shù)不確定性引起的模擬誤差外,影響模型模擬精度的因素還包括模型類型及選擇、網(wǎng)格種類、垂向坐標系統(tǒng)、方程離散方法、初始條件及邊界條件等。
選擇合適的模型是建模前應明確的重要環(huán)節(jié),需要參考模擬對象的特性、精度要求和計算效率等。當模擬對象具有干濕交替或洪泛區(qū)特性時,應選用較為精細的水動力模型[58];若模擬深水湖庫水質(zhì),需要選用垂向二維或三維模型,考察垂向溫度和水質(zhì)變化[59];若模擬重金屬或有毒物質(zhì)遷移轉(zhuǎn)換過程,關(guān)鍵是構(gòu)建并校驗泥沙模塊,其余影響因素和模塊可進行一定程度的概化處理[60]。由于復雜的物理、化學、生物過程無法在模型中詳述,所以模型選擇不恰當可能會浪費計算時間,且不一定能達到預期目標。因此,模型的選擇并非越復雜越好,在明確模擬目標后,選擇合適的模型即可,這個環(huán)節(jié)需要一定的實踐經(jīng)驗。
網(wǎng)格種類及坐標系統(tǒng)選擇也直接影響模擬精度。常見的平面網(wǎng)格有矩形網(wǎng)格、正交曲線網(wǎng)格、三角形網(wǎng)格及四叉樹網(wǎng)格。矩形網(wǎng)格基于直角坐標系便于組織數(shù)據(jù)結(jié)構(gòu),計算效率高,但不適合處理復雜的邊界且不易調(diào)控網(wǎng)格密度。正交曲線網(wǎng)格是一種基于曲線坐標系統(tǒng)的有結(jié)構(gòu)網(wǎng)格,可以適應不規(guī)則邊界,但處理過于復雜的邊界時效果不佳。三角形網(wǎng)格利于研究復雜地形和邊界問題,易于控制網(wǎng)格密度,但計算效率較低。具有樹狀結(jié)構(gòu)的四叉樹網(wǎng)格能高度擬合復雜的自然水體,易實現(xiàn)水動力及物質(zhì)輸運數(shù)值模擬[61]。在邊界和地形較復雜的位置宜采用三角形網(wǎng)格,在計算域內(nèi)部和地形變化不大的地方宜采用矩形網(wǎng)格或者正交曲線網(wǎng)格。
垂向坐標系統(tǒng)一般分為平面(z)坐標、等密度(ρ)坐標和地形擬合(σ)坐標,分別對應不同的網(wǎng)格。z坐標模式方程簡單,易于數(shù)值離散,適用于具有準水平運動特點的水體,但不易處理底部邊界,在淺水區(qū)域難以滿足必要的垂直分辨率。ρ坐標模式常用于密度流模擬,但在混合層、非層化水體內(nèi)和底部邊界層的分辨率較低。σ坐標模式實現(xiàn)了垂直相對分層,能夠有效擬合底部地形,但斜壓梯度力會有較大截斷誤差。對于地形復雜的水域可采用混合坐標模式,例如σ-z坐標系統(tǒng)適用于局部陡峭深水區(qū),能夠降低水平壓力梯度誤差的影響。在實際應用中,應根據(jù)岸線形狀、水下地形數(shù)據(jù)和模擬精度的要求選擇合適的網(wǎng)格及垂向坐標系統(tǒng)。
精細的分辨率能降低模擬誤差,但會增加不必要的模擬成本。若部分區(qū)域的模擬精度要求較高,可進行局部網(wǎng)格加密處理,既保證擬合計算的合理性,又提高運行效率。
為確保計算的穩(wěn)定性和收斂性,時間步長應該足夠小,通常須將時間步長減小到幾分或幾秒,與水動力過程模擬所需的時間步長差不多,可重現(xiàn)泥沙輸送等水質(zhì)動力學過程。目前基于有限差分法、有限體積法和有限單元法的高效數(shù)值離散計算方法、并行計算方法、集群計算方法等同時滿足了精度和效度的要求[62]。有限差分法根據(jù)時間和空間步長對定解區(qū)域進行網(wǎng)格劃分,用差商代替導數(shù),求解方法簡單但不易處理復雜邊界問題。有限體積法可以根據(jù)實際問題的物理特點對任意形狀網(wǎng)格體進行積分,且不會影響計算精度和守恒性,但是對于質(zhì)量差的網(wǎng)格,離散的過程會產(chǎn)生更大的誤差。有限單元法根據(jù)實際問題的物理特點對求解區(qū)域進行單元剖分,能夠滿足一定的精度要求,但對于二維和三維問題需要建立許多人為的節(jié)點,且求解精度過于依賴網(wǎng)格劃分,適合分析連續(xù)變形問題。
初始條件包括初始水位、流速、溫度等,初始條件設置是否準確,對不同的水體影響不同。模擬湖庫水質(zhì)時,若只依賴邊界條件驅(qū)動,運行至合理的初始狀態(tài)需要一定的時間且對模擬的結(jié)果影響較大。因此,模擬湖庫水質(zhì)時,如沒有準確的實測值作為初始條件,則通常將模型運行一段時間后的結(jié)果作為初始態(tài),稱為“預熱”。對于河流、河口等水動力過程較劇烈的水體,邊界條件的驅(qū)動將很快覆蓋初始條件,初始條件的影響很小。
邊界條件包括大氣邊界、出入流邊界、開邊界的作用力、水工建筑物、取水退水邊界,其中水質(zhì)模塊的邊界條件還需考慮各種出入水體邊界的水質(zhì)變量、大氣干濕沉降、農(nóng)田面源污染、地表徑流、內(nèi)源釋放、地下水等。在確定污染源的過程中,大氣干濕沉降對水域整體的貢獻不可忽略,隨降雨進入水體的污染物質(zhì)往往是水質(zhì)惡化的關(guān)鍵因素[63]。邊界條件還可分為垂向和水平邊界條件,比如氣溫和風速作為模型垂向邊界條件,雖然不參與直接計算模擬,但它們影響了潮流、混合和熱傳輸?shù)人畡恿^程。合適的邊界條件可以有效避免模擬誤差。
目前,地表水環(huán)境數(shù)學模型正朝著系統(tǒng)化、綜合化、法規(guī)化方向發(fā)展,模型涉及的要素越來越復雜,與許多新興技術(shù)的結(jié)合,使得地表水環(huán)境數(shù)學模型的發(fā)展充滿了機遇與挑戰(zhàn)。
a. 模型系統(tǒng)的系統(tǒng)化、綜合化和平臺化。水環(huán)境模型逐步涉及水文、水動力、生物、地理等多領(lǐng)域,以單獨模塊的形式集成一體。模擬元素和模擬情景的增加拓展了模型的研究領(lǐng)域,比如藥品及個人護理用品(PPCPs)、微塑料等新型污染物將成為新的模擬對象,藻類競爭、熱量平衡等過程也逐步引入水生態(tài)模型。適用于大型流域的三維水生態(tài)模型開發(fā),以及利用統(tǒng)一平臺構(gòu)建促使數(shù)值模擬成為流域控制與規(guī)劃決策的支撐技術(shù),將成為未來重要的發(fā)展趨勢。這些將會是今后地表水環(huán)境數(shù)學模型研究的熱點。
b. 與新興技術(shù)的結(jié)合。大數(shù)據(jù)平臺的參與能夠有效解決邊界條件、初始條件輸入缺失難題;超算、云計算以及人工智能算法的應用將會大幅度提升地表水環(huán)境數(shù)學模型的計算效率和精度,比如遺傳算法、模擬退火算法能夠強化參數(shù)識別;VR技術(shù)豐富了模型結(jié)果的輸出展現(xiàn)方式,實現(xiàn)了人性化的設計理念。在與物聯(lián)網(wǎng)、“互聯(lián)網(wǎng)+”、云技術(shù)的碰撞下,“3S”技術(shù)與水環(huán)境模擬集成模塊的開發(fā)促進了環(huán)境因素之間的相互作用和動態(tài)變化的確定性分析,基于遙感技術(shù)的水質(zhì)反演模型促進了“天地一體化”水環(huán)境監(jiān)測系統(tǒng)的建立,能夠?qū)崿F(xiàn)數(shù)字預警、識別黑臭水體、考察海域污染,體現(xiàn)實時性、大尺度、高速度、動態(tài)性等優(yōu)勢。在GIS-Lab技術(shù)的引入結(jié)合下,風場、地形、床層等空間因素的考慮會更加全面,不確定性分析正從全局概化向區(qū)域性精細化發(fā)展,模擬過程更加真實,結(jié)果更加可靠。
c. 法規(guī)化趨勢。基于GIS技術(shù),模型庫管理系統(tǒng)的建立促使水環(huán)境模型走上法規(guī)化管理道路。目前,美國國家環(huán)境保護局(EPA)已將97種地表水環(huán)境數(shù)學模型列入模型信息庫,澳大利亞政府也針對模型選擇、參數(shù)率定、敏感性分析等給出了系統(tǒng)推薦,我國生態(tài)環(huán)境部也正式發(fā)布了HJ2.3—2018《環(huán)境影響評價技術(shù)導則 地表水環(huán)境》,并推薦了適用于河流、湖泊、河口和海洋的數(shù)值模型。今后的發(fā)展將基于已有的模型信息庫進行拓展豐富,有望實現(xiàn)模型運用流程的規(guī)范化操作、指導、監(jiān)督。