基于Hydrus-2D的紅壤區(qū)涌泉根灌自由入滲土壤水分運移數(shù)值模擬

代智光

(南昌工程學(xué)院水利與生態(tài)工程學(xué)院，江西 南昌 330099)

紅壤主要分布在我國南方丘陵區(qū)，是我國主要土壤資源之一。紅壤區(qū)光熱資源充足，氣候條件優(yōu)越，特別適合柑、橘等果樹的生長。隨著農(nóng)業(yè)的發(fā)展，紅壤區(qū)季節(jié)性干旱以及紅壤酸黏貧瘠等問題日益突出，人們往往通過過量灌溉來提高果樹產(chǎn)量，這不僅導(dǎo)致灌溉水利用效率低下，而且造成當(dāng)?shù)赝寥狼治g和養(yǎng)分流失問題突出，發(fā)展科學(xué)合理的灌溉技術(shù)迫在眉睫。

涌泉根灌是一種微灌新技術(shù)，類似于地下滴灌，它通過埋在地下的灌水器將水分直接輸送到作物根區(qū)進(jìn)行局部灌溉，提高了水分利用效率。相比其它灌溉方式，涌泉根灌出流方式為三維柱狀出流，與土壤接觸面積大，特別適宜于紅壤等質(zhì)地較為黏重的土壤。近些年，相關(guān)學(xué)者對涌泉根灌技術(shù)進(jìn)行了大量研究，如Dai等[1-2]通過棗樹的涌泉根灌水氮耦合試驗研究表明，涌泉根灌能有效提高果樹的水肥利用效率，特別適宜于果樹灌溉；劉風(fēng)華等[3]研究了容重對紅壤區(qū)涌泉根灌土壤水分的運移特性影響，結(jié)果表明容重不同土壤水分運移特性存在顯著差異；代智光等[4]研究肥液濃度對紅壤區(qū)涌泉根灌水氮運移特性的影響，并建立了土壤累計入滲量及濕潤鋒運移距離與肥液濃度的關(guān)系模型；費良軍等[5]研究了黃土高原地區(qū)不同肥液濃度對涌泉根灌土壤濕潤體的影響，揭示了涌泉根灌土壤入滲量和濕潤鋒運移距離隨肥液濃度的變化關(guān)系；吳恒卿等[6]對雙點源涌泉根灌交匯入滲土壤濕潤體特性進(jìn)行了試驗研究，結(jié)果表明灌水器間距對濕潤體的運移、交匯時間、土壤含水率有明顯影響。近些年，國內(nèi)外對土壤水分運移的數(shù)值模擬開展了大量研究工作，如美國鹽分實驗室開發(fā)的Hydrus1D/2D/3D軟件可以模擬飽和非飽和土壤水、熱、鹽、單組分和多組分溶質(zhì)運移轉(zhuǎn)化問題；李耀剛等[7]采用Hydrus-3D軟件模擬了陜北黃綿土涌泉根灌土壤水分的運移，結(jié)果較好；Kandelous等[8]利用Hydrus-2D軟件對地下滴灌的土壤水和溶質(zhì)運移進(jìn)行了模擬，精度滿足要求；Cai等[9]利用Hydrus-2D軟件對陶瓷根灌土壤水分運移進(jìn)行了模擬，結(jié)果滿足要求；李久生等[10]基于土壤水分運動的動力學(xué)方程和溶質(zhì)運移的對流-彌散方程，建立了不同土壤中地表滴灌水肥運移的Hydrus-2D模型，并對模型進(jìn)行了求解；趙新宇等[11]采用Hydrus-3D軟件對涌泉灌的水分運移進(jìn)行了模擬，結(jié)果較好。

目前，有關(guān)紅壤區(qū)涌泉根灌土壤水分運移的研究還主要集中在室內(nèi)試驗，局限性較強，而對土壤水分運移的模擬研究還少見報道?；诖耍疚耐ㄟ^紅壤區(qū)涌泉根灌土壤水分運動的數(shù)值模擬，并與室內(nèi)試驗結(jié)果相對比，驗證Hydrus軟件在紅壤區(qū)涌泉根灌中的適用性，以期為涌泉根灌技術(shù)在紅壤區(qū)的推廣應(yīng)用提供參考。

1 材料與方法

1.1 試驗土壤

試驗于2019年5月在南昌工程學(xué)院灌溉排水實驗室進(jìn)行，供試土壤取自江西省進(jìn)賢縣田間表層0～50 cm熟土，土壤經(jīng)自然風(fēng)干后進(jìn)行研磨，過2 mm篩備用。土壤顆粒組成用激光粒度分析儀(MS2000型，馬爾文，英國)測定，粒徑為0 mm

1.2 試驗裝置及試驗設(shè)計

試驗裝置主要由馬氏瓶、土箱和灌水器組成(圖1)。馬氏瓶用于提供恒定水頭，其內(nèi)徑為14 cm，高為100 cm，通過調(diào)整出口開關(guān)來控制出流流量；土箱材質(zhì)為有機(jī)玻璃，長×寬×高為50 cm×50 cm×50 cm；灌水器長度為20 cm，內(nèi)徑為4 cm，在其下部10 cm范圍內(nèi)均勻開孔，開孔率為25%，取1/4用玻璃膠粘于土箱一角，并用紗布包裹以避免土壤顆粒堵塞灌水器。試驗土壤初始含水率為10.20%(體積含水率)，設(shè)定土壤容重為1.35 g·cm-3，灌水歷時180 min。灌水器初始流量設(shè)2個梯度：2.0、3.0 L·h-1，其中2.0 L·h-1用于參數(shù)率定，3.0 L·h-1用于模型驗證。

1.馬氏瓶 Markow bottle; 2.橡膠管 Rubber tube; 3. 土箱 Test box;4. 濕潤體 Wet body; 5. 灌水器細(xì)部 Detail of emitter 圖1 試驗裝置示意圖及灌水器細(xì)部Fig.1 Schematic diagram of the test device anddetails of the irrigator

將試驗土壤按容重分層(每層5 cm)裝入土箱，層間打毛，自然沉降24 h后開始試驗，按先密后疏的原則設(shè)定時間間隔，觀測馬氏瓶讀數(shù)，同時以灌水器底部為起點，用直尺測量濕潤鋒在豎直向上、豎直向下和水平方向的運移距離。試驗結(jié)束時，在距離灌水器5、10、15 cm處使用土鉆取土，每隔10 cm 取樣1次，取土深度為50 cm。每組試驗重復(fù)3次，取其平均值。

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

2.1 土壤水分運動基本方程

涌泉根灌土壤水分入滲為三維柱狀水分入滲。假定土壤各向同性，則土壤水分為軸對稱分布，因此本文選用Hydrus-2D模型進(jìn)行模擬研究。本次室內(nèi)試驗時長為180 min，由于試驗時間較短，因此地表蒸發(fā)可忽略。不考慮空氣、溫度及滯后效應(yīng)對土壤水分的影響，則土壤水分運動用Richards方程描述為：

(1)

式中,θ為土壤體積含水率(cm3·cm-3)；t為入滲歷時(min)；K(h)為非飽和導(dǎo)水率(cm·min-1)；h為土壤水頭(cm)；x，z為坐標(biāo)軸，其中z軸向上為正(cm)。

2.2 定解條件

2.2.1 初始條件 土壤水分運動方程的初始條件為：

h(x,y,z,t)=h0(x,y,z)
0≤x≤X，0≤z≤Z，t=0

(2)

式中,h0(x,y,z)為土壤初始負(fù)壓水頭(cm)；X，Z為模擬計算區(qū)域橫向和豎向的最大長度，與室內(nèi)試驗土箱一致，即X=Z=50 cm。

圖2 數(shù)值模擬控制體Fig.2 Numerical simulation control body

2.2.2 邊界條件 圖1為涌泉根灌土壤水分運動數(shù)值模擬控制體示意圖，O點為坐標(biāo)系原點。由于室內(nèi)試驗不考慮蒸發(fā)，因此上邊界DE段：?h/?z=0 (2 cm≤x≤50 cm，z=50 cm，t>0 min)；右邊界為土箱邊壁，因此EF為零通量面，即?h/?z=0 (x=50 cm，0 cm≤z≤50 cm，t>0 min)；下邊界為土箱底部，因此FO也為零通量面，即?h/?x=0 (0 cm≤x≤50 cm，z=0 cm，t>0 min)；左邊界OA段為土箱邊壁，則有?h/?z=0 (x=0，0 cm≤z≤30 cm，t>0 min)；CD段為灌水器迷宮流道部分，因此?h/?z=0 (x=2 cm，40 cm≤z≤50 cm，t>0 min)。

AB和BC段為灌水器出流部分，由于涌泉根灌流量較大，灌水器出流量大于土壤入滲量，因此套管中形成積水，積水深度隨著時間不斷變化，開始階段逐漸增加，一定時間后達(dá)到穩(wěn)定，因此按變水頭邊界處理：h(x,y,z,t)=H(0 cm≤x≤2 cm，z=30 cm，t>0 min；x=2 cm，30 cm≤z≤40 cm，t>0 min)，其中H為灌水器內(nèi)積水深度，通過試驗過程中實測可得。

2.3 數(shù)值求解方法

利用Galerkin有限單元法完成對數(shù)學(xué)模型的數(shù)值計算，時間采用隱式差分法離散[13-14]。計算區(qū)內(nèi)網(wǎng)格劃分采用三角形單元，網(wǎng)格尺寸劃分見圖2，三角形邊長為3 cm，灌水器位置入滲界面處適當(dāng)加密網(wǎng)格，以提高計算精度。

2.4 模型參數(shù)率定

在Hydrus-2D模型中，土壤水分特征曲線采用Van Genuchten (VG)模型，其表達(dá)式為：

(3)

式中,θ(h)為土壤體積含水率，與土壤水吸力有關(guān)(cm3·cm-3)；θr為土壤滯留含水率(cm3·cm-3)；θs為土壤飽和含水率(cm3·cm-3)；h為土壤吸力(cm)，當(dāng)土壤處于飽和狀態(tài)時，h=0；α為土壤尺度函數(shù)，其值為土壤進(jìn)氣吸力值的倒數(shù)(cm-1)，容重較大的土壤，其對應(yīng)的進(jìn)氣吸力值較大，而土壤容重較小時，對應(yīng)的進(jìn)氣吸力值也較??；參數(shù)n與土壤水分特征曲線的坡度有關(guān)，對于砂質(zhì)土而言，其土壤特征曲線的坡度較緩，對應(yīng)的參數(shù)n就較大，對于粘性土壤而言，其土壤特征曲線的坡度較陡，其對應(yīng)的參數(shù)值n就較小；參數(shù)m=1-1/n(n>1)。

土壤非飽和導(dǎo)水率公式為：

K(θ)=Ks×S0.5[1-(1-S1/m)m]2

(4)

(5)

式中,S為飽和度；Ks為飽和導(dǎo)水率(cm·min-1)。θs、θr、Ks通過室內(nèi)試驗測定。利用RETC軟件獲取土壤水分特征曲線初始參數(shù)α、n，采用馬氏瓶流量為2.0 L·h-1時的實測數(shù)據(jù)進(jìn)行率定。

運行模型并在一定范圍內(nèi)調(diào)整參數(shù)，使模擬值盡可能接近實測值(表1)。利用SPSS 21.0軟件對濕潤鋒運移距離及不同時刻含水率的模擬值和實測值分別進(jìn)行對比，當(dāng)決定系數(shù)(R2)均大于0.05水平對應(yīng)的臨界決定系數(shù)，土壤濕潤鋒和含水率的均方根誤差分別小于1.5 cm和0.05 cm3·cm-3、F檢驗P值>0.05時，認(rèn)為精度滿足要求，確定模型參數(shù)。

3 模型準(zhǔn)確性評價

3.1 濕潤鋒運移距離模擬準(zhǔn)確性評價

圖3為涌泉根灌自由入滲實測和模擬對比。土壤各向濕潤鋒運移距離的模擬值和實測值隨入滲歷時變化見圖4。由圖4可知，在灌水結(jié)束時，土壤各向濕潤鋒運移距離的模擬值和實測值基本一致，均表現(xiàn)為隨著入滲歷時的增長先快速增加，而后緩慢增加的變化趨勢，這是因為在開始階段土壤含水率還較低，土壤空氣含量較高，因此水分運移速度也較快，隨著入滲歷時的逐漸增長，土壤空氣含量逐漸減小，因此土壤水分運移速度也逐漸變慢。

在灌水結(jié)束時，在豎直向上方向濕潤鋒均達(dá)到土箱表層；在豎直向下方向模擬值和實測值分別為22.2、21.1 cm，相差1.1 cm，相對誤差均為5.21%；在水平方向，模擬值和實測值分別為19.1、20.6 cm，相對誤差僅為-7.28%，可以看出，在灌水結(jié)束時，各向濕潤鋒的相對誤差均不超過±10%。

3.2 土壤剖面含水率模擬準(zhǔn)確性評價

圖5為灌水結(jié)束時，距離灌水器5、10、15 cm處土壤水分在0～50 cm土層內(nèi)的分布?？梢钥闯?，在灌水結(jié)束時，距離灌水器不同距離處土壤含水率剖面分布的模擬值和實測值基本一致，均表現(xiàn)為隨著土層深度的增加而先增大后減小，這是因為灌水器埋置于20 cm處，距離灌水器越近，土壤水分就越容易到達(dá)，因此相應(yīng)的土壤含水率就越大。

圖3 涌泉根灌自由入滲實測和模擬對比Fig.3 Comparison of measured and simulated of freeinfiltration of surge-root irrigation

表2 濕潤鋒運移距離模擬誤差分析

通過對模擬誤差的分析表明，在距離灌水器5 cm處，40～50 cm土層的模擬誤差最大，其值為9.92%；在距離灌水器10 cm處，20～30 cm土層的模擬誤差最大，其值為7.39%；在距離灌水器15 cm處，30～40 cm土層的模擬誤差最大，其值為-8.37%，可以看出在距離灌水器不同位置處，誤差的分布比較隨機(jī)，且均在±10%以內(nèi)。

圖5 距離灌水器不同距離處土壤含水率模擬值和實測值對比Fig.5 Comparison of simulated and measured values of soilwater content at different distances from the irrigator

表3 土壤含水率模擬誤差分析

4 結(jié) 論

根據(jù)土壤水動力學(xué)理論結(jié)合土壤水分的運動特征，建立了紅壤區(qū)涌泉根灌條件下土壤水分運動的數(shù)學(xué)模型，按照試驗條件設(shè)定數(shù)值模擬邊界條件，并采用Hydrus-2D軟件對所建模型進(jìn)行求解。通過涌泉根灌土壤水分隨時間的變化規(guī)律及空間分布規(guī)律的模擬值與實測值的對比可知，土壤濕潤鋒運移距離和土壤含水率隨時間變化的模擬值與實測值非常吻合，相對誤差均在±10%以內(nèi)，濕潤鋒運移距離和含水率模擬值和實測值之間RMSE分別小于1.300 cm和0.030 cm3·cm-3，且差異均未達(dá)到顯著水平，這說明Hydrus-2D軟件可以用于紅壤區(qū)涌泉根灌土壤水分運動的模擬。