河流潛流帶氮素遷移轉化數值模擬研究進展

蔡 奕，邢婧文，阮西科，周念清，黃若堯，乙東澤

(1.同濟大學土木工程學院，上海 200092； 2.長江水環境教育部重點實驗室，上海 200092)

隨著工農業生產的快速發展，大量含氮生活污水、工業廢水以及農業面源污染等通過多種途徑進入河流，破壞了原有的營養鹽平衡，導致水質下降，從而威脅河流生態系統健康；而且河流中超負荷承載的氮素會隨水流向下游輸移，進入河口后會誘發近海水體富營養化，嚴重時會導致大量魚類和貝類等生物的死亡[1-2]。如何有效降低或去除河流中超標的氮素以恢復生態平衡，已成為當前水環境保護亟須解決的問題[3]。

河流潛流帶(以下簡稱“潛流帶”)是河床以下并延伸至兩側岸灘區的水分飽和沉積物含水層，地表水和地下水在此混合發生轉化，其內部環境條件復雜多變[4]。潛流帶中水體交換引發了碳、氮、磷、氧等物質交換和能量傳遞，即潛流交換，為潛流帶生物群落提供了適宜生存的條件。潛流帶生物地球化學反應強烈，對調節河流的生態具有重要作用[5]。潛流帶沉積物中碳源物質、氧化還原環境、生物群落等因素對氮素賦存形態和循環模式產生重要影響[6]。因此，要解決河流中氮素污染問題，必須對潛流帶氮素遷移轉化機理和運移規律進行深入研究。

研究潛流帶中氮素遷移轉化機理和規律的方法主要有現場監測、室內試驗和數值模擬等方法[7-9]。相比于現場監測和室內試驗，數值模擬不受點位布設和尺寸效應的影響，可定量描述氮素在潛流帶地下水中遷移轉化的連續行為，還可對其進行反演和預測，因而得到廣泛應用。由于潛流帶氮素遷移轉化過程復雜，影響因素多，精準量化存在一定的困難[10-11]。本文系統總結了潛流帶氮素遷移轉化模擬的研究現狀與存在的不足，并提出今后的發展方向，可為河流污染防治和水資源保護提供參考。

1 潛流帶氮素遷移轉化模擬常用數學模型

1.1 潛流帶水動力模型

1.2 潛流帶溶質運移模型

潛流帶溶質運移模擬常用對流-彌散模型。該模型基于連續性方程和質量守恒定律建立，包含對流、彌散及其他源(匯)項，通過水流速度與水動力方程耦合求解[19]。常見的源(匯)項有吸附、解吸等物理作用項及自然降解、微生物催化等生物地球化學反應項[20-21]。

①有氧呼吸；②硝化作用；③反硝化作用；④微生物同化作用；⑤植物同化作用；⑥氨化作用；⑦厭氧氨氧化作用圖1 潛流帶中氮素循環示意圖Fig.1 Schematic diagram of nitrogen cycle in hyporheic zone

1.3 潛流帶微生物生長模型

Hampton等[31-32]研究表明，好氧菌、硝化菌、反硝化菌等微生物的生長與死亡會影響潛流帶中氮素遷移轉化。潛流帶中微生物生長是一個自限性的過程：隨著微生物量(生物膜)的增長，沉積物孔隙空間逐漸被填充，生物堵塞導致底物擴散通量和沉積物滲透系數降低，使得物質和能量供給下降，造成微生物的死亡，而隨著生物堵塞程度的降低，營養物質得以補充，微生物又開始新一輪的生長繁殖[33-34]。Ping等[35-37]建立了微生物生物量與沉積物孔隙度(或滲透系數)之間的定量關系，涉及微生物生長速率、死亡速率、細胞密度、初始孔隙度等參數。將計算獲得的微生物量及孔隙度(或滲透系數)代入溶質運移模型和水動力模型，可實現潛流帶滲流場、化學場和生物場的信息互饋。

2 耦合模型主要參數及其影響因素

2.1 水動力參數

滲透系數和水動力彌散系數是反映氮素在地下水中對流和彌散運動的重要參數。滲透系數的大小主要受沉積物組成、沉積條件及飽和度等因素的影響。潛流帶沉積物是在流水中以機械方式沉積的碎屑物，通常具有沉積韻律性和流水成因的沉積構造特征，分選性較好，成層性較清晰，非均質性特征顯著。潛流帶既是河水和地下水的交換區，又是大多數無脊椎動物和微生物的生境，所以潛流帶沉積物孔隙結構會因水流沖淤、生物行為、顆粒物質溶解與沉淀等作用影響而呈現動態變化[38-39]。此外，河岸潛流帶是變飽和區，滲透系數與飽和度和含水量之間存在一定的關系。潛流帶的非均質性和動態性導致了沉積物滲透系數的時空差異性，其差異可達到幾個數量級[40-41]。由于潛流帶水-土-生物環境復雜且多變，滲透系數變化隨機性大，目前難以建立定量關系準確表征滲透系數的空間分布特征和動態變化特征。因此，用于模擬計算的滲透系數主要是通過室內試驗、野外現場測定或經驗估算確定的。作為水文地質的重要參數，滲透系數的取值會直接關系到滲流速度的計算并影響氮素遷移轉化模擬結果[42]?？梢?，滲透系數的準確測定是潛流帶氮素遷移轉化模擬的重要基礎。

水動力彌散系數包括機械彌散系數和分子擴散系數，通常采用室內或現場的彌散試驗確定。Godoy等[43]研究表明，室內試驗存在尺寸效應，其測定值往往比野外測量值小，兩者之間可能有數量級上的差異。大量的研究成果顯示影響水動力彌散系數的因素有很多，如孔隙水流速、介質粒徑大小、溶質性質、溫度、遷移距離等。在低流速條件下，分子擴散作用不可忽略，水動力彌散系數與流速之間并非正比關系[44]?？紫督橘|粒徑的增大會導致水動力彌散系數的增大[45]。飽和土的水動力彌散系數與孔隙水流速存在顯著的相關性，非飽和土的水動力彌散系數與含水量有關且與濃度無關，無論是飽和土還是非飽和土，吸附性溶質的水動力彌散系數通常要大于非吸附性溶質[46]。此外，溫度的升高和遷移距離的增大也會導致水動力彌散系數的增大[47]。受河水、地下水和大氣溫度差異的影響，潛流帶的溫度場時刻處于變化狀態?？梢?，潛流帶的非均質性和動態性會導致水動力彌散系數的時空差異性。在已有的相關數值模擬研究中，水動力彌散系數通常采用定值，其時空和溶質差異性并沒有得到很好的體現，這在一定程度上也影響了氮素遷移轉化模擬的精度。

2.2 氨氮吸附-解吸行為參數

2.3 氮素生物地球化學反應參數

潛流帶中氮素生物地球化學作用受多環境因子控制，包括氧化還原電位、pH值、溫度、濃度梯度、水流條件、微生物群落等[53]。Pescimoro等[54-56]研究表明，DO和DOC是潛流帶氮素生化反應類型和進程的重要控制因素。在沉積物環境下，反硝化細菌酶在25～35℃范圍內活性高，反硝化反應速率快，氮素去除效率高[57]。當氧化還原電位范圍處于100～350 mV時，反硝化速率隨著氧化還原電位增大而減小[58]。異養反硝化菌適宜的pH值為5.5～8.0，此時的反硝化速率相對較高，pH值小于5.0的酸性環境會抑制反硝化作用[59]。以硝化作用為主的生化反應會隨著DO的消耗而轉變為以反硝化作用為主，常通過建立DO消耗與河水滯留時間之間的函數關系實現對這兩類反應的動態識別[60-61]?？梢?，與水動力參數和吸附-解吸參數相似，潛流帶氮素生化反應參數并非為常數，而是會隨著滲流運動和環境變化而發生改變。不過，在目前的模擬研究中并未考慮生化反應參數與環境因子和水流條件之間的關系。

3 數值模擬的不確定性問題

因受觀測數據和控制方程的限制，潛流帶氮素遷移轉化模擬結果往往與真實值存在一定的偏差。建模時應充分考慮信息數據、數學模型及模型參數對模擬結果的影響，以降低過程模擬的不確定性。

3.1 數據信息獲取

潛流帶是地表以下的黑暗生境，其滲流場、化學場、氧化還原環境、生物菌類活性等信息獲取很容易受人為干擾而失真，所以潛流帶原生環境信息的準確獲取是高精度數值模擬的前提。目前，針對潛流帶的研究主要借助于觀測井，通過在觀測井中放置水位、溫度、電導率、pH值、氧化還原電位等傳感器或采集水、土樣進行分析，以獲取數值模擬所需數據[62-63]。然而，潛流帶孔隙水水質、細菌分布特征等相關數據會隨采樣方式或條件的不同而發生變化。例如，觀測井與大氣連通會影響水中DO含量，可能會導致原先的還原環境轉變為氧化環境，使得氮素存在形態發生改變；氧氣、光照、溫度等采樣條件控制不當會造成沉積物土樣中微生物生物量測定失真；采集頻率低可能會導致多場信息變化幅度難以捕捉。因此，需要結合高分辨率高靈敏度的原位傳感監測、非擾動測試及情景模擬試驗等技術獲取原生環境信息數據，這對氮素遷移轉化模擬研究具有重要意義。

3.2 耦合模型構建

潛流帶氮素遷移轉化行為十分復雜，是滲流場、溫度場、化學場、生物場等多場相互作用的結果，故其數值模擬的關鍵在于描述多場時空演變過程及其相互之間的耦合聯動關系。目前的數值模擬主要關注水動力、溶質運移、微生物生長等過程的耦合計算，其中溶質運移計算通?？紤]對流項、彌散項、反應項。由于土顆粒對氨氮具有一定的吸附作用，氨氮運移模型應增加吸附項，否則獲得的氨氮時空分布規律與真實情況不符。朱靜思等[64-65]研究表明，溫度變化會影響氮素的水動力彌散系數、吸附速率和平衡常數及微生物催化反應速率等，進而控制氮素的循環過程。潛流帶溫度場常伴隨著水量交換發生顯著的動態變化[66-67]，但在現有的耦合模型中卻極少考慮溫度的時空演變。為了對潛流帶氮素遷移轉化有更全面的認識，有必要建立模型參數與溫度、pH值等環境因子之間的定量關系，通過水動力、熱傳導、溶質運移及微生物生長模型的耦合計算，實現多過程之間的動態聯動，如圖2所示。此外，目前大多數的耦合模型采用一、二維數學模型[68-69]，然而，潛流帶具有復雜的三維空間非均質各向異性結構，采用三維模型更符合潛流帶氮素遷移轉化過程的定量化表達。

圖2 潛流帶氮素遷移轉化模型中不同過程耦合示意圖Fig.2 Schematic diagram of coupling different processes involved in the model of nitrogen migration and transformation in hyporheic zones of rivers

3.3 模型參數識別

潛流帶氮素遷移轉化過程耦合模型結構復雜，涉及參數多，如滲透系數、水動力彌散系數、反應物最大生化反應速率等，其中一些參數在維度、位置或時間上有差異，參數不確定問題是客觀存在的?；诮涷灩浪慊蛴^測數據優化所獲取的參數通常不能保證模型在各種可能情景下的模擬精度。而且，因不同算法在參數收斂軌跡上存在差異，優化結果往往也不盡相同，甚至無法判斷其是否達到全局最優[70]。鑒于優化方法在復雜模型的參數識別上有一定的局限性，參數敏感性研究越來越受到關注。目前，潛流帶氮素遷移轉化模型的參數敏感性研究主要采用基于“最佳”估計值的擾動分析方法，其研究結果表明敏感參數主要為滲透系數、水動力彌散系數、最大生化反應速率及半飽和常數等[71-73]。此外，多參數識別通常要比單參數要復雜得多，這是因為復雜模型結構中不同參數之間可能會存在某種相關性，單獨考慮各參數的“最佳”估計值未必達到最優的模型效率[74]?？梢?，潛流帶氮素遷移轉化模型的參數識別需要考慮參數時空分布的復雜性、靈敏度及各參數之間的相關性，而且合理的算法應充分體現參數不確定性而不是參數最優化。然而，這些參數識別中所要關注的問題現有研究成果還較少，還有待進一步研究。

4 結論與展望

b.深入研究潛流帶氮素遷移轉化中多過程間的耦合聯動關系。通過試驗研究探討潛流帶滲流場、溫度場、化學場、生物場的時空演變規律，厘清潛流帶水動力、熱傳導、溶質遷移及微生物活動等多過程間的耦合聯動關系，建立主要模型參數與溫度、pH值等環境因子間的定量關系，在此基礎上對現有耦合模型進一步完善，如耦合熱傳導過程，考慮氨氮吸附-解吸行為和其他氮素生化反應等。

c.開展多維度多過程耦合模型優化計算研究。深入研究潛流帶氮素遷移轉化模型參數的復雜性、敏感性及相關性，提出基于不確定性分析的參數識別算法。分析不同過程的持續時間和涉及范圍，合理設置各過程的時間步長和計算網格，以實現多過程間的同步與異步模擬。此外，三維的多過程耦合會消耗大量計算資源，需加強模擬軟件的開發與集成，優化空間網格，強化并行計算，節省運行時間，以提升模擬精度和效率。