基于HYDRUS的深施作業(yè)液肥水分入滲數(shù)值模擬

張澤龍，何建煒，張彬，張穎興，吳子彬，劉海軍，曹鵬，薛秀云

（華南農(nóng)業(yè)大學(xué)電子工程學(xué)院、人工智能學(xué)院，廣州510642）

0 引言

液肥憑借其施肥方式方便、節(jié)肥增產(chǎn)效果顯著等特點(diǎn)被廣泛應(yīng)用。液肥中所包含的水分和養(yǎng)分（鹽分）是作物生長的最基礎(chǔ)條件（李俊峰，2017）。土壤中水分跟養(yǎng)分的實(shí)際分布與作物根系最大化地吸收營養(yǎng)具有密切的聯(lián)系（Wang，2018）。因此充分了解不同條件下液肥的水分遷移規(guī)律，對提高作物肥料利用率的積極性，進(jìn)一步為制定合理的灌溉施肥制度具有重大意義。

數(shù)值模擬計(jì)算是研究液肥入滲特性的一個重要方法，在數(shù)值模擬中HYDRUS模型應(yīng)用最為廣泛（J.Si-munek等人，2012）、Simunek等人（2010）結(jié)合室內(nèi)試驗(yàn)和田間試驗(yàn)，驗(yàn)證了HYDRUS-2D模擬地下滴灌土壤水分分布具有較高可靠性。Bufon等人（2012）利用HYDRUS-2D進(jìn)行滴灌條件下的棉花含水量模擬，結(jié)果表明模擬值與測量值數(shù)據(jù)擬合程度較好。HYDRUS模型可對不同條件下水鹽運(yùn)移規(guī)律進(jìn)行模擬計(jì)算，為區(qū)域田間水肥管理提供理論依據(jù)和參考價(jià)值。

當(dāng)前我國大部分地區(qū)仍采用施肥效率低下、化肥損耗嚴(yán)重的施肥機(jī)械或非機(jī)械化施肥方式。為此，提出了一種能夠定量精準(zhǔn)均勻施肥的液肥變量深施控制系統(tǒng)。為給液肥深施系統(tǒng)精準(zhǔn)深施效果提供科學(xué)依據(jù)，進(jìn)一步提高系統(tǒng)深施效率，對液肥深施系統(tǒng)深施作業(yè)條件下的注肥壓力、注肥深度、土壤質(zhì)地3個因素進(jìn)行HYDRUS建模，探究各因素對水分入滲的影響。

1 試驗(yàn)方法

試驗(yàn)采用單因素試驗(yàn)法進(jìn)行研究，不考慮各個因素間的交互作用。試驗(yàn)中包括了三個主要影響因素（A系統(tǒng)注肥壓力、B系統(tǒng)注肥深度、C土壤質(zhì)地類型），其中分別對系統(tǒng)壓力和注肥深度設(shè)置了3個水平，土壤質(zhì)地設(shè)置了2個水平。為更準(zhǔn)確地觀測土壤濕潤體特征值變化規(guī)律，系統(tǒng)壓力控制水平為：0.03MPa、0.04MPa、0.08MPa。由于成年的南方柑橘樹的根系主要集中分布在地表下15～55cm處，故系統(tǒng)注肥深度水平控制為：25cm、35cm、45cm（薛秀云等人，2020）；通過查閱中國科學(xué)引文數(shù)據(jù)庫相關(guān)文獻(xiàn)，了解華南地區(qū)主要土壤質(zhì)地類型，壤質(zhì)地控制水平為砂壤土、壤土。在數(shù)值模擬計(jì)算試驗(yàn)中，對以上3因素組合，注肥時(shí)間均設(shè)定為10min，并在注肥開始后記錄垂直剖面下不同時(shí)刻9個方向土壤濕潤鋒運(yùn)移規(guī)律、濕潤體形態(tài)大小，如圖1所示。對模擬試驗(yàn)得到的仿真數(shù)據(jù)利用Origin-Pro 2018軟件進(jìn)行數(shù)據(jù)處理及圖表繪制，以此來分析不同因素條件下濕潤體時(shí)空變化規(guī)律。

圖1 垂直觀測面下濕潤鋒記錄示意圖

2 數(shù)值模擬模型

選擇HYDRUS（2D/3D）對液肥深施系統(tǒng)深施作業(yè)下土壤水分入滲規(guī)律進(jìn)行動態(tài)模擬。HYDRUS模型是由美國鹽土改良實(shí)驗(yàn)室（US Salinity Laboratory）研發(fā)的一款用于模擬飽和-非飽和多孔介質(zhì)中水分、能量及溶質(zhì)運(yùn)移的數(shù)值模型。該模型適用于不同初始條件與各種恒定或可變邊界條件，能較好地模擬土壤中水分、熱以及溶質(zhì)的時(shí)空變化、空間分布及運(yùn)移規(guī)律（遲卉等人，2014）。

2.1 土壤水分運(yùn)動方程

本研究所探討的基于液肥變量深施系統(tǒng)的注肥深施入滲為恒定流量條件下的三維柱狀面源水分入滲。假定土壤為均質(zhì)，固相骨架不變形的剛性多孔介質(zhì)，忽略土壤水分運(yùn)動中的滯后效應(yīng)，同時(shí)不考慮溫度及空氣阻力對土壤水分運(yùn)動的影響，土壤水分運(yùn)動可用Richards方程來描述：

式中：θ為土壤體積含水率，cm3/cm3；ψ為土壤負(fù)壓水頭，cm；K（θ）為非飽和土壤導(dǎo)水率，cm/min；t為時(shí)間，min；x、y為水平面坐標(biāo)，cm；z為垂直面坐標(biāo)，cm。

2.2 模擬區(qū)域

將試驗(yàn)土壤區(qū)域的幾何模型規(guī)格設(shè)置為80cm×80cm×100cm（長×寬×高），注肥槍由直徑2cm的圓柱槍體和圓錐形槍頭組成，系統(tǒng)注肥處設(shè)置為模型計(jì)算區(qū)域中央，模型入滲面是三維柱狀出流邊界。采用垂直剖面法觀測豎直平面內(nèi)濕潤體形態(tài)與特征值，efgh垂直剖面設(shè)定為YZ平面，垂直觀測面設(shè)為距注肥槍口0.1cm處。模型示意圖如圖2。

圖2 模型示意圖

2.3 初始條件

本研究模型計(jì)算的初始條件為初始條件為初始含水率，因此求解土壤水分運(yùn)動的初始條件為：

式中：H0（x,y,z）為土壤初始負(fù)壓，假設(shè)計(jì)算區(qū)域內(nèi)土壤初始負(fù)壓處處相等，設(shè)置為-100000cm（對應(yīng)不同土壤質(zhì)地初始含水率為：砂壤土—0.06cm3/cm3；壤土—0.09cm3/cm3）土壤初始含水量由土壤水分特性曲線轉(zhuǎn)換得到。

2.4 邊界條件

如圖3所示，本研究模型的上邊界A’B’C’D’為土壤表層，為大氣邊界條件，垂直通量為零；前后左右四面邊界AA’B’B、CC’D’D、AA’D’D、BBC’C為自由排水邊界。底面邊界ABCD，其垂直通量為零。而對于注肥槍末端O及槍管圓柱側(cè)面OE為出流邊界，是一個比較復(fù)雜的運(yùn)動邊界。試驗(yàn)設(shè)置的系統(tǒng)注肥壓力較大，注肥槍出口流量遠(yuǎn)大于土壤入滲速率。因此要將系統(tǒng)注肥壓力轉(zhuǎn)化為模型能夠設(shè)置的邊界條件，即為可變壓力水頭邊界處理。對于注肥槍入滲面的可變壓力水頭按因素組合處理去不同數(shù)值還實(shí)現(xiàn)。通過系統(tǒng)壓力與水柱高度轉(zhuǎn)換公式將系統(tǒng)壓力轉(zhuǎn)化為水柱高度：已知一個大氣壓A0=0.1MPa等于h0=10m高水柱，故利用公式（3）將系統(tǒng)壓力轉(zhuǎn)化為壓力水頭。計(jì)算得出壓力水頭見表1。

表1 三種系統(tǒng)壓力對應(yīng)壓力水頭

式中：A0為標(biāo)準(zhǔn)大氣壓，0.1MPa；h0為一個大氣壓下水柱高度，1000cm；A為系統(tǒng)壓力，MPa；h為壓力水頭，cm。

2.5 數(shù)值求解方法

HYDRUS-3D流體建模采用Galerkin有限單元分析法求解三維土壤水分運(yùn)移問題，時(shí)間采用隱式差分法，迭代參數(shù)和殘差參數(shù)均按默認(rèn)值設(shè)置。模型模擬土壤水流區(qū)域選用Triangular Prism（四面體）網(wǎng)格進(jìn)行劃分。模型邊界離散化采用軟件自動生成的邊界密度節(jié)點(diǎn)。

3 結(jié)果與分析

3.1 系統(tǒng)壓力對水分入滲的影響

圖3為注肥深度為35cm，土壤質(zhì)地為砂壤土條件下，不同系統(tǒng)壓力下濕潤鋒動態(tài)運(yùn)移過程。由圖3可知，系統(tǒng)壓力對土壤濕潤體形狀大小的影響較為明顯，濕潤鋒垂直向上擴(kuò)散距離、濕潤鋒水平擴(kuò)散距離、濕潤鋒垂直向下擴(kuò)散距離均隨注肥時(shí)間的增加而增加，但各方向濕潤鋒運(yùn)移速率在注肥后期明顯逐漸減小。在整個注肥過程中，由于液肥深施系統(tǒng)施肥為大壓力施肥，供水強(qiáng)度大于土壤入滲能力，又因在豎直方向上水分受重力勢的影響，導(dǎo)致濕潤鋒垂直向下運(yùn)移速率略大于水平運(yùn)移速率。而在豎直向上方向土壤孔隙相對較大，水分在運(yùn)移過程中受到阻力相對較小，且在土壤水分入滲過程中土壤基質(zhì)勢起到重要作用，所以垂直向上濕潤鋒運(yùn)移速率略大于垂直向下濕潤鋒。至注肥結(jié)束，三個方向的濕潤鋒運(yùn)移距離差異很小，所以深施作業(yè)下濕潤體形狀近似“球體”，隨著系統(tǒng)壓力的增大，濕潤體大小也隨著增大。由于液肥深施系統(tǒng)系統(tǒng)壓力決定了系統(tǒng)輸出流量，當(dāng)系統(tǒng)壓力越大時(shí)，系統(tǒng)輸出流量越大，濕潤體形狀因而隨著增大。

圖3 不同系統(tǒng)壓力條件下濕潤鋒運(yùn)移動態(tài)變化圖

對在垂直向上、水平擴(kuò)散距離和垂直向下入滲深度三個方向的濕潤鋒運(yùn)移進(jìn)行擬合，如圖4所示。從中可以發(fā)現(xiàn)不同系統(tǒng)壓力條件下垂直剖面下的垂直向上濕潤鋒運(yùn)移距離X、水平濕潤鋒運(yùn)移距離R和垂直向下濕潤鋒運(yùn)移距離H與注肥時(shí)間t均存在顯著的冪函數(shù)關(guān)系，R2均大于0.98，表明擬合結(jié)果與模擬值吻合度較高。其中，對于同一濕潤鋒運(yùn)移方向上冪函數(shù)指數(shù)差異性較小，其平均值分別為（a）0.2487、（b）0.33、（c）0.3213。而其冪函數(shù)系數(shù)均隨系統(tǒng)壓力的增大而增大，這表明各方向濕潤鋒運(yùn)移距離受系統(tǒng)注肥壓力影響，均呈現(xiàn)正相關(guān)。對于同一系統(tǒng)壓力條件下不同方向濕潤鋒冪函數(shù)系數(shù)均呈現(xiàn)相同規(guī)律，具體表現(xiàn)為：垂直向上濕潤鋒的冪函數(shù)系數(shù)>水平濕潤鋒的冪函數(shù)系數(shù)>垂直向下濕潤鋒的冪函數(shù)系數(shù)。

圖4 不同系統(tǒng)壓條件力下濕潤鋒隨時(shí)間運(yùn)移距離

3.2 注肥深度對水分入滲的影響

圖5為注肥深度為35cm，土壤質(zhì)地為砂壤土條件下時(shí)，不同注肥深度濕潤鋒動態(tài)運(yùn)移過程。從圖5可知，注肥深度對濕潤體形態(tài)和大小影響較小，不同處理下深施作下土壤濕潤體垂直剖面均近似于圓形且形態(tài)大小差異很小，濕潤鋒水平方向最大運(yùn)移距離略高于注肥槍注肥槍口所在深度。在整個注肥過程中，各處理垂直向上、水平方向、垂直向下濕潤鋒運(yùn)移距離接近，至注肥結(jié)束，各方向濕潤鋒運(yùn)移距離差距始終不大，所以導(dǎo)致了深施模式下土壤濕潤體形狀最終近似為圓形。雖然不同注肥深度處理對土壤濕潤體形態(tài)大小影響較小，但對土壤濕潤體的空間分布影響很大，隨著注肥深度的增加，土壤濕潤體分布所處深度越大。因此，針對不同作物種類制定合理的灌溉制度時(shí)，要根據(jù)作物種類根部生長狀況來確定合適的注肥深度，以達(dá)到在有效防止地表蒸發(fā)和深層滲漏前提下滿足作物水肥需求。

圖5 不同注肥深度條件下濕潤鋒運(yùn)移動態(tài)變化圖

土壤濕潤體的垂直方向（向上和向下）濕潤距離和最大水平濕潤半徑是衡量液肥入滲效果的關(guān)鍵指標(biāo)。對注肥深度分別為25cm、35cm和45cm三種梯度條件下各方向濕潤鋒隨時(shí)間的變化過程進(jìn)行擬合。從表2中可知，在相同注肥時(shí)間內(nèi)，土壤濕潤鋒垂直向上、水平和垂直向下運(yùn)移距離隨注肥深度的增加不發(fā)生明顯變化，說明注肥深度對濕潤鋒運(yùn)移影響較小。

3.3 注肥深度對水分入滲的影響

圖6為施肥深度為35cm，系統(tǒng)壓力為0.03MPa條件下的不同土壤質(zhì)地濕潤鋒運(yùn)移動態(tài)變化圖。由圖6可知，對于不同土壤質(zhì)地條件下的液肥深施作業(yè)而言，土壤濕潤體在注肥結(jié)束時(shí)均近似為“圓形”，濕潤體形狀差異不大；而濕潤鋒運(yùn)移距離卻差異較大，各方向濕潤鋒運(yùn)移距離在砂質(zhì)壤土條件下均大于在壤土條件下，由此，導(dǎo)致了土壤濕潤鋒形狀在大小上差異較大。這主要是因?yàn)橥寥浪秩霛B主要受土壤孔隙度的影響，土壤孔隙度越多，氣相比例越多，土壤內(nèi)部連通性越好，越有利于土壤水分入滲的進(jìn)行，刑杰等人（2013）研究表明砂質(zhì)壤土黏粒含量約為10.43%，壤土黏粒含量約為18.63%，孔隙度隨土壤黏粒含量的增加而減小。因此，土壤濕潤體體積在砂質(zhì)壤土條件下大于在壤土條件下。

圖6 不同土壤質(zhì)地條件下濕潤鋒運(yùn)移動態(tài)變化圖

表2 不同注肥深度各方向濕潤鋒擬合結(jié)果

對濕潤鋒隨時(shí)間的變化過程進(jìn)行處理與分析，發(fā)現(xiàn)垂直向上濕潤鋒運(yùn)移距離X，水平方向濕潤鋒運(yùn)移距離R和垂直向下濕潤鋒運(yùn)移距離與注肥時(shí)間t均呈較好的冪函數(shù)關(guān)系，即：

式中：a、c、e分別為垂直向下、水平、垂直向上的入滲系數(shù)，b、d、f分別為下、水平、垂直向上的入滲指數(shù)，數(shù)據(jù)擬合結(jié)果見表3。

表3 各方向濕潤鋒運(yùn)移距離與注肥時(shí)間擬合結(jié)果

從表3中可以看出，無論垂直向下，水平方向還是垂直向上濕潤鋒的入滲系數(shù)均受土壤質(zhì)地類型影響較大，隨土壤黏粒含量的增加而下降。而入滲系數(shù)的物理意義表征的是土壤入滲開始后第一單位時(shí)間內(nèi)濕潤鋒入滲的能力，這也說明了在不同土壤質(zhì)地在第一單位時(shí)間段末濕潤鋒運(yùn)移距離差異較大。

4 結(jié)語

利用HYDRSUS-3D模型建立了系統(tǒng)深施作業(yè)下土壤水分運(yùn)動模型，對濕潤體特征值進(jìn)行觀測，分析了系統(tǒng)壓力、注肥深度、土壤質(zhì)地對土壤水分入滲的影響，得出以下結(jié)論：

（1）系統(tǒng)壓力對濕潤體大小，均隨系統(tǒng)壓力的增大而增大；注肥深度增大對濕潤體形態(tài)大小無顯著影響，對濕潤體空間分布位置有著顯著影響，隨著注肥深度的下移，土壤濕潤體空間分布也隨之下移；不同土壤質(zhì)地條件對濕潤體形態(tài)影響較小，濕潤體形態(tài)均近似為圓形，而對對濕潤體大小移影響較大，土壤質(zhì)地越黏重，濕潤體越小。

（2）深施作業(yè)下土壤濕潤體在形狀上差異較小，濕潤鋒運(yùn)移距離隨系統(tǒng)壓力的增大而增大，隨土壤黏重程度的增大而減小，而與注肥深度無顯著關(guān)系。不同方向土壤濕潤鋒運(yùn)移距離與注肥時(shí)間都呈顯著冪函數(shù)關(guān)系，且在相同注肥時(shí)間內(nèi)均符合垂直向上>垂直向下>水平方向的規(guī)律。

通過以上研究，可為液肥深施系統(tǒng)的深施作業(yè)提供理論依據(jù)，但未考慮實(shí)際情況下土壤水力特性的差異性，因而仍需通過室內(nèi)土箱試驗(yàn)進(jìn)一步研究系統(tǒng)壓力、注肥深度、土壤質(zhì)地對土壤水分入滲的影響。