姚毓香,高喜杰,高鵬洋,朱瑞祥,黃玉祥
(西北農林科技大學 機械與電子工程學院,陜西 楊凌 712100)
適宜的土壤水分入滲能力對于保持農田良性水文循環尤為重要[1]。深松耕作為消除土壤板結、打破犁底層、提高土壤水分入滲能力的有效措施,在農業生產中得到了普遍應用[2-3]。
深松鏟橫向間距(鏟距)是影響土壤擾動和深松耕效果的主要因素之一,對土壤水分入滲特性具有重要影響[4-5]。Franzluebbers等研究表明:相較于免耕,深松后土壤容重降低了12%,土壤水分入滲速率提高27%[6]。Hamza等分析了不同鏟距和耕深組合對土壤容重、硬度、土壤含水量及土壤水分入滲速率的影響[7]。深松有利于土壤大孔隙形成,而土壤大孔隙引起的優先流能夠有效提高水分入滲能力,是影響土壤水分入滲速率的主要原因[8]。許迪等研究發現,在10~30cm土層內,深松后直徑大于50μm的土壤孔隙體積明顯增多而直徑小于10μm的孔隙體積卻顯著減少[9]。Lipiec等的研究進一步發現,孔隙分布和水分入滲之間存在很強的對數關系,3h內的累積入滲量與土壤表層0~10cm的孔隙度密切相關[10]。目前,國內外學者主要開展了深松對土壤物理性質及水分入滲特性的影響研究,而對深松鏟鏟距與土壤水分入滲關系的研究較少。
本研究以箭型深松鏟為對象,通過離散元仿真提取不同鏟距條件下的土壤剖面輪廓,將其導入Hydrus2D軟件進行數值模擬,并結合田間試驗,研究不同鏟距(30、35、40、45、50cm)對土壤垂直剖面內水分入滲特性的影響,為深松鏟橫向間距的選擇提供參考,有助于提高深松土壤水分的利用效率和效益。
試驗區位于陜西省咸陽市楊凌區(108°03'E,34°18'N),試驗區土壤為壤質塿土,粘性較大。常年采用旋耕機進行淺層作業,地勢平坦。試驗前測定0~15cm、15~30cm耕層內的土壤緊實度分別為0~1 800kPa和1 800~6 000kPa,土壤體積含水量分別為17.3%、19.2%。
深松過程中作業機具選用約翰迪爾904拖拉機(長×寬×高=4 320mm×2 150mm×2 880mm),機器作業速度為3km/h。試驗區耕作層厚度介于18~20cm之間,為保證完全打破犁底層,設置耕深為30cm,每組試驗的小區面積為60m×3m,3次重復。深松后隨機選取試驗區土壤坑形、壟形擾動輪廓各3組進行測繪,并計算土壤擾動面積,如圖1所示。

圖1 土壤擾動輪廓測量
深松后在土壤垂直剖面內的壟形中心線位置(位置A)和壟間中心線位置(位置B)埋設Trime測管,并測量Trime測管0~10cm、10~20cm、20~30cm、30~40cm、40~60cm、60~100cm等6個土層深度的土壤體積含水量,如圖2所示。試驗時,入滲區域設置在A位置,通過內徑32.5cm 和外徑62.5cm 的2個金屬環控制,將兩環以Trime 測管為圓心埋入土中約10cm深度處,然后內外環同時注水,水深控制在5cm左右,每次試驗持續時間為120min。根據試驗耗水量,使用Trime-IPH手持式TDR測量A、B 位置各土層深度的土壤體積含水量[11]。

圖2 深松后土壤體積含水量觀測點布設位置

圖3 離散元土壤耕層結構模型
1.3.1 深松鏟和土壤顆粒邊界模型構建
本研究以國標(JBFF 9788.1999《深松鏟和深松鏟柄》)中所規定的箭型深松鏟鏟尖和圓弧形深松鏟鏟柄為對象,利用SoildWorks2013三維軟件進行圖形繪制,構建深松鏟邊界模型。
相關研究表明,實際土壤顆粒形狀是不規則的,因此在離散元仿真中不能直接用球形顆粒統一替代[12]。根據現有研究并結合楊凌地區的土壤結構,以離散元中球形單元顆粒為基礎,對土壤顆粒進行分類建模,主要包括球形塊狀顆粒、核狀顆粒及柱狀顆粒[13]。
1.3.2 仿真參數設置
離散元仿真過程中,模型參數主要包括材料參數及接觸參數兩大類。其中,材料參數主要包括深松鏟和土壤顆粒的泊松比、密度及剪切模量等。本文中土壤密度通過實際測量獲得,為2 600kg/cm3。其余參數主要根據現有文獻獲取,設置土壤泊松比為0.3,剪切模量為1.0×106G/Pa、深松鏟的泊松比為0.3,密度為7 865kg/cm3、剪切強度為8.19×1010G/Pa[14]。
接觸參數主要包括土壤-土壤、土壤-深松鏟間的靜摩擦因數、動摩擦因數以及恢復系數,同樣可以通過現有文獻[15-16]獲取。其中,土壤-土壤靜摩擦因數、動摩擦因數及恢復系數分別為0.33、0.17、0.6;土壤-深松鏟的接觸參數則分別為0.46、0.11、0.6。
1.3.3 離散元仿真模型
在離散元模型中,土壤與土壤顆粒之間存在粘結鍵。根據現有研究,設置土壤顆粒之間的接觸模型為Hertz-mindlin Bonding模型[17]。在進行建模時,首先建立長、寬、高分別為1 500mm×1 500mm×500mm的土槽模型;其次,設置深松鏟耕深為30cm,前進速度為3km/h(0.83m/s),時間步長為20%,仿真時間為10s。為模擬顆粒分布的實際情況,設置土壤顆粒大小呈正態分布。在土槽模型中依次生成心土層(柱狀顆粒)、犁底層(核狀顆粒)、耕作層(球形塊狀和球形單元顆粒)的土壤顆粒,各層的土壤顆粒數分別為80 000、75 000、40 000、40 000。
1.4.1 模型建立
Hydrus2D是一個可用來模擬土壤水分運動、溶質運移、熱量傳遞的有限元計算機模型。該模型的水流控制方程采用Richards方程,利用Galerkin有限元法對土壤剖面進行求解。將深松后0~100cm的土壤垂直剖面劃分為5層,通過離散元仿真提取土壤擾動輪廓,并將其導入Hydrus2D模型中。土壤擾動輪廓的主要參數如表1所示。設置土壤水分運動模擬區域,如圖4所示。利用Hydrus2D模型對不同鏟距條件下的土壤水分入滲過程進行數值求解。

表1 土壤擾動輪廓分析

Sp.鏟距(cm) d.耕深(cm) X、S.土壤體積含水量觀測點
1.4.2 模型求解的初始和邊界條件
1)初始條件:在Hydrus2D模型中,土壤垂直剖面內0~100cm的土壤初始體積含水量按照田間實測數據進行分層設置,并假定各層的土壤初始體積含水量在水平方向均勻分布[18]。
2)邊界條件:模擬試驗中,在壟形中心線位置設置恒定壓力水頭為5cm,上邊界其余位置按大氣邊界條件設置,左右邊界設為零通量邊界,下邊界無水量交換,模擬時不影響土壤水分入滲過程,設為自由排水邊界。
1.4.3 模型參數設置
深松后,以地表為基準面,用環刀取擾動區域和未擾動區域0~10cm、10~20cm、20~40cm、40~60cm和60~100cm土層的土樣,每層重復取3個,測定土壤容重。通過MS2000型激光粒度儀測定土壤顆粒組成,CR21G高速恒溫冷凍離心機測定土壤水分特征曲線。利用Hydrus2D模型中的Rosetta模塊預測土壤水力特性參數,如表2所示。

表2 試驗土壤的水力特性參數
1.4.4 模型網格劃分
在Hydrus2D模型中通過FE-Mesh模塊對求解區域進行網格劃分,如圖5所示。為提高模擬精度,設置擾動區域三角形網格間隔為1cm,離散為15 640個網格;未擾動區域的三角形網格間隔為5cm,離散為7 952個網格。設定模擬時間為120min,初始時間步長為0.01min,最小時間步長為0.001min,最大時間步長為5min。土壤含水量允許偏差為0.001,壓力水頭允許偏差為0.1cm,最大迭代次數為10次。

圖5 模型區域網格劃分
1.4.5 模型驗證
為求解模型并驗證所建模型的合理性,采用 Hydrus2D模型對不同鏟距條件下的土壤水分入滲過程進行數值模擬。將模擬得到的土壤體積含水量與試驗結果進行對比分析,采用均方根誤差(RMSE)和決定系數(R2)檢驗土壤水分運動模型的合理性[19]。
由表3可知:土壤擾動面積的仿真值均小于實際值,土壤未擾動面積的仿真值均大于試驗值。其中,鏟距為30cm時,土壤未擾動面積的相對誤差最大(13.2%);鏟距為50cm時土壤擾動面積的相對誤差最小(0.6%),均滿足相關試驗的誤差要求[14]。試驗結果表明:基于EDEM法提取深松土壤擾動輪廓具有可行性。

表3 EDEM仿真與田間試驗結果對比
入滲結束后,對比分析土壤垂直剖面內的體積含水量分布狀況。如表4所示,土壤體積含水量模擬值和實測值的RMSE值在0.025~0.059cm3/cm3之間,R2均大于0.826。模擬值和實測值之間的對應關系良好,表明可以通過Hydrus2D模型來模擬不同鏟距條件下的土壤水分入滲過程。

表4 土壤體積含水量模擬值和實測值對比
由圖6可知:入滲過程結束后,在垂直方向上不同鏟距條件下的土壤含水量變化一致,整體呈現先增加、后減小的趨勢;隨著土壤水分入滲過程的進行,各層土壤水分均得到補給,但補給量隨入滲深度的增加逐漸減少。在不同鏟距條件下,土壤垂直剖面內0~10cm、10~20cm、20~30cm、30~40cm、40~60cm、60~100cm土層的土壤含水量分別增加了141%、155%、149%、78.1%、5.5%、4.3%。試驗表明:在土壤水分入滲過程中,土壤水分優先補給土壤表層,深層土壤的含水量變化不大。

圖6 不同鏟距下的土壤水分分布
鏟距對土壤水分分布有重要影響,隨著鏟距的增加,土壤水分富集區逐漸向淺層移動。當鏟距為50cm時,土壤水分富集區位置最淺,位于13~21cm土層內;當鏟距為30cm時,土壤水分富集區位置最深,位于32~38cm土層內,且鏟距每增加10cm,土壤水分富集區位置平均向淺層移動約9cm,通過調整鏟距可以改變土壤水分富集區位置。因此,在農業生產中,對淺根系作物宜選用大鏟距作業,對深根系作物,宜選用小鏟距作業。
入滲速率是表征土壤水分入滲能力的指標之一,其快慢可以反映土壤接收外界降水或灌溉水速度的能力。如圖7所示:在整個入滲過程中,不同鏟距條件下的土壤水分入滲速率均呈現先減小、后穩定的趨勢。初始入滲時,隨著鏟距的增加,土壤水分入滲速率先增大、后減小;當鏟距為30、35、40、45、50cm時,土壤水分初始入滲速率依次為0.532、0.554、0.656、0.636、0.616cm/min,且隨著入滲時間的增加,土壤水分入滲速率逐漸減小并趨于穩定。

圖7 土壤水分入滲速率
隨著鏟距的增加,土壤水分穩定入滲速率逐漸增大。當鏟距分別為30、35、40、45、50cm時,土壤水分穩定入滲速率依次為:0.084、0.085、0.088、0.090、0.093cm/min。通過Spss21.0對不同鏟距條件下的土壤擾動面積和土壤水分穩定入滲速率進行相關性分析,結果表明:土壤水分穩定入滲速率與土壤擾動面積之間的相關性系數為0.996,二者之間存在極顯著的正相關關系(P<0.01),土壤水分穩定入滲速率隨著土壤擾動面積的增加逐漸增大。這是由于隨著鏟距的增加,被推移到地表的土壤沿深松鏟前進方向向兩側抬升的范圍增大,使土壤擾動面積增加,土壤破碎范圍隨之增大,進而影響土壤水分穩定入滲速率[20]。
濕潤鋒垂直運移距離反映了土壤水分的下滲能力。由圖8可知:隨著入滲時間的增加,濕潤鋒垂直運移距離的增長速率呈現先減小、后穩定的趨勢。初始入滲時,濕潤鋒垂直運移距離的增長速率為2.41cm/min,入滲至20min時逐漸穩定,濕潤鋒垂直運移距離以0.69cm/min的速率穩定增加。

圖8 濕潤鋒垂直運移距離
隨著鏟距的增加,濕潤鋒垂直運移距離略有波動,整體呈減小趨勢。當鏟距分別為30、35、40、45、50cm時,濕潤鋒垂直運移距離依次為:41.12、41.66、39.11、39.34、35.30cm。結合鏟距對土壤的擾動情況可知,當鏟距為30、35、40、45、50cm時,土壤未擾動面積分別為100.00、156.25、225.00、306.25、400.00cm2,濕潤鋒垂直運移距離與雙鏟間未擾動面積的相關性系數為-0.914,表明二者之間存在顯著的負相關關系(P=0.03)。這是由于在深松過程中,相鄰兩鏟之間存在協同作用力,位于雙鏟間的土壤會從不同方向被推擠到同一位置,進入到同一位置的土壤相互擠壓,使土壤發生破碎[20]。另外,隨著鏟距的增加,受雙鏟間交互作用影響被推擠到同一位置的土壤減少,導致土壤未擾動面積增加,土壤破碎程度降低,從而降低了土壤水分入滲深度。
1)離散元仿真與試驗獲取的土壤剖面輪廓形狀基本吻合。Hydrus2D模型的土壤體積含水量模擬值與實測值的R2均大于0.826,擬合效果較好。基于EDEM-Hydrus2D相結合的方法,可以用于研究不同鏟距條件下的土壤水分入滲過程。
2)鏟距對土壤水分富集區產生重要影響,隨著鏟距的增加,土壤水分的富集區位置向上移動。在農業生產中,可根據作物需水特性選擇合適的鏟距,以提高深松土壤水分的利用效率和效益。
3)土壤水分穩定入滲速率隨鏟距的增加而增大,濕潤鋒垂直運移距離隨鏟距的增加略有波動,整體呈減小趨勢。在滿足作物對土壤水分需求的前提下,適當增加鏟距,有利于提高土壤水分的入滲速率,減少耕層水分的滲漏損失。