基于HYDRUS?1D模型的平原河網(wǎng)區(qū)農(nóng)田氮素遷移特征研究

許雪婷，班如雪，2，金有杰，張毅敏*，高月香，徐豪杰，2

（1.生態(tài)環(huán)境部南京環(huán)境科學(xué)研究所，南京 210042；2.河海大學(xué)環(huán)境學(xué)院，南京 210098；3.水利部南京水利水文自動化研究所，南京 210012）

太湖平原河網(wǎng)區(qū)農(nóng)業(yè)面源氮素流失已成為水質(zhì)惡化和水體富營養(yǎng)化的重要來源之一，太湖流域農(nóng)業(yè)面源污染物氨氮和總氮排放量分別占總量的43.3%和51.3%。其中，來自農(nóng)業(yè)活動的氮素輸入是導(dǎo)致太湖流域水質(zhì)惡化的一個重要因素，江蘇省2016年化肥施用強度為407.1 kg·hm，高于全國的平均水平359.1 kg·hm，接近國家制定的化肥污染的安全使用上限（225.1 kg·hm）的2倍。因此，選擇太湖流域平原河網(wǎng)農(nóng)田深入研究氮素流失特征具有重要意義。

基于此，本研究選擇HYDRUS?1D 模型進(jìn)行相關(guān)研究，該模型在一維土壤水動力模型基礎(chǔ)上，動態(tài)連接水分、溶質(zhì)、作物生長等多個模塊，充分考慮水分和溶質(zhì)在土壤包氣帶?地下水復(fù)合介質(zhì)中的遷移轉(zhuǎn)化過程，動態(tài)模擬氮素在土壤剖面的遷移轉(zhuǎn)化特征。本文以常州市武進(jìn)區(qū)為研究區(qū)域，該地區(qū)是典型的平原河網(wǎng)區(qū)，地勢平坦、河網(wǎng)密集，氣候溫和濕潤，雨量充沛，區(qū)域內(nèi)主要土地利用類型為農(nóng)業(yè)種植，除水域面積外，農(nóng)業(yè)種植面積占比高達(dá)60%。選擇區(qū)域典型的水稻土作為供試土壤，通過室外監(jiān)測和室內(nèi)淋濾實驗相結(jié)合的方式構(gòu)建本地化的HYDRUS?1D 水氮數(shù)值模型，實現(xiàn)對農(nóng)田土壤水氮運移分布模擬，并探究農(nóng)業(yè)面源污染氮素流失對水環(huán)境的影響，以期為太湖流域農(nóng)業(yè)面源污染控管提供科學(xué)依據(jù)。

1 材料與方法

1.1 室內(nèi)土柱淋濾試驗

1.1.1 樣品采集與分析

平原河網(wǎng)區(qū)土壤主要為水稻土，涉及土種包括黃泥土、白土、紅黃土和灰斑黃泥土等。本研究采樣點（120°6′50.64″E，31°46′32.66″N）位于常州市武進(jìn)區(qū)，選擇區(qū)域內(nèi)分布最為廣泛的黃泥土作為土壤樣品，借助GPS 全球定位系統(tǒng)，結(jié)合區(qū)域土地利用開展土壤取樣，設(shè)置3 個平行取樣點。每個取樣點在清除地表植被后，挖出1 m×1 m×0.5 m的正方形深坑，分別在0~10 cm 和10~20 cm 處取樣，利用環(huán)刀取原狀土用于測容重，利用鉛盒取土樣用于測土壤含水量，再取土樣若干于自封袋中用于土壤理化性質(zhì)測定和土柱制備。

將采集的土壤樣品自然風(fēng)干，去除植物殘體、根須及礫石等雜質(zhì)，混勻研磨后過1 mm 篩用于測定土壤背景值及分析土壤性質(zhì)。本文采用傳統(tǒng)的環(huán)刀法對土壤容重進(jìn)行測定，采用烘干法測定土壤的含水率，采用重鉻酸鉀容量法?外加熱法測定土壤有機質(zhì)。同時采用甲種比重計法對機械組成進(jìn)行測定。上述分析方法具體步驟參考《土壤分析技術(shù)指南》，測定結(jié)果見表1。

表1 土壤基本物理性質(zhì)Table 1 Basic physical properties of soil

1.1.2 實驗裝置與實驗設(shè)計

本研究設(shè)計的土柱淋濾裝置，其外徑8 cm，內(nèi)徑7.5 cm，高度55 cm，管內(nèi)土樣采用分層填土的方式裝填。裝填前，在柱子底部墊4 層紗布，再鋪上5 cm 經(jīng)鹽酸浸泡的石英砂作為承托層和過濾層，隨后依次在管內(nèi)鋪上10~20 cm 和0~10 cm 土層的土壤自然樣品各20 cm，0~10 cm土樣質(zhì)量為1 995.86 g，10~20 cm土樣質(zhì)量為2 172.49 g，并在土壤表層鋪5 cm 經(jīng)鹽酸浸泡過的石英砂。從石英砂層開始計算，每10 cm 設(shè)置一個取樣口，總共5個取樣口（圖1）。

圖1 土柱裝置示意圖Figure 1 The experimental device of indoor soil column leaching

土柱制備完成后，自上而下用純凈水濕潤土柱，充分淋洗。研究區(qū)年均降水量1 074 mm，灌水量300 m，施氮量在360~1 500 kg·hm，為模擬常規(guī)水肥條件下的氮素遷移過程，土柱模擬水量設(shè)置為400 mL，施用尿素總量設(shè)置為0.4 g。將稱取的尿素與土柱表層5 cm石英砂拌勻，回填進(jìn)土柱，灌水400 mL至土柱表面。按照設(shè)計的時間間隔（60、120、180、240、400、720 h），從土柱不同高度取土壤樣品及水樣，用于監(jiān)測土柱中氮素的剖面分布和淋失情況，并根據(jù)土柱實測數(shù)據(jù)對模型進(jìn)行校驗。

1.2 HYDRUS?1D模型構(gòu)建

HYDRUS?1D 軟件是由美國農(nóng)業(yè)部國家鹽土實驗室研發(fā)，可用于模擬飽和條件下水分、能量、溶質(zhì)在土壤中的遷移轉(zhuǎn)化規(guī)律和分布情況。該模型采用修改過的Richards方程為水流控制方程，可靈活設(shè)置并處理各邊界條件，且充分考慮植物根系吸水和土壤持水特性等多種因素，在模擬土壤水分、溶質(zhì)及農(nóng)藥的運移方面得到廣泛運用。自2000年引入我國，經(jīng)眾多學(xué)者開發(fā)應(yīng)用，模型得到不斷完善，廣泛應(yīng)用于飽和?非飽和帶污染物遷移相關(guān)領(lǐng)域的研究。

1.2.1 基本方程

（1）土壤水分運動基本方程

采用VanGenuchten 土壤水力模型模擬預(yù)測土壤水分運動，其模型如公式（1）~公式（4）所示：

式中：為土壤含水量，cm·cm；為壓力水頭，cm；為土壤飽和含水量，cm·cm；為土壤殘余含水量，cm·cm；為土壤水力傳導(dǎo)系數(shù)，cm·d；為有效飽和度，cm·cm；、、為經(jīng)驗系數(shù)，的值為進(jìn)氣值的倒數(shù)，cm；為孔隙分布指數(shù)；為土壤孔隙連通性參數(shù)，通常取值為0.5。

（2）土壤溶質(zhì)運動基本方程

HYDRUS?1D 的溶質(zhì)運移模型主要采用經(jīng)典的對流?彌散方程描述飽和非飽和孔隙介質(zhì)中的一維溶質(zhì)運移，見公式（5）：

式中：為溶質(zhì)的液相濃度，g·cm；為溶質(zhì)的固相濃度，g·g；為彌散系數(shù)，反映土壤水中有效分子擴散和機械彌散機制；為垂向水分通量，cm·s；為源匯項，g·cm·s。

式中：為土壤溶液中NO?N的濃度，mg·L；為源匯項，g·cm·s，主要包括有機氮的礦化、反硝化作用以及微生物固定等；為空間坐標(biāo)。

式中：、分別為氨氮的硝化速率和硝態(tài)氮的反硝化速率，h；為根系吸收硝態(tài)氮濃度，mg·L；為根系吸水量，cm·cm·h。

氨氮運移方程見公式（8）：

式中：為土壤溶液中氨氮的濃度，mg·L；為線性吸附平衡常數(shù)；表示土壤干容重，g·cm。

式中：為有機氮的礦化速率，h；為土壤有機氮含量，mg·L；為根系吸收硝態(tài)氮濃度，mg·L；為根系吸水量，cm·cm·h。

1.2.2初始條件及模型邊界

（1）土壤水分運動初始條件邊界與邊界條件

本次模擬水分初始條件見公式（10）：

式中：（）為灌水前不同深度土壤的初始含水量，cm·cm；（0，）為土柱的高度范圍，cm。

由于室內(nèi)土柱實驗?zāi)M的是包氣帶非飽和流的連續(xù)灌溉情況，因此水流模型的上邊界選用有表面積水的大氣邊界條件，其最大積水厚度為9 cm；下邊界選擇自由排水邊界條件。

（2）土壤溶質(zhì)運動初始條件與邊界條件

式中：c為初始氮素濃度的分布，mg·mL；為已知的土壤表層或隨入滲水下滲的溶質(zhì)濃度，mg·mL；為輸入氮肥的脈沖時間，h。

模型假設(shè)初始時刻土壤表層為飽和狀態(tài)，此時模型頂部的壓力水頭設(shè)置為0 cm，根據(jù)介質(zhì)剖分的節(jié)點數(shù)對土柱進(jìn)行離散，土柱底部的水力水頭設(shè)置為?50 cm。

1.2.3網(wǎng)格剖分和時間步長

根據(jù)室內(nèi)土柱淋濾實驗，均等劃分模擬土柱，并設(shè)置觀測點，模擬時長720 h，模擬的初始時間步長0.024 h，最小時間步長0.000 24 h，最大時間步長1 h。迭代控制參數(shù)使用默認(rèn)值。

1.2.4模型參數(shù)確定

（1）土壤水力參數(shù)

土壤水力參數(shù)是在土壤機械組成和容重基礎(chǔ)上，利用神經(jīng)網(wǎng)絡(luò)Rosetta模型進(jìn)行測算，再結(jié)合土柱淋濾實測數(shù)據(jù)校驗得出的，各土層土壤水力學(xué)數(shù)值見表2。

表2 土壤水力學(xué)參數(shù)Table 2 Soil hydraulic parameters

（2）土壤溶質(zhì)運移轉(zhuǎn)化參數(shù)

該模型涉及土壤溶質(zhì)運移轉(zhuǎn)化的多個參數(shù)，在參數(shù)敏感性檢驗和前人對土壤溶質(zhì)運移以及土壤溶質(zhì)轉(zhuǎn)化的研究的基礎(chǔ)上，利用土柱淋溶試驗測定并校驗，得出的模型土壤溶質(zhì)運移和轉(zhuǎn)化參數(shù)取值如表3和表4所示。

表3 土壤溶質(zhì)運移參數(shù)[14]Table 3 Soil nitrogen migration parameters[14]

表4 土壤溶質(zhì)轉(zhuǎn)化參數(shù)[23]Table 4 Soil nitrogen conversion parameters[23]

1.2.5模型校驗

本研究采用相關(guān)系數(shù)和均方根誤差（RMSE）評價指標(biāo)來定量評估模型模擬結(jié)果。計算公式如下：

式中：s和o為第個樣本的土壤含水量及氨氮、硝態(tài)氮濃度的模擬值和實測值；為樣本個數(shù)。

代表模型的擬合程度，數(shù)值越接近于1，表示模擬數(shù)值越接近實際情況。RMSE 反映模擬值與實測值之間的平均絕對誤差程度，數(shù)值范圍是0 到正無窮，數(shù)值越小，表示模擬越接近真實條件。

2 結(jié)果與討論

2.1 模型率定與驗證

模型率定中，選取兩種常規(guī)灌水施肥條件，即灌水量400 mL、施尿素0.4 g時和灌水量800 mL、施尿素0.4 g 時，15 cm 和35 cm 處土壤剖面土壤含水量、氨氮和硝態(tài)氮濃度監(jiān)測值與模擬值情況，結(jié)果如圖2、圖3、表5所示。

圖2 400 mL灌水條件下土壤含水量、氨氮及硝態(tài)氮濃度率定結(jié)果Figure 2 Comparison of measured values and simulated values of soil water content，ammonia nitrogen and nitrate nitrogen concentration under 400 mL irrigation conditions（fixed）

圖3 800 mL灌水條件下土壤含水量、氨氮及硝態(tài)氮濃度驗證結(jié)果Figure 3 Comparison of measured values and simulated values of soil water content，ammonia nitrogen and nitrate nitrogen concentration under 800 mL irrigation conditions（verification）

表5 兩種率定條件下模型模擬土壤含水量、水溶液中氨氮和硝態(tài)氮的評價結(jié)果Table 5 The evaluation results of soil moisture content，ammonia nitrogen and nitrate nitrogen in aqueous solution were simulated by the model under two calibration conditions

在灌水量為400 mL 的率定條件下，土壤剖面RMSE整體可控，均大于0.8，說明該模型結(jié)果較好，與實際情況較為一致，但RMSE＞RMSE，說明模型在深層的模擬效果更佳，不同形態(tài)氮素中，氨氮效果最好，硝態(tài)氮次之。

在灌水量達(dá)到800 mL的驗證條件下，RMSE均小于3，均大于0.85，說明模型模擬效果較好。綜合來看，模型邊界合適，參數(shù)選擇合理，在率定和驗證條件下，均大于0.8，模型模擬數(shù)據(jù)與實測數(shù)值較為一致，該模型可用于模擬平原河網(wǎng)區(qū)農(nóng)田包氣帶水氮運移。

2.2 農(nóng)田包氣帶氮素遷移轉(zhuǎn)化時空分布

土壤包氣帶中氮素主要包括有機氮、氨氮和硝態(tài)氮3種形態(tài)，圖4為不同時刻尿素態(tài)氮、硝態(tài)氮和氨氮隨土壤剖面遷移特征以及在不同剖面的濃度分布模擬結(jié)果。農(nóng)田施肥后，尿素會隨土壤深度增加而逐步遞減，在30 cm 土壤深度時濃度基本降解為0；而氨氮和硝態(tài)氮波動趨勢較為一致，均是隨著土壤深度的增加，濃度呈現(xiàn)先增加后緩慢減少的趨勢，其中氨氮在120 h 時達(dá)到濃度峰值，此時濃度為37.56 mg·L；而硝態(tài)氮在240 h 時可達(dá)到濃度最高值，為51.02 mg·L。這與氮素在土壤中的遷移轉(zhuǎn)化相關(guān)，尿素施用初期大量水解為氨氮，尿素態(tài)氮濃度下降而氨氮濃度持續(xù)增加，同時伴隨著硝化作用的發(fā)生，水解的氨氮在硝化細(xì)菌的作用下進(jìn)一步轉(zhuǎn)化為硝態(tài)氮，導(dǎo)致硝態(tài)氮的濃度增加，然而受到平原河網(wǎng)地下水位較淺的影響，硝態(tài)氮極易向深層土層遷移，污染地下水。

初始條件下，充分淋溶土柱直至檢測不出氮素，此時尿素態(tài)氮、氨氮和硝態(tài)氮濃度均為0。由圖4 氮素剖面分布情況可知，隨著尿素的加入，農(nóng)田中各形態(tài)氮素發(fā)生動態(tài)轉(zhuǎn)化。首先是尿素發(fā)生水解轉(zhuǎn)化為氨氮，由于尿素態(tài)氮是有機氮的一種形態(tài)，遷移能力較弱，該過程主要發(fā)生在0~5 cm 土壤表層，且反應(yīng)速率較快，400 h 后尿素態(tài)氮含量接近于0，與此同時氨氮濃度不斷增加，在180 h 時達(dá)到最高值。土壤剖面5、15、25 cm 和35 cm 處氨氮平均濃度依次為：10.97、8.94、2.54 mg·L和0.51 mg·L，由此可知，氨氮轉(zhuǎn)化過程主要集中在0~15 cm 土層。由于土壤膠粒與硝態(tài)氮均帶負(fù)電荷，土壤中的硝態(tài)氮易隨土壤水分遷移，硝態(tài)氮在5、15、25 cm 和35 cm 土層處平均濃度依次為：27.31、27.55、15.89 mg·L和5.88 mg·L，相較于其他形態(tài)氮素，硝態(tài)氮在整個土壤剖面中的濃度均較高，在15 cm 處濃度達(dá)到峰值，并隨著土壤水分逐漸向深層土壤中遷移。研究表明，以硝酸鹽形式存在的溶解性氮是地下水環(huán)境中最為常見的污染物，尤其在平原河網(wǎng)區(qū)地下水水位埋藏較淺，平均水位在50~300 cm 間，不考慮其他外源污染輸入等條件，15 cm 以下的土壤剖面每增加10 cm 硝態(tài)氮會減少10 mg·L，由此可見，平原河網(wǎng)區(qū)硝態(tài)氮極易進(jìn)入地下水污染水體。

圖4 不同時刻尿素態(tài)氮、氨氮和硝態(tài)氮遷移及剖面分布特征Figure 4 Migration and profile distribution characteristics of urea nitrogen，ammonia nitrogen and nitrate nitrogen at different times

2.3 農(nóng)業(yè)灌溉對氮素遷移轉(zhuǎn)化的影響

2.3.1 不同形態(tài)氮素在土壤剖面的濃度分布

本研究設(shè)置了200、400、600 mL 和800 mL 4 種農(nóng)業(yè)灌溉情景，探究農(nóng)業(yè)灌溉對不同形態(tài)氮素遷移轉(zhuǎn)化和流失的影響。總體來看，隨著灌溉量的增加，尿素、氨氮和硝態(tài)氮在各土壤剖面的波動趨勢一致（圖5）。其中，尿素態(tài)氮濃度隨土壤深度增加而遞減，其降解主要集中在表層土壤中，濃度最高達(dá)到40 mg·L，10~30 cm 土層中尿素急劇減少，到40 cm 處濃度基本為0。氨氮和硝態(tài)氮濃度都是呈現(xiàn)先增大后減少的趨勢，10 cm處氨氮濃度基本可達(dá)到最高值，200、400、600 mL 和800 mL 灌溉量下最高氨氮濃度分別為：10.02、13.54、15.64 mg·L和18.67 mg·L。伴隨著硝化反應(yīng)，氨氮濃度逐漸降低，同時硝態(tài)氮濃度迅速上升，10~20 cm 處硝態(tài)氮濃度達(dá)到峰值，200、400、600 mL 和800 mL 灌溉量下最高硝態(tài)氮濃度分別為：20.27、25.11、27.62 mg·L和28.37 mg·L。隨后硝態(tài)氮的濃度緩慢降低，在50 cm 處，當(dāng)灌溉量為200 mL時，硝態(tài)氮的濃度接近0，而灌溉量為400、600 mL 和800 mL時，硝態(tài)氮的濃度依次為1.78、6.41、15.27 mg·L?？傮w來看，隨著灌溉量的增加，各形態(tài)氮素流失量增加，受電荷影響，硝態(tài)氮較氨氮更易隨土壤水分向土壤深層遷移。但當(dāng)灌溉量高于400 mL 時，氨氮濃度在40~50 cm 深層土壤中出現(xiàn)增加趨勢，這是由于隨著剖面中硝態(tài)氮所占比例的升高，氨氮轉(zhuǎn)為硝態(tài)氮的速率降低，且深層土壤中氧氣不足，硝化反應(yīng)受到抑制，隨著灌溉量的增加，抑制作用愈加明顯，因此當(dāng)灌溉量高于400 mL，且深層土壤中硝態(tài)氮與氨氮比例達(dá)到9.8∶1時，土壤中的硝化作用受到抑制。

圖5 不同灌水條件下各形態(tài)氮素遷移轉(zhuǎn)化分布Figure 5 The average nitrogen concentration of the profile changes with time under different irrigation conditions

2.3.2 不同形態(tài)氮素在土壤剖面的累積淋失量

不同灌水條件下各形態(tài)氮素累積淋失量見圖6。尿素態(tài)氮在表層土壤中完成降解，僅在800 mL 灌溉條件下有少量流失，其值小于5 mg·cm，而氨氮和硝態(tài)氮隨土壤水分存在明顯垂向遷移的趨勢，尤其是硝態(tài)氮流失量較大，在200、400、600 mL 和800 mL 灌溉量條件下累積淋失量分別為0.08、3.18、35.27、142.89 mg·cm，相同條件下氨氮的累積流失量分別為0.01、0.71、9.61、35.33 mg·cm。由此可見，氮素在土壤剖面中主要以硝態(tài)氮和氨氮兩種形態(tài)流失，且隨著灌溉量的增加，流失量呈指數(shù)增加。灌溉量由200 mL 增加至800 mL 時，硝態(tài)氮流失量增加了142.81 mg·cm，增加了1 785 倍；氨氮流失量增加了35.32 mg·cm，增加倍數(shù)高達(dá)3 532 倍。相關(guān)研究表明，硝態(tài)氮是地下水污染防控重點關(guān)注的污染物質(zhì)，但土壤剖面中氨氮流失對灌溉量更為敏感，同樣也是地下水環(huán)境中需要關(guān)注的污染物質(zhì)。

圖6 不同灌水條件下各形態(tài)氮素累積淋失量Figure 6 Accumulated leaching amount of different forms of nitrogen under different irrigation conditions

平原河網(wǎng)區(qū)地下水水位埋藏較淺，平均水位在50~300 cm，根據(jù)《地下水質(zhì)量標(biāo)準(zhǔn)》（GB/T 14848—2017）中的Ⅲ類飲用水標(biāo)準(zhǔn)，硝酸鹽（NO）含量不超過20 mg·L（即硝態(tài)氮含量不超過4.5 mg·L），氨氮含量不超過0.5 mg·L。研究區(qū)域常規(guī)灌溉量為400 mL，在此灌溉條件下硝態(tài)氮和氨氮流失濃度分別為2.84 mg·L和0.34 mg·L，在不考慮土壤背景值和其他外源污染輸入的情況下均未超過Ⅲ類水質(zhì)標(biāo)準(zhǔn)。但區(qū)域降雨等因素會導(dǎo)致區(qū)域水量遠(yuǎn)高于400 mL，在800 mL 灌溉量時，硝態(tài)氮和氨氮的濃度分別達(dá)到30.62 mg·L和6.12 mg·L，分別超過標(biāo)準(zhǔn)5.8 倍和11.2 倍，綜合考慮實際情況，平原河網(wǎng)區(qū)進(jìn)入地下水水環(huán)境中的氮素將遠(yuǎn)超于Ⅲ類水質(zhì)標(biāo)準(zhǔn)。

3 結(jié)論

（1）本研究利用HYDRUS?1D 構(gòu)建了太湖流域平原河網(wǎng)區(qū)水氮運移模型，模型的RMSE 可控，大于0.85，表明模擬值與實測值擬合度較好，可應(yīng)用于區(qū)域土壤包氣帶水氮遷移轉(zhuǎn)化模擬。

（2）研究區(qū)常規(guī)灌水施肥條件下，隨土壤剖面的加深，尿素態(tài)氮濃度遞減，氨氮和硝態(tài)氮濃度呈現(xiàn)先增大后減小的變化趨勢。尿素態(tài)氮主要分布在0~5 cm 處，氨氮主要分布在0~15 cm 處，硝態(tài)氮主要分布在0~35 cm處。

（3）土壤包氣帶中氨氮流失對灌溉量更為敏感。當(dāng)灌溉量高于400 mL 時，氮素流失量呈指數(shù)增加趨勢，當(dāng)深層土壤中硝態(tài)氮與氨氮比例達(dá)到9.8∶1時，土壤中硝化作用受到抑制。從流失量和流失機理來看，硝態(tài)氮更易進(jìn)入土壤深層，污染地下水環(huán)境?？紤]到自然降雨、外源污染輸入和土壤背景值，平原河網(wǎng)區(qū)進(jìn)入地下水環(huán)境中的氮素將遠(yuǎn)超于Ⅲ類水質(zhì)標(biāo)準(zhǔn)。