(中國(guó)地質(zhì)大學(xué)(武漢) 環(huán)境學(xué)院，湖北 武漢 430074)

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

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

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

1 物理試驗(yàn)?zāi)Ｐ?/h2>

1.1 土柱試驗(yàn)裝置及試驗(yàn)方法步驟

1.1.1試驗(yàn)裝置

土柱淋濾裝置由用于穩(wěn)定水頭的馬氏瓶、土柱和取樣瓶組成。土柱主體為有機(jī)玻璃柱，內(nèi)徑 8 cm，高度50 cm，側(cè)壁36 cm和48 cm 深度處設(shè)有取樣口。實(shí)驗(yàn)裝置簡(jiǎn)圖如圖1所示。

1.1.2土樣性質(zhì)及裝填

本次試驗(yàn)土樣采集于常德市辰陽(yáng)村 0～20 cm 耕作層，土樣采集后經(jīng)自然風(fēng)干過(guò)2 mm篩，按《土壤農(nóng)化分析》測(cè)定土樣基本理化性質(zhì)。經(jīng)測(cè)定，土樣基本理化性質(zhì)如下所示。

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

裝填土樣時(shí)先在柱子底部裝入2 cm 厚的石英砂防止土壤沖散堵塞出水口；之后將風(fēng)干的土樣按實(shí)測(cè)容重?fù)Q算稱(chēng)重后裝入柱內(nèi)，隔1 cm分層裝填，使土樣盡量壓實(shí)均勻；最后再在土柱的頂部裝入2 cm厚的石英砂，防止淋濾液進(jìn)入土柱時(shí)沖散土壤，使之均勻進(jìn)水。

1.1.3淋濾液的配置

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

1.1.4淋濾試驗(yàn)及取樣分析

首先從土柱底部注入去離子水，排除土柱內(nèi)空氣。當(dāng)土柱頂部有水溢出時(shí)，改從土柱頂部進(jìn)水。然后通過(guò)控制裝置使土柱中流場(chǎng)保持穩(wěn)定，改換配置好的淋濾液進(jìn)行淋濾。

由于土樣滲透系數(shù)較小，每隔3 d從取樣口取樣，測(cè)定銨態(tài)氮、亞硝態(tài)氮、硝態(tài)氮含量。

采用納氏試劑分光光度法檢測(cè)銨態(tài)氮，采用α-萘胺光度法檢測(cè)亞硝態(tài)氮，采用紫外分光光度法檢測(cè)硝態(tài)氮。

1.2 土柱試驗(yàn)結(jié)果

圖2 硝態(tài)氮濃度隨時(shí)間變化曲線Fig.2 The temporal variation of concentrations at different soil depths

2 數(shù)學(xué)模型

2.1 水流運(yùn)動(dòng)基本方程

在不考慮土壤水平和側(cè)向水流運(yùn)動(dòng)時(shí)，即在一維垂向運(yùn)移中，水流運(yùn)動(dòng)可以用 Richards 方程來(lái)描述[10-11]，方程如下：

(1)

式中，θ為土壤體積含水率；h為土壤壓力水頭；t為模擬時(shí)間；α為流向與垂直向夾角，根據(jù)上述物理試驗(yàn)?zāi)Ｐ?，其水流為一維垂向滲流，即α=0；K為非飽和滲透系數(shù)；S為源匯項(xiàng)。

2.2 水分運(yùn)移方程

土壤水分運(yùn)移模型可用來(lái)描述水分在土壤中的運(yùn)移過(guò)程。HYDRUS-1D 軟件水流模型中包括單孔介質(zhì)模型、雙孔隙/雙滲透介質(zhì)模型等多種土壤水分運(yùn)移模型。本文模擬時(shí)采用 Van Genuchten-Mualem 提出的土壤水力模型來(lái)進(jìn)行模擬預(yù)測(cè)，且在模擬中不考慮水流滯后的現(xiàn)象[12]，方程為

(2)

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

(3)

(4)

(5)

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

2.2 溶質(zhì)運(yùn)移方程

土壤溶質(zhì)運(yùn)移主要有3個(gè)過(guò)程，分別為對(duì)流、分子擴(kuò)散和機(jī)械彌散。本文模擬時(shí)采用經(jīng)典的對(duì)流-彌散方程描述飽和-非飽和孔隙介質(zhì)中的一維溶質(zhì)運(yùn)移[13]。

(6)

式中，c為溶液液相濃度;ρ為土壤容重;s為溶質(zhì)固相濃度；D為綜合彌散系數(shù)；q為體積流動(dòng)通量密度；S為源匯項(xiàng)。

3 數(shù)值模擬

3.1 網(wǎng)格剖分和時(shí)間步長(zhǎng)

根據(jù)土柱淋濾試驗(yàn)，模擬時(shí)將模擬土柱均等化分，進(jìn)行1 cm等距剖分，并在36 cm和48 cm設(shè)置觀測(cè)點(diǎn)。初始時(shí)間步長(zhǎng)設(shè)為0.001 d，最小時(shí)間步長(zhǎng)設(shè)為0.001 d，最大時(shí)間步長(zhǎng)設(shè)為0.1 d，迭代控制參數(shù)則使用默認(rèn)值。

3.2 初始條件與邊界條件

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

3.3 參數(shù)確定

數(shù)值模擬溶質(zhì)運(yùn)移的主要參數(shù)包含土壤水力特征參數(shù)、溶質(zhì)運(yùn)移特征參數(shù)和溶質(zhì)反應(yīng)特征參數(shù)[14]。

3.3.1土壤水力特征參數(shù)確定

將過(guò)2 mm篩的風(fēng)干土樣采用馬爾文激光儀進(jìn)行粒度分析試驗(yàn)，試驗(yàn)結(jié)果顯示該土樣黏粒(0.01～2.00 μm)所占百分比為3.60%，粉粒(2.00～50.00 μm)所占百分比為96.4%。利用模擬軟件水流運(yùn)動(dòng)參數(shù)界面中的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)功能(Neural Network Prediction )，根據(jù)土壤質(zhì)地和容重來(lái)進(jìn)行監(jiān)測(cè) VG 型水分曲線的水力學(xué)參數(shù)的預(yù)測(cè)。輸入土壤砂粒、粉粒及黏粒的百分比組成及土壤容重，即可以得出曲線中所需要的參數(shù)。輸入數(shù)據(jù)并根據(jù)模擬結(jié)果調(diào)整參數(shù)得出土樣各水力特征參數(shù)如表2所示。

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

3.3.2反硝化反應(yīng)速率測(cè)定

表2 反硝化反應(yīng)轉(zhuǎn)化量及轉(zhuǎn)化率Tab.2 Denitrification reaction conversion and conversion rate

由實(shí)驗(yàn)數(shù)據(jù)可知，0～1 d時(shí)的反硝化作用較為微弱，反硝化作用對(duì)于硝酸氮濃度的影響主要發(fā)生在1～4 d，4 d后反硝化反應(yīng)速率緩慢，所以以 1～4 d的硝酸氮濃度變化曲線計(jì)算反硝化動(dòng)力學(xué)方程。經(jīng)過(guò)擬合發(fā)現(xiàn)，1～4 d反硝化反應(yīng)符合一級(jí)反應(yīng)速率方程，速率常數(shù)為 0.399，硝態(tài)氮濃度呈指數(shù)形式下降。

圖3 硝態(tài)氮濃度隨時(shí)間變化曲線Fig.3 The temporal changes of concentrations at different time

表3 氮素遷移轉(zhuǎn)化參數(shù)Tab.3 N transport and transformation parameters for the soil profile

3.4 模型校正及檢驗(yàn)

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

圖4 實(shí)測(cè)硝態(tài)氮的濃度與模擬值對(duì)比曲線Fig. 4 Comparison of the measured and

4 不同灌溉強(qiáng)度硝態(tài)氮濃度的變化特征

常德市雖然雨量充沛，但降雨較集中，缺水季節(jié)需要對(duì)農(nóng)作物進(jìn)行灌溉。利用校正后的模型對(duì)不同灌溉強(qiáng)度下硝態(tài)氮濃度變化進(jìn)行模擬。該地區(qū)土壤滲透系數(shù)約為 4 cm/d，將灌溉強(qiáng)度分別設(shè)為1，2，4，8，12 cm/d，運(yùn)行模型。模擬結(jié)果如圖 5 所示。

圖5 不同灌溉強(qiáng)度硝態(tài)氮濃度Fig.5 The values of concentrations in different irrigation intensity scenarios

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

5 結(jié) 論

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