地表水環境數學模型研究進展

李一平，施媛媛，姜 龍，朱向宇，龔 然

(1.河海大學環境學院，江蘇 南京 210098； 2.清水源(上海)環?？萍加邢薰?，上海 201107；3.南京工程學院環境工程學院，江蘇 南京 211167)

水環境保護是國家水安全保障的重中之重，掌握新形勢下水環境演變過程是開展水環境保護工作的關鍵[1]。研究水環境問題的手段主要包括野外觀測、室內試驗和數值模擬，其中水環境數學模型扮演著越來越重要的角色。地表水環境數學模型(surface water environment numerical models，SWENM)主要分為水動力學模型、水質模型和水生態模型，基本原理是將氣象條件、水動力條件、邊界條件等因素進行定量化約束，通過求解方程組，獲得污染物的時空分布特征及遷移轉化規律，分析和判別各環境因子間的相互關系，實現模擬與預測等功能。

目前，已研發的地表水環境數學模型有數百種，算法及應用的差異給水環境決策帶來了一定的困難。現有綜述論文對3類模型的介紹不夠全面，水文領域更關注水動力模型，環境領域又更注重水質模型，水生態模型因其復雜性及涵蓋面廣難以進行整合歸納。此外，關于模型模擬精度影響因素的綜合探討比較欠缺。因此，歸納比較各種模型的適用范圍、強化模型敏感性和不確定分析、辨析模型精度的影響因素等均有利于地表水環境數學模型的精準優化發展。鑒于傳統模型應用的時空局限性，“3S”、云計算、人工智能等新技術的引入以及模型綜合化和法規化，將成為水環境模擬與預測領域未來的研究熱點。

1 地表水環境數學模型的研究現狀

1.1 水動力學模型

近代水動力學起步于Navier-Stokes方程(簡稱N-S方程)，隨后，Saint-Venant于1871年首次提出了計算一維河道及河網水流的一維水流運動基本方程——Saint-Venant方程。由此，水動力學模型逐步應用于河流、河網；20世紀60—70年代，隨著Saint-Venant方程的廣泛應用，各種求解方法的數值穩定性和精度問題得到深入研究。平原地區河道交錯，流向流速易受潮汐、水利調度等影響而催生了涉及面更廣的河網模型[2]。在湖庫方面，Hesen于1956年最早提出淺水平面二維水動力學模型。之后的模型側重于水流運動，包括風生環流[3-4]和吞吐流，以及在深水湖庫還因溫度分層而存在垂向密度流[5-6]。

隨著計算機技術引入水動力學，計算水動力學在計算方法、網絡技術、紊流模型、大渦模擬等方面的新理論和新方法也為水動力學模型領域提供了有效借鑒[7]。直接數值模擬(direct numerical simulation,DNS)、Reynolds平均模型(Reynolds average Navier-Stokes,RANS)和大渦模擬(large eddy simulation，LES)等方法得以應用，進一步修正了N-S方程，典型的紊流模型有k-ε模型和v2f模型等[8-9]。此外，水動力學模型在洋流領域的應用研究推動了紊流機理的深入探索，洋流模型由最初的普林斯頓海洋模型(Princeton ocean model，POM)發展至淺海三維水動力模型(3D estuarine coastaland ocean model，ECOM)，直至如今更為完善的有限體積海岸海洋模型(finite volume coastal ocean model，FVCOM)。在河口海岸區域，水動力學模型解決了潮汐作用下的紊流混合[10]、因鹽度差而形成的密度流以及因徑流和鹽度密度流產生的入海口滯留問題[11]。

當前，水動力學模型與其他模型(流域水文模型、波浪模型和泥沙模型等)的耦合應用實現了空間多維度的模擬。例如，外夾河流域的耦合模型以水文模型計算的流量作為水動力學模型的邊界條件，成功實現了丘陵、平原混合地區的洪水預報[12]；太湖透明度模型[13]基于淺水波浪數值模型(simulating waves nearshore，SWAN)，耦合了波浪模塊和湖流三維模型，動態模擬了太湖波浪和湖流的生消過程；三維水動力模型也常常應用于探究泥沙輸移規律[14]及魚道模擬[15]。

1.2 水質模型

自1925年Streeter和Phelps建立了BOD-DO耦合模型(S-P模型)以來，水質模型已經發展了90多年。水質模型基于水動力學模型，根據物質守恒原理概化污染物在水中發生的物理、化學、生物化學變化過程。主要分為4個發展階段[16]：

第一階段(1925—1965年)：簡單的氧平衡模型階段。集中研究氧平衡，部分涉及非耗氧物質，不斷修正S-P模型。

第二階段(1966—1985年)：水質模型迅速發展階段。考慮污染物不同形態的影響機制，20世紀70年代初，美國的研發機構開始推出QUAL-Ⅰ、WASP(water quality analysis simulation program modeling system)等綜合水質模型軟件，后續列入其他污染源、底泥、邊界的作用，研發的QUAL系列模型應用于河流綜合水質規劃管理[17]。

第三階段(1986—1998年)：水質模型深入研究、廣泛應用階段。更多復雜的水質過程被納入模型系統，例如，考慮沉積物“土-水”界面動態過程的沉積成巖模型[18]；模擬對象涵蓋河流、湖泊(水庫)、河口、海岸帶等。

第四階段(1999年至今)：水質模型集成化、設計人性化階段。MIKE、EFDC(environmental fluid dynamics code)、DELFT-3D等模型集成水動力學、水質、泥沙、生態等模塊，并提供網格生成和前后處理工具。新興技術也逐步被引入水質模型，人工智能提高了水質模型的預測水平[19]，遺傳算法、模擬退火算法強化了參數識別[20]，神經網絡明晰了河網物理結構[21]。

1.3 水生態模型

水生態模型是描述水生生態系統中生物個體或種群間的內在變化機制，及構建水文、水質、氣象等因素連接的復雜模型，主要用于研究水體富營養化、生物富集及水域系統食物網[25]。

20世紀70年代初期誕生的簡單總磷模型成為水生態模型的基石；20世紀80年代，美國和日本開發了第一批三維生態數學模型[25]。現代水生態模型考慮了自然界中多因素相互作用及時空變化。例如，Delft3D BLOOM/GEM模型能成功模擬4種不同海域的水生系統[26]；三維ERSEM(european regional seas ecosystem model)能夠用于研究紅海的水動力和生化動力[27]。

近年來，水生植物、魚類遷移及生物棲息地等模塊得到不斷開發。江志超等[28]建立了硅藻中肋骨條藻與氮、磷關系的非線性動力學模型，發現光衰變率和營養鹽是影響中肋骨條藻赤潮生消過程的關鍵因子。RIVER2D模型已被用于識別在研究范圍內最小水深和寬度處成年大鱗大馬哈魚[29]?；贖ABITAT模型建立的生物棲息地評價模型，能夠反映瀘沽湖水質變化對寧蒗裂腹魚的影響[30]。此外“水生態足跡”概念的提出衍生出了水生態足跡計算模型、水生態承載力計算模型和水生態赤字或盈余計算模型[31]。

當前國內水環境模型研究主要關注水動力-水質-水生態耦合模型的應用。例如，太湖富營養化機理模型耦合了太湖三維風生湖流模型、垂向平均的二維水質模型和富營養化模型，考慮了水溫、總氮、總磷和太陽輻射等因子對藻類生長的影響，模擬了藻類生消過程以及其隨風生流遷移的規律[32-33]；王生愿等[34]以湖泊水動力水質響應機制為研究背景，構建了水-生態-底泥耦合的湖泊水動力水生態模型。

2 常用地表水環境數學模型

近30年來，先后涌現出許多高品質的地表水環境數學模型，例如WASP模型素有“萬能水質模型”之稱；還有近年應用廣泛的EFDC(environmental fluid dynamics code)模型[35]，集成了水動力、水質、風浪、泥沙、重金屬及有毒物質、沉積成巖和水生植物等模塊。近些年國內也開發出了諸如CJK3D、IWIND等商用軟件。表1列出了國內外常用的18個地表水環境數學模型。

3 地表水環境數學模型的不確定性和敏感性

模型不確定性和敏感性分析是數值建模研究的熱點。不確定性分析是將模型輸出中的不確定性進行量化評定[36]，而敏感性分析是研究模型輸入因素和輸出變化的響應關系。

水質模型不確定性的研究始于20世紀70年代，O’Neill等[37]指出單純尋找模型最優化參數沒有意義，需要確定參數的分布；之后，水質模型的不確定性得到了系統性歸納[38]，采用不同的判定標準能夠確定水質模型不確定性來源及分類[39]；生態模型和統計模型開啟了水環境數學模型不確定性分析的時代[40]；各種不確定性分析方法基本原理的提出，為模型參數排序及決策評估提供了有效依據[36]。

模擬結果的不確定性主要來自參數、輸入數據和模型結構不確定性3個方面，其中參數不確定性指參數估計存在誤差，輸入數據不確定性指模型邊界條件和初始條件的不確定性，而模型結構不確定性是由于人類對復雜環境系統認識的局限性，在系統建模過程中常常對一些現象和變化過程進行抽象和概化。

常用的不確定性分析的數學表達方法主要有區間數學法、模糊理論法以及概率分析法。區間數學法用于計算測量和參數估值誤差引起的不確定性；模糊理論法解決具有模糊性的系統不確定問題，但難以實現定量評估；概率分析法常用于描述物理系統的不確定性，根據模型輸入的概率分布來確定模型輸出的概率分布，最終以概率分布的形式來表達不確定性，如蒙特卡羅(Monte Carlo)法[41]、拉丁超立方抽樣法(Latin hypercube sampling，LHS)[42]、普適似然不確定估計方法(generalized likelihood uncertainty estimation，GLUE)[43]以及單純多邊形進化算法(shuffled complex evolution algorithm，SCE-UA)[44]等。

表1 常用地表水環境數學模型

敏感性分析方法主要分為局部分析方法和全局分析方法。常采用的局部分析方法是檢驗單個參數的變化對模型結果的影響程度(one factor at a time，OAT)。局部分析方法簡單、易于實施，但響應結果較片面，無法解決“異參同效”問題。全局分析方法克服了局部分析方法的缺點，能夠反映整體參數組合對結果輸出不確定性的影響。全局敏感性分析方法有很多[45]：①基于回歸或相關分析技術的方法，如多元回歸法、響應曲面方法(response surface methodology，RSM)；②全局篩選法，如LH-OAT方法、Morris方法；③基于方差理論的方法，如傅里葉振幅敏感性檢驗法(Fourier amplitude sensitivity test，FAST)、Sobol方法和擴展傅里葉振幅敏感性檢驗法(extend FAST)；④因子設計實驗[46]；⑤摩爾斯分析法[47]；⑥取樣分析法[48]等。

敏感性分析常用的蒙特卡羅法屬于取樣分析法，能簡單有效地評價多個參數對模型輸出結果不確定性的貢獻。LHS克服了蒙特卡羅法計算成本高的缺點，抽取的樣本能更精確地反映輸入概率函數的分布，不僅高度控制抽樣值，又為它們留有變化的余地。基于LHS的思想，運行模擬的次數由輸入變量數決定，最少可為隨機變量數的1.5倍，一般為幾百次[49]。

當前國內外已有對模型敏感性和不確定性研究的成功案例。李一平等[50-51]利用EFDC模型和LHS抽樣方法分析了水動力模塊中的5個重要參數以及4個重要的外部輸入條件對太湖水位和流場分布的敏感性和不確定性，并借助原位觀測簡化了模型參數的率定驗證；Gong等[52]基于SWAT(soil and water assessment tool)模型采用GLUE不確定性分析方法分析模型參數對模型不確定性的影響，其中針對不同的研究目標需要選取不同的似然函數[53]；Reder等[54]采用LHS和全局敏感性分析法(global sensitivity analysis，GSA)在42個參數中篩選出4個最敏感的參數以提高模型性能。

除了傳統分析方法以外，基于熵的不確定性敏感性分析方法也被提出[55]，并用以識別不確定變量對隨機變量和模糊變量組合系統的影響[56]。全局靈敏度測算方法(sobol)是一種兩階段不確定性量化方法，簡化了不確定度量化計算過程[57]。風場、邊界及底部地形的空間變化導致研究區域內的參數分布也有所差異，因此結合GIS-Lab衍生的參數空間不確定分析將成為新的研究方向。

參數不確定性和敏感性分析已成為模型構建的必要工作，通過分析可表征輸出結果的不確定性對輸入參數不確定性的依賴程度。對于特定的模擬目標，篩選出敏感參數、參數合理分布范圍以及各個參數對模型結果不確定性的貢獻率，可大幅度減少模型后續率定驗證工作，同時可指導關鍵參數的監測工作。

4 地表水環境數學模型模擬精度關鍵影響因素

除由參數不確定性引起的模擬誤差外，影響模型模擬精度的因素還包括模型類型及選擇、網格種類、垂向坐標系統、方程離散方法、初始條件及邊界條件等。

4.1 模型類型及選擇

選擇合適的模型是建模前應明確的重要環節，需要參考模擬對象的特性、精度要求和計算效率等。當模擬對象具有干濕交替或洪泛區特性時，應選用較為精細的水動力模型[58]；若模擬深水湖庫水質，需要選用垂向二維或三維模型，考察垂向溫度和水質變化[59]；若模擬重金屬或有毒物質遷移轉換過程，關鍵是構建并校驗泥沙模塊，其余影響因素和模塊可進行一定程度的概化處理[60]。由于復雜的物理、化學、生物過程無法在模型中詳述，所以模型選擇不恰當可能會浪費計算時間，且不一定能達到預期目標。因此，模型的選擇并非越復雜越好，在明確模擬目標后，選擇合適的模型即可，這個環節需要一定的實踐經驗。

4.2 網格種類及垂向坐標系統

網格種類及坐標系統選擇也直接影響模擬精度。常見的平面網格有矩形網格、正交曲線網格、三角形網格及四叉樹網格。矩形網格基于直角坐標系便于組織數據結構，計算效率高，但不適合處理復雜的邊界且不易調控網格密度。正交曲線網格是一種基于曲線坐標系統的有結構網格，可以適應不規則邊界，但處理過于復雜的邊界時效果不佳。三角形網格利于研究復雜地形和邊界問題，易于控制網格密度，但計算效率較低。具有樹狀結構的四叉樹網格能高度擬合復雜的自然水體，易實現水動力及物質輸運數值模擬[61]。在邊界和地形較復雜的位置宜采用三角形網格，在計算域內部和地形變化不大的地方宜采用矩形網格或者正交曲線網格。

垂向坐標系統一般分為平面(z)坐標、等密度(ρ)坐標和地形擬合(σ)坐標，分別對應不同的網格。z坐標模式方程簡單，易于數值離散，適用于具有準水平運動特點的水體，但不易處理底部邊界，在淺水區域難以滿足必要的垂直分辨率。ρ坐標模式常用于密度流模擬，但在混合層、非層化水體內和底部邊界層的分辨率較低。σ坐標模式實現了垂直相對分層，能夠有效擬合底部地形，但斜壓梯度力會有較大截斷誤差。對于地形復雜的水域可采用混合坐標模式，例如σ-z坐標系統適用于局部陡峭深水區，能夠降低水平壓力梯度誤差的影響。在實際應用中，應根據岸線形狀、水下地形數據和模擬精度的要求選擇合適的網格及垂向坐標系統。

精細的分辨率能降低模擬誤差，但會增加不必要的模擬成本。若部分區域的模擬精度要求較高，可進行局部網格加密處理，既保證擬合計算的合理性，又提高運行效率。

4.3 方程離散方法

為確保計算的穩定性和收斂性，時間步長應該足夠小，通常須將時間步長減小到幾分或幾秒，與水動力過程模擬所需的時間步長差不多，可重現泥沙輸送等水質動力學過程。目前基于有限差分法、有限體積法和有限單元法的高效數值離散計算方法、并行計算方法、集群計算方法等同時滿足了精度和效度的要求[62]。有限差分法根據時間和空間步長對定解區域進行網格劃分，用差商代替導數，求解方法簡單但不易處理復雜邊界問題。有限體積法可以根據實際問題的物理特點對任意形狀網格體進行積分，且不會影響計算精度和守恒性，但是對于質量差的網格，離散的過程會產生更大的誤差。有限單元法根據實際問題的物理特點對求解區域進行單元剖分，能夠滿足一定的精度要求，但對于二維和三維問題需要建立許多人為的節點，且求解精度過于依賴網格劃分，適合分析連續變形問題。

4.4 初始條件及邊界條件

初始條件包括初始水位、流速、溫度等，初始條件設置是否準確，對不同的水體影響不同。模擬湖庫水質時，若只依賴邊界條件驅動，運行至合理的初始狀態需要一定的時間且對模擬的結果影響較大。因此，模擬湖庫水質時，如沒有準確的實測值作為初始條件，則通常將模型運行一段時間后的結果作為初始態，稱為“預熱”。對于河流、河口等水動力過程較劇烈的水體，邊界條件的驅動將很快覆蓋初始條件，初始條件的影響很小。

邊界條件包括大氣邊界、出入流邊界、開邊界的作用力、水工建筑物、取水退水邊界，其中水質模塊的邊界條件還需考慮各種出入水體邊界的水質變量、大氣干濕沉降、農田面源污染、地表徑流、內源釋放、地下水等。在確定污染源的過程中，大氣干濕沉降對水域整體的貢獻不可忽略，隨降雨進入水體的污染物質往往是水質惡化的關鍵因素[63]。邊界條件還可分為垂向和水平邊界條件，比如氣溫和風速作為模型垂向邊界條件，雖然不參與直接計算模擬，但它們影響了潮流、混合和熱傳輸等水動力過程。合適的邊界條件可以有效避免模擬誤差。

5 地表水環境數學模型發展趨勢

目前，地表水環境數學模型正朝著系統化、綜合化、法規化方向發展，模型涉及的要素越來越復雜，與許多新興技術的結合，使得地表水環境數學模型的發展充滿了機遇與挑戰。

a. 模型系統的系統化、綜合化和平臺化。水環境模型逐步涉及水文、水動力、生物、地理等多領域，以單獨模塊的形式集成一體。模擬元素和模擬情景的增加拓展了模型的研究領域，比如藥品及個人護理用品(PPCPs)、微塑料等新型污染物將成為新的模擬對象，藻類競爭、熱量平衡等過程也逐步引入水生態模型。適用于大型流域的三維水生態模型開發，以及利用統一平臺構建促使數值模擬成為流域控制與規劃決策的支撐技術，將成為未來重要的發展趨勢。這些將會是今后地表水環境數學模型研究的熱點。

b. 與新興技術的結合。大數據平臺的參與能夠有效解決邊界條件、初始條件輸入缺失難題；超算、云計算以及人工智能算法的應用將會大幅度提升地表水環境數學模型的計算效率和精度，比如遺傳算法、模擬退火算法能夠強化參數識別；VR技術豐富了模型結果的輸出展現方式，實現了人性化的設計理念。在與物聯網、“互聯網+”、云技術的碰撞下，“3S”技術與水環境模擬集成模塊的開發促進了環境因素之間的相互作用和動態變化的確定性分析，基于遙感技術的水質反演模型促進了“天地一體化”水環境監測系統的建立，能夠實現數字預警、識別黑臭水體、考察海域污染，體現實時性、大尺度、高速度、動態性等優勢。在GIS-Lab技術的引入結合下，風場、地形、床層等空間因素的考慮會更加全面，不確定性分析正從全局概化向區域性精細化發展，模擬過程更加真實，結果更加可靠。

c. 法規化趨勢?；贕IS技術，模型庫管理系統的建立促使水環境模型走上法規化管理道路。目前，美國國家環境保護局(EPA)已將97種地表水環境數學模型列入模型信息庫，澳大利亞政府也針對模型選擇、參數率定、敏感性分析等給出了系統推薦，我國生態環境部也正式發布了HJ2.3—2018《環境影響評價技術導則 地表水環境》，并推薦了適用于河流、湖泊、河口和海洋的數值模型。今后的發展將基于已有的模型信息庫進行拓展豐富，有望實現模型運用流程的規范化操作、指導、監督。