河流潛流帶氮素遷移轉(zhuǎn)化數(shù)值模擬研究進(jìn)展

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

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

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

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

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

1 潛流帶氮素遷移轉(zhuǎn)化模擬常用數(shù)學(xué)模型

1.1 潛流帶水動力模型

1.2 潛流帶溶質(zhì)運移模型

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

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

1.3 潛流帶微生物生長模型

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

2 耦合模型主要參數(shù)及其影響因素

2.1 水動力參數(shù)

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

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

2.2 氨氮吸附-解吸行為參數(shù)

2.3 氮素生物地球化學(xué)反應(yīng)參數(shù)

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

3 數(shù)值模擬的不確定性問題

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

3.1 數(shù)據(jù)信息獲取

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

3.2 耦合模型構(gòu)建

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

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

3.3 模型參數(shù)識別

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

4 結(jié)論與展望

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

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