滴灌施肥條件下土壤水氮運移數(shù)值模擬

黎會仙，王文娥，胡笑濤

(西北農林科技大學旱區(qū)農業(yè)水土工程教育部重點實驗室,陜西 楊凌 712100)

水、肥是影響作物生長的2個重要因子，在實際生產中如果灌溉和施肥方式不合理，將造成水、肥浪費及地下水污染。滴灌施肥是一種節(jié)水節(jié)肥的先進技術[1]，目前滴灌隨水施肥過程中土壤水分、養(yǎng)分的運移分布規(guī)律大部分是基于試驗獲取，而試驗方法費時費力。計算機數(shù)值模擬具有簡便、靈活等特點，能實現(xiàn)對各種灌溉情況的預測[2]，建立適當?shù)臄?shù)學模型可以更高效準確地研究滴灌施肥條件下水分、養(yǎng)分的運移轉化過程。近年來，相關學者建立了不同滴灌條件下土壤水分的運移模型[3-4]，并且滴灌條件下土壤水和溶質運移模型都已取得發(fā)展[5-6]。脫云飛等[7]建立了膜孔肥液自由入滲的土壤水氮運移模型，較好地反映了土壤水氮運移轉化規(guī)律。有學者利用 HYDRUS-1D 軟件對不同情景下土壤水氮的遷移變化進行動態(tài)模擬，均獲得了擬合度較高的水氮模擬模型[8-9]。EU-Rotate_N和HYDRUS-2D 模型均能模擬氮的淋溶損失，較精確地預測氮含量，在模擬氮的吸收及運移方面各自具有優(yōu)勢[10]。Tournebize等[11]應用SWMS模型模擬了植被覆蓋條件下氮的運移。G?rden?s等[12]模擬了不同微灌系統(tǒng)、灌溉施肥措施和土壤類型對硝酸鹽淋失的影響。熱合木等[13]模擬了膜下滴灌土壤水、鹽分運移，獲得較好的效果。

以上研究者應用不同的模型對土壤水、溶質在不同灌溉施肥條件下的運移進行研究，均獲得較好的效果，但現(xiàn)有的土壤水、溶質運移研究大多僅考慮單一溶質的運移。滴灌施肥條件下土壤氮的動力學過程非常復雜，與土壤溫度、pH、含水率、微生物、土壤類型等多種因素有關。

尿素是廣泛使用的肥料，且在土壤中易分解轉化，由于尿素態(tài)氮在土壤中運移轉化的復雜性，目前對滴灌施肥條件下土壤尿素態(tài)氮的運移轉化模擬較少，因此，有必要對尿素態(tài)氮在土壤中的運移變化進行模擬計算，以便對滴灌系統(tǒng)水肥的運移進行預測。本研究基于試驗基礎，利用HYDRUS軟件對滴灌施肥條件下土壤水氮(尿素)運移及再分布過程進行模擬分析，并利用滴灌水肥入滲試驗對模型進行驗證，以確定模型對滴灌施肥條件下土壤水氮運移分布模擬的精度，為滴灌施肥系統(tǒng)的優(yōu)化設計提供參考依據(jù)，提高水肥利用效率，促進增產增收。

1 試驗設計與方法

1.1 土壤物理性質

供試土壤取自西北農林科技大學試驗田0～40 cm土層，土壤基本物理參數(shù)見表1。

1.2 試驗裝置與方法

試驗裝置包括馬氏瓶、圓土柱(高50 cm，半徑35 cm)、滴頭、EM50傳感器及數(shù)據(jù)采集器，如圖1所示。試驗所用肥料為普通尿素(CON2H4)。

試驗設置滴頭流量分別為 2、3、4 L·h-1，灌水量為21 L，施肥濃度為300 mg·L-1。試驗前將風干土過2 mm篩，土壤初始體積含水率為0.057 cm3·cm-3，初始土壤銨態(tài)氮、硝態(tài)氮含量分別為6.5、25.05 mg·kg-1。試驗土柱分層裝土，裝土時進行攪拌，各層之間進行打毛、平整，以保證土柱內土壤分布均勻，各向同性均質。裝土過程中埋設傳感器，每個土柱埋設12個EM50傳感器，分層布置在土體中，根據(jù)滴灌土壤水分入滲特點[14-15]及傳感器測試精度和范圍，傳感器的布置方式見圖1。裝土完畢后靜置24 h，使土體初始含水率均勻穩(wěn)定。灌水前在土柱表面覆蓋一層薄膜。試驗前對所有儀器進行標定。

表1 土壤基本物理特性

注:圖中數(shù)字所示為傳感器分布位置。Note: The numbers in the figure show the distributed location of sensor.圖1 試驗裝置及傳感器布置Fig.1 Experimental equipment and sensor placement

試驗開始后通過數(shù)據(jù)采集器監(jiān)測土壤含水率隨時間的變化情況，滴灌過程中每5 min監(jiān)測一次，滴灌結束后每10 min監(jiān)測一次。在滴灌結束24、48、72、96、120 h分別在5個縱向半剖面上用土鉆按傳感器布置的位置點取土，取回的土樣用紫外分光儀測定土壤硝態(tài)氮、銨態(tài)氮含量。

2 滴灌施肥土壤水氮運移基本方程及模型

2.1 土壤水分運動基本方程

根據(jù)達西定律和質量守恒定律[16]及土壤水動力學原理，假定土壤均質和各向同性，忽略空氣、溫度和滯后效應對水分運動的影響，滴灌三維土壤水分入滲過程可用Richard方程來描述：

(1)

式中,φm為基質勢(cm)；k(θ)為非飽和土壤導水率(cm·h-1)；θ為土壤體積含水率(cm3·cm-3)；t為時間(h)。

2.2 土壤氮素運移基本方程

尿素態(tài)氮運移方程

(2)

銨態(tài)氮運移方程

+k0ρ+k1θC1-k2θC2-k2ρs

(3)

硝態(tài)氮運移方程

+k2θC2+k2ρs-k3θC3

(4)

式中,x為空間坐標(角標i=1，2，3，j=1，2，3，表示空間3個方向)；C1、C2、C3分別為尿素態(tài)氮、銨態(tài)氮、硝態(tài)氮的質量濃度(mg·cm-3)；s為固相氨態(tài)氮吸附量，s=kdC2，kd為銨態(tài)氮吸附速率(cm3·g-1)；qi為土壤水通量(cm·h-1)；k0為有機質礦化速率(h-1)；k1為尿素水解速率(h-1)；k2為銨態(tài)氮硝化速率(h-1)；k3為硝態(tài)氮反硝化速率(h-1)；ρ為土壤干容重(mg·cm-3)；Dij為水動力彌散系數(shù)張量的分量(cm2·h-1)，包含縱向、橫向彌散系數(shù)DL和Dr(cm)，以及自由水中的分子擴散系數(shù)DW(cm2·h-1)。運移及轉化參數(shù)取值見表2?；谝陨贤寥浪?、溶質動力學原理，利用HYDRUS軟件建立幾何模型并進行數(shù)值計算。

2.3 土壤水力參數(shù)

本文土壤水分特征曲線θ(h)利用高速離心機測定，土壤飽和導水率K(h)用定水頭法進行測定，土壤水分特征曲線采用Van Genuchten模型擬合[17-18]。

(5)

(6)

其中

Se=(θ-θr)/(θs-θr)

(7)

m=1-1/n,n>1

(8)

式中,θs為土壤飽和導水率(cm3·cm-3)；θr為土壤殘余含水率(cm3·cm-3)；α,n,m為土壤水分特征曲線實測值(土壤含水率、土壤水吸力)擬合得到的參數(shù)(cm-1)；h為土壤負壓水頭(cm)；Ks為滲透系數(shù)(cm·min-1)；Se為土壤有效含水率(cm3·cm-3)；l為空隙連通性參數(shù)，一般情況取為0.5。擬合得到土壤水分特性的參數(shù)值為θr=0.03 cm3·cm-3；θs=0.48 cm3·cm-3；α=0.036 cm-1；n=1.56；Ks=0.125 cm·min-1。

2.4 土壤溶質運移參數(shù)

模型中氮素運移轉化參數(shù)的確定根據(jù)國內外文獻[19-23]，并根據(jù)實測試驗數(shù)據(jù)以流量分別為2、3 L·h-1對各參數(shù)進行率定確定，模型中氮素轉化參數(shù)的取值見表2。

2.5 定解條件

2.5.1 初始條件 假定土壤初始含水率和氮濃度在研究區(qū)域內均勻分布，初始條件可表示為：

θ(x,y,z,t)=θ0
(0≤x≤X，0≤y≤Y，0≤z≤Z，t=0)

(9)

C(x,y,z,t)=C0
(0≤x≤X，0≤y≤Y，0≤z≤Z，t=0)

(10)

式中,θ和C分別為土壤含水率，土壤銨態(tài)氮或硝態(tài)氮質量濃度(mg·cm-3)；θ0和C0分別為初始土壤含水率,初始土壤銨態(tài)氮或硝態(tài)氮質量濃度(mg·cm-3)；t為時間(h)，初始時刻時間為0 h。X，Y和Z為模擬區(qū)域邊界(裝置物理邊界)的坐標。

2.5.2 邊界條件 根據(jù)試驗設計模擬研究區(qū)域為圓柱體(高50 cm，半徑35 cm)，滴頭在圓柱軸心的正上方，取圓柱體的一半為模擬研究區(qū)域。根據(jù)試驗監(jiān)測，滴灌開始時，土壤飽和區(qū)發(fā)展很快，隨著時間的推移，飽和區(qū)擴展的速度逐漸減小，一定時間后飽和區(qū)半徑(Rs)趨于定值為7.3 cm，所以模擬時忽略飽和區(qū)半徑隨時間的變化。飽和區(qū)在滴灌過程中為不隨時間變化的定流量邊界，不灌水時為零通量邊界。上邊界其余覆膜區(qū)域為零通量邊界。故模擬區(qū)域內將上邊界分為兩種情況：

(1)飽和區(qū)流量和溶質濃度邊界條件為：

(11)

(12)

(2)非飽和區(qū)流量和溶質濃度邊界條件為：

(13)

(14)

式中,Ca為肥料溶液尿素態(tài)氮的質量濃度(mg·cm-3)；T為灌水歷時(h)；C1，xi，Dij,qi與上式(4)中的含義相同；K(θ)為非飽和土壤導水率(cm·h-1)；σ(t)為灌水期間滴頭處定流量邊界通量(cm·h-1)；t>T時，σ(t)=0，z為豎軸方向的坐標。

模擬計算區(qū)域邊壁為零通量邊界；下邊界為自由排水邊界。

2.6 數(shù)值求解

利用HYDRUS軟件對上述模型進行求解，用Galerkin有限單元法對數(shù)學模型進行數(shù)值計算[24]；時間采用隱式差分法離散[25]，空間離散采用伽遼金有限單元法。計算區(qū)內網(wǎng)格劃分采用三棱柱單元，滴頭飽和區(qū)半徑內適當加密網(wǎng)格，以提高計算精度。根據(jù)試驗設置，在模擬區(qū)設置與試驗布置相對應的12個觀測點，以便于對模擬結果進行驗證，以滴頭正下方土表為坐標原點，觀測點坐標分別為(水平，縱向，垂向)：點1(5,0,-5)，點2(15,0,-5)，點3(25,0,-5)，點4(35,0,-5)，點5(5,0,-15)，點6(15,0,-15)，點7(25,0,-15)，點8(5,0,-25)，點9(15,0,-25)，點10(25,0,-25)，點11(5,0,-35)，點12(15,0,-35)。模擬的網(wǎng)格劃分及觀測點設置如圖2所示半圓柱體。

表2 土壤尿素態(tài)氮運移轉化參數(shù)

3 模型驗證及結果分析

3.1 土壤含水率分布規(guī)律

土壤含水率的分布對作物的生長有直接影響，并影響土壤溶質等因素的變化，基于試驗研究，用數(shù)值模擬的方法對滴灌及其再分布過程土壤含水率分布進行分析，以流量4 L·h-1土體剖面內各觀測點土壤含水率對模型進行驗證，對比了滴頭流量為4 L·h-1時各觀測點的模擬值與實測值。模擬與實測的土壤含水率分布用Kriging插值法進行插值計算，在Surfer軟件中繪制。由圖3可知不同時刻模擬與實測值的土壤含水率剖面分布基本一致。滴灌結束時，土壤含水率實測均值與模擬值在土層5 cm相對誤差為4%，土層15 cm相對誤差為10%，土層25 cm相對誤差為7%。再分布過程，120 h部分觀測點的相對誤差較大。

注:圖中數(shù)字所示為觀測點的分布位置。Note: The numbers in the figure show the distribution of observation points.圖2 模擬研究區(qū)域網(wǎng)格劃分及觀測點布置Fig.2 Flow domain gridding and layout of observation points in simulation region

注：圖中土壤含水率單位為cm3·cm-3。Note: The unit of soil moisture in the figure is cm3·cm-3.圖3 土壤含水率分布剖面模擬與實測對比Fig.3 Comparison of simulation and measured values of soil water content

在濕潤土體內，滴灌結束時土壤含水率的分布隨著離滴頭距離的增加而減??；經再分布72 h，土體內水分向下層運移，土層25～30 cm的土壤含水率增大到0.2 cm3·cm-3。120 h土體剖面內的土壤含水率均有所下降，較滴灌結束時下降了18%，且土壤剖面內的土壤含水率分布趨于均勻穩(wěn)定。對滴灌后土壤含水率的再分布過程模擬與實測值均能較好地吻合，各觀測點相對誤差的絕對值均在10%以內，模擬能較好地反映土壤含水率的分布規(guī)律，可用數(shù)值模擬的方法對不同滴灌系統(tǒng)條件下的土壤含水率分布情況進行預測。

3.2 土壤氮分布特征

氮肥是植物生長所必需的營養(yǎng)元素，有機態(tài)氮不能被植物直接吸收利用，銨態(tài)氮、硝態(tài)氮是植物能直接吸收利用的營養(yǎng)離子[27-28];尿素是目前廣泛使用的氮肥，在土壤中能直接水解轉化為銨態(tài)氮、硝態(tài)氮，研究滴灌施肥后土壤硝態(tài)氮、銨態(tài)氮量的時空分布有助于確定不同作物氮肥的次施用量及濃度。

3.2.1 土壤銨態(tài)氮時空分布特征 土壤銨態(tài)氮在土體中易被吸附，在土壤剖面中的分布決定于作物生長、灌溉方式、土壤性質等因素[29]。基于試驗與模擬，以滴頭流量4 L·h-1的土壤銨態(tài)氮動態(tài)分布進行模擬值與實測值的對比分析，模型在與試驗布置相同位置坐標點處分別設置12個觀測點，根據(jù)觀測點的布置選取了有代表性的3個觀測點(點1、6、8)進行分析(觀測點位置見圖2)。

由圖4可知不同觀測點由于硝化作用土壤銨態(tài)氮質量濃度均隨時間的推移逐漸減小，且模擬值與實測值較接近。對滴灌結束后24～120 h的土壤銨態(tài)氮質量濃度模擬值與實測值進行對比分析，距滴頭較近的觀測點1(5，-5)土壤銨態(tài)氮質量濃度在土體中較大，72 h之前隨時間的變化幅度較大，隨后變化較小，其24、48、72、96、120 h各模擬值與實測值的相對誤差分別為17%、20%、20%、17%、14%；觀測點6(15,-15)、點8(5,-25)的銨態(tài)氮質量濃度在土體中隨時間的變化幅度較小，24～120 h內模擬值與實測值的相對誤差變化范圍分別為14%～20%，12%～20%。觀測點5(5,-15)的相對誤差變化范圍為1%～14%，對其它各觀測點進行統(tǒng)計分析其相對誤差最小值為2%，均在20%以內。由于實際土壤中微生物的作用不同，且土壤環(huán)境復雜，同一觀測點處不同時間土壤銨態(tài)氮質量濃度模擬值與實測值相對誤差存在差異。

圖5見16頁中土壤銨態(tài)氮質量濃度空間分布表明：滴灌后土壤銨態(tài)氮質量濃度主要分布于距滴頭20 cm的范圍內，隨著離滴頭距離的增加而逐漸減小，24 h土壤銨態(tài)氮質量濃度最大，到120 h時減少了40%。滴灌后不同觀測點的銨態(tài)氮濃度不同，且硝化作用存在差異，120 h部分觀測點的土壤銨態(tài)氮質量濃度低于本底值。模擬值與試驗值的相對誤差在合理范圍內，土壤銨態(tài)氮分布結果與Breve[23]模擬的結果相近，可用該模型模擬土壤銨態(tài)氮的動態(tài)分布。

圖4 不同觀測點滴灌后24～120 h土壤銨態(tài)氮質量濃度變化Fig.4 The change process of concentration at different observation points 24～120 h after drip irrigation

3.2.2 土壤硝態(tài)氮時空分布特征 滴灌施肥后土壤硝態(tài)氮分布因灌溉方式、施肥種類等多種因素而存在差異。結合試驗，利用數(shù)值模型對滴灌施尿素后土壤硝態(tài)氮的時空分布進行模擬，并對比分析實測值與模擬值。由圖6可知，不同觀測點土壤硝態(tài)氮質量濃度均隨著時間的推移而逐漸增大，且模擬值與實測值較吻合。觀測點1(5，-5)土壤硝態(tài)氮質量濃度隨時間的增長較快，72 h之后變化較小，24、48、72、96 h和120 h內土壤硝態(tài)氮質量濃度模擬值與實測值的最大相對誤差分別為18%、19%、16%、19%、20%；觀測點6(15,-15)、點8(5,-25)隨時間的變化幅度較小，24～120 h土壤硝態(tài)氮模擬值與實測值的相對誤差變化范圍分別為1%～20%，19%～20%。觀測點5(5,-15)模擬值與實測值的相對誤差變化范圍為14%～20%，對其它各觀測點進行統(tǒng)計分析其相對誤差最小值為13%，均在20%以內。

24 h時觀測點1，3，6，8，5的土壤硝態(tài)氮質量濃度分別為：0.437，0.467，0.451，0.482 mg·cm-3和0.447 mg·cm-3，相互間的差值小于0.05 mg·cm-3，各觀測點在空間分布上差異不大，48 h之后土體中觀測點1(5，-5)的土壤硝態(tài)氮質量濃度較大，120 h時為1.2 mg·cm-3。土壤銨態(tài)氮在硝化作用下逐漸轉化為土壤硝態(tài)氮，由于滴頭附近的土壤銨態(tài)氮質量濃度較大，土壤硝態(tài)氮質量濃度在滴頭附近增加較多。由圖6和圖7可知：滴灌施肥后土壤硝態(tài)氮在整個土體內均有分布，48 h后土壤硝態(tài)氮質量濃度在離滴頭20 cm的范圍值較大，且空間分布上隨距滴頭距離的增加而減小，垂直方向上從距滴頭5 cm的觀測點1(5，-5)到距滴頭25 cm的觀測點8(5，-25)減少了53%。該空間分布特征與Breve[23]模擬的結果相近。土壤硝態(tài)氮的模擬值與實測值不同觀測點、不同時間其相對誤差存在差異，但其平均相對誤差在18%以內，基本反映出土壤硝態(tài)氮分布的動態(tài)變化過程，可用該模型模擬滴灌施肥土壤硝態(tài)氮的動態(tài)分布。

圖6 不同觀測點滴灌后24～120 h土壤硝態(tài)氮質量濃度變化Fig.6 The change process of soil nitrate nitrogen concentration at different observation points 24～120 h after drip irrigation

圖5 滴灌結束后72 h、120 h土壤銨態(tài)氮質量濃度分布Fig.5 The concentration distribution of in soil 72 h and 120 h after drip irrigation

圖7 滴灌結束后72 h、120 h土壤硝態(tài)氮質量濃度分布Fig.7 The concentration distribution of nitrate nitrogen in soil 72 h and 120 h after drip irrigation

4 結論與討論

尿素施入土壤后，會以分子形態(tài)被土壤膠體吸附, 但在合適的水熱條件和豐富的有機物料供應情況下，會很快水解[30]。李世清[31]研究表明溫度、土壤含水率、尿素起始濃度、土壤質地、外加脲酶對尿素水解速率均有影響；汪建飛[32]研究了培養(yǎng)時間、尿素濃度、土壤溫度、土壤含水量對土壤外源尿素(施入土壤中的尿素)水解速率的影響，其影響大小順序為：培養(yǎng)時間>尿素濃度>土壤溫度>土壤含水量。樊慶魯[33]指出在pH>5.6 的土壤中, 溫度為 25℃時，尿素3 d可完全分解。王旭洋[34]研究了滴灌施肥(尿素)條件下速效氮在土壤中的運移，表明滴灌結束1 d，尿素水解，銨態(tài)氮含量達到最大值，隨后隨時間的推移逐漸減?。缓蚣t雨[35]研究表明尿素施入土壤6～24 h內，剖面中尿素的含量普遍降低了約66.7%。以上研究表明施入土壤中的尿素在不同的條件下，水解轉化為銨態(tài)氮和硝態(tài)氮的速率不同。

本文采用試驗與數(shù)值模擬的方法對滴灌施肥后土壤水肥的運移變化進行研究，得出滴灌結束24 h，尿素水解，銨態(tài)氮含量達到最大值。數(shù)值模擬中根據(jù)試驗設計設置參數(shù)，并對參數(shù)值率定，對比分析了模擬與試驗結果，模型能夠反映土壤水氮運移轉化過程，模擬所得土壤含水率剖面分布與實測值基本一致，各觀測點的土壤硝態(tài)氮、銨態(tài)氮質量濃度的模擬值與實測值相對誤差均在20%以內，不同觀測點、不同時間的相對誤差均不同。滴灌結束時，土壤含水率隨離滴頭距離的增加而減小。滴灌施肥后土壤銨態(tài)氮質量濃度主要分布于離滴頭20 cm的范圍，且隨著時間的推移逐漸減小；24～120 h土壤硝態(tài)氮質量濃度各觀測點隨著時間的推移逐漸增大，且硝態(tài)氮質量濃度在滴頭20 cm的范圍內由0.442 mg·cm-3增加到1.2 mg·cm-3，48 h后土體內土壤硝態(tài)氮質量濃度空間分布上隨著離滴頭距離的增加呈減小的趨勢。由于實際土壤中尿素態(tài)氮的轉化受多種因素的影響，且文中數(shù)值模型進行了簡化假設，忽略了銨態(tài)氮的揮發(fā)，土壤氮素的模擬值與實測值存在一定誤差，該模型還有待進一步研究完善。