土壤中硝態氮遷移轉化的數值模擬研究

(中國地質大學(武漢) 環境學院，湖北 武漢 430074)

近年來農業生產中氮肥的施用量大大增加，氮素在土壤中的存在形態可以分為有機氮、硝態氮、銨態氮和亞硝態氮[1]。由于硝態氮不易被土壤膠體吸附[2]，易隨降雨和灌溉水從土壤表層向土壤深層遷移，進而淋失到地下水中，造成地下水污染。所以，深入研究硝態氮在土壤中的淋失過程，對合理施用氮肥和控制地下水污染均有重要意義。

在影響硝態氮淋失的眾多因素中，灌溉量和施氮量是極其重要的兩個因素[3]。研究土壤中硝態氮的方法主要為田間試驗和土柱試驗。魯艷紅[4]以兩系雜交稻為研究對象，研究了不同施氮量對水稻產量形成及產量構成因子的影響，分析了不同施氮水平下氮素吸收利用效率的差異。馮紹元[5]通過在北京順義區進行模擬降雨的田間試驗，研究了不同降雨與施肥水平對玉米土壤硝態氮分布與累積的影響。張海濤[6]通過開展室內土柱實驗研究了碳源對潛流帶中氮素遷移轉化的影響。隨著計算機技術的發展，現階段對硝態氮的研究已優化為田間試驗、室內試驗與數值模擬相結合的方法。何佳吉[7]通過對宜興市梅林流域土壤中水分運動以及硝態氮淋濾過程進行動態模擬分析，揭示了氮肥施用與硝態氮淋失之間的關系。Arash Tafteh[8]通過數值模擬，分析評價了施用不同氮肥對土壤中硝酸鹽的浸出率和變化量的影響。Li[9]則運用試驗與數值模擬相結合的方法研究了農田地表徑流中硝態氮的淋失規律。

洞庭湖地區是我國重要的糧油生產和養殖業生產基地，常德市位于洞庭湖西部，獨特的氣候和豐富的水土資源極利于農作物的生長，該地區糧食、棉花總產值居湖南之首。然而由于該地區長期施用大量氮肥，有些未被植物吸收利用的氮肥經土壤入滲到含水層中造成地下水污染。本文以常德市典型淺層土壤為研究對象，通過室內靜態試驗和動態土柱淋濾試驗獲得硝態氮靜態條件和淋濾條件下的遷移轉化規律及相關參數。在此基礎上，利用數值模擬軟件建立數值模型，并通過校正好的模型模擬不同灌溉強度下硝態氮濃度的變化特征，探討分析灌溉強度對硝態氮濃度的影響，從而為該地區氮肥施用和地下水保護提供科學依據。

1 物理試驗模型

1.1 土柱試驗裝置及試驗方法步驟

1.1.1試驗裝置

土柱淋濾裝置由用于穩定水頭的馬氏瓶、土柱和取樣瓶組成。土柱主體為有機玻璃柱，內徑 8 cm，高度50 cm，側壁36 cm和48 cm 深度處設有取樣口。實驗裝置簡圖如圖1所示。

1.1.2土樣性質及裝填

本次試驗土樣采集于常德市辰陽村 0～20 cm 耕作層，土樣采集后經自然風干過2 mm篩，按《土壤農化分析》測定土樣基本理化性質。經測定，土樣基本理化性質如下所示。

項目粉質黏土容重1.41g/cm3Eh136mVpH6.54NH+4-N1.98mg/kgNO-3-N5.92mg/kgNO-2-N0.36mg/kg

裝填土樣時先在柱子底部裝入2 cm 厚的石英砂防止土壤沖散堵塞出水口；之后將風干的土樣按實測容重換算稱重后裝入柱內，隔1 cm分層裝填，使土樣盡量壓實均勻；最后再在土柱的頂部裝入2 cm厚的石英砂，防止淋濾液進入土柱時沖散土壤，使之均勻進水。

1.1.3淋濾液的配置

根據湖南省統計年鑒，2015年湖南省年降雨量為1 580.50 mm，占30%循環地表徑流；耕地面積415.35萬hm2，全省氮肥使用量為2 016 259 t。為便于計算，配置40 mg/L 的硝酸鈉溶液(以N計)進行淋濾。

1.1.4淋濾試驗及取樣分析

首先從土柱底部注入去離子水，排除土柱內空氣。當土柱頂部有水溢出時，改從土柱頂部進水。然后通過控制裝置使土柱中流場保持穩定，改換配置好的淋濾液進行淋濾。

由于土樣滲透系數較小，每隔3 d從取樣口取樣，測定銨態氮、亞硝態氮、硝態氮含量。

采用納氏試劑分光光度法檢測銨態氮，采用α-萘胺光度法檢測亞硝態氮，采用紫外分光光度法檢測硝態氮。

1.2 土柱試驗結果

圖2 硝態氮濃度隨時間變化曲線Fig.2 The temporal variation of concentrations at different soil depths

2 數學模型

2.1 水流運動基本方程

在不考慮土壤水平和側向水流運動時，即在一維垂向運移中，水流運動可以用 Richards 方程來描述[10-11]，方程如下：

(1)

式中，θ為土壤體積含水率；h為土壤壓力水頭；t為模擬時間；α為流向與垂直向夾角，根據上述物理試驗模型，其水流為一維垂向滲流，即α=0；K為非飽和滲透系數；S為源匯項。

2.2 水分運移方程

土壤水分運移模型可用來描述水分在土壤中的運移過程。HYDRUS-1D 軟件水流模型中包括單孔介質模型、雙孔隙/雙滲透介質模型等多種土壤水分運移模型。本文模擬時采用 Van Genuchten-Mualem 提出的土壤水力模型來進行模擬預測，且在模擬中不考慮水流滯后的現象[12]，方程為

(2)

K(h)=KsSel[1-(1-Se1/m)n]2

(3)

(4)

(5)

式中，θr為土壤殘余含水率；θs為土壤飽和含水率；Se為有效飽和度；α為冒泡壓力；n為土壤孔隙大小分配指數；Ks為飽和水力傳導系數；l為土壤孔隙連通性參數，通常取0.5。

2.2 溶質運移方程

土壤溶質運移主要有3個過程，分別為對流、分子擴散和機械彌散。本文模擬時采用經典的對流-彌散方程描述飽和-非飽和孔隙介質中的一維溶質運移[13]。

(6)

式中，c為溶液液相濃度;ρ為土壤容重;s為溶質固相濃度；D為綜合彌散系數；q為體積流動通量密度；S為源匯項。

3 數值模擬

3.1 網格剖分和時間步長

根據土柱淋濾試驗，模擬時將模擬土柱均等化分，進行1 cm等距剖分，并在36 cm和48 cm設置觀測點。初始時間步長設為0.001 d，最小時間步長設為0.001 d，最大時間步長設為0.1 d，迭代控制參數則使用默認值。

3.2 初始條件與邊界條件

由于氮的轉化過程十分復雜，本文建立模型時只考慮銨態氮、亞硝態氮和硝態氮之間的轉化。進行土柱淋濾試驗時采用馬氏瓶從土柱上端連續進水保持進水壓力水頭不變，因此水流模型上邊界選用恒定壓力水頭邊界(Constant Pressure Head)，水流模型的下邊界概化為包氣帶底部(即潛水面)，所以下邊界選用自由排水邊界(Free Drainage)；土柱試驗采用40 mg/L 的硝酸鈉持續淋濾，所以溶質運移模型的上邊界設定為恒定濃度邊界(Constant Concentration Boundary Condition)，下邊界設定為零濃度梯度邊界(Zero Concentration Gradient Boundary)，代表初始狀態為液相0濃度狀態。設置初始水頭時，軟件根據介質剖分的節點數對初始水頭進行離散。將模型頂部的壓力水頭設為0 cm，表示初始時刻土壤表層處于近飽和狀態，而土柱底部的壓力水頭設為-48 cm，初始壓力水頭自上而下為0～-48 cm 均勻分布。

3.3 參數確定

數值模擬溶質運移的主要參數包含土壤水力特征參數、溶質運移特征參數和溶質反應特征參數[14]。

3.3.1土壤水力特征參數確定

將過2 mm篩的風干土樣采用馬爾文激光儀進行粒度分析試驗，試驗結果顯示該土樣黏粒(0.01～2.00 μm)所占百分比為3.60%，粉粒(2.00～50.00 μm)所占百分比為96.4%。利用模擬軟件水流運動參數界面中的神經網絡預測功能(Neural Network Prediction )，根據土壤質地和容重來進行監測 VG 型水分曲線的水力學參數的預測。輸入土壤砂粒、粉粒及黏粒的百分比組成及土壤容重，即可以得出曲線中所需要的參數。輸入數據并根據模擬結果調整參數得出土樣各水力特征參數如表2所示。

表1 土壤水力特征參數Tab.1 The hydraulic characteristic parameters of the soil

3.3.2反硝化反應速率測定

表2 反硝化反應轉化量及轉化率Tab.2 Denitrification reaction conversion and conversion rate

由實驗數據可知，0～1 d時的反硝化作用較為微弱，反硝化作用對于硝酸氮濃度的影響主要發生在1～4 d，4 d后反硝化反應速率緩慢，所以以 1～4 d的硝酸氮濃度變化曲線計算反硝化動力學方程。經過擬合發現，1～4 d反硝化反應符合一級反應速率方程，速率常數為 0.399，硝態氮濃度呈指數形式下降。

圖3 硝態氮濃度隨時間變化曲線Fig.3 The temporal changes of concentrations at different time

表3 氮素遷移轉化參數Tab.3 N transport and transformation parameters for the soil profile

3.4 模型校正及檢驗

輸入土壤性質參數驗證模型的可靠性，以土柱不同深度觀測點的硝態氮含量實測值及其淋濾液濃度為模型的初始邊界條件，模擬計算硝態氮的變化情況。將實測硝態氮的濃度與模擬值進行對比，圖4 是硝態氮模擬值與實測值隨時間變化的比較結果。由圖可知：對土柱 2 個深度處的觀測點進行取樣監測，測得土壤溶液中硝態氮濃度模擬曲線與實測曲線變化形態基本一致，模擬值與實測值相差不大。1～4 d模擬 48 cm 處的硝態氮濃度擬合要優于 36 cm 處，這主要是因為硝化作用主要發生于 1～4 d，由于土柱底部與外界空氣相連通，48 cm 處較 36 cm 溶解氧多，抑制了反硝化作用。而模擬時則認為模型各點反硝化反應均為一級反應且反應常數相同，輸入模型的靜態反硝化反應速率常數測定環境更近似于 48 cm 土柱處的環境，所以 1～4 d土柱 36 cm 處硝態氮濃度模擬值高于實測值。而在試驗第 10 天 48cm 處實測土壤硝態氮含量有明顯升高，與模擬值偏離較大，這可能是因為取樣點變化造成的土壤背景值差異所導致。通過 SPSS軟件進行模擬值與實測值的距離相關分析(以 Pearson 相關系數為距離)，36 cm 和48 cm 硝態氮含量實測值與模擬值相關系數分別為 0.998 和 0.995 ，皆達到顯著水平，說明模擬效果較為理想，模型可以較好反映硝態氮在土壤水中的運移情況。

圖4 實測硝態氮的濃度與模擬值對比曲線Fig. 4 Comparison of the measured and

4 不同灌溉強度硝態氮濃度的變化特征

常德市雖然雨量充沛，但降雨較集中，缺水季節需要對農作物進行灌溉。利用校正后的模型對不同灌溉強度下硝態氮濃度變化進行模擬。該地區土壤滲透系數約為 4 cm/d，將灌溉強度分別設為1，2，4，8，12 cm/d，運行模型。模擬結果如圖 5 所示。

圖5 不同灌溉強度硝態氮濃度Fig.5 The values of concentrations in different irrigation intensity scenarios

各組灌溉強度下硝態氮的濃度最大值相差不大，硝態氮濃度總體變化為先升高后降低。不同灌溉強度影響入滲到土壤中的水量，引起土壤剖面含水率變化，進而影響土壤硝態氮濃度變化。當灌溉強度大于 4 cm/d 時，超過土壤的下滲能力，灌溉水按下滲能力下滲，多余的水會形成地面積水，進而形成地表徑流，所以模擬 10 d后各組土柱底部硝態氮的濃度差異不大；當灌溉強度小于等于4 cm/d 時，灌溉水全部下滲到土壤中，而硝態氮易溶于水中，所以模擬10 d后土壤中硝態氮濃度隨灌溉強度的增加而增大。

5 結 論

土壤中硝態氮的靜態反硝化作用主要發生于實驗開始后 1～4 d，在這期間硝態氮濃度呈指數形式下降，反硝化反應近似一級反應，反應速率常數為0.399/d，試驗4 d后反硝化速率緩慢，反硝化作用對硝態氮濃度影響比較微弱。由于硝態氮易溶于水中，易隨灌溉水下滲到土壤中，而灌溉水的入滲量受土壤下滲能力制約：當灌溉強度大于土壤下滲能力時，灌溉水按下滲能力下滲，多余的水會形成地面積水，進而形成地表徑流，硝態氮濃度的變化趨勢基本一致；當灌溉強度小于等于土壤下滲能力時，灌溉水全部下滲到土壤中，硝態氮濃度隨灌溉強度的增加而增大。通過控制灌溉強度，使其約為 4 cm/d，既可以提高水的利用，也可以減少硝態氮向土壤深層大量流失，控制硝態氮污染土壤和地下水。