精準滴灌農田排水過程數值模擬

麥榮幸， 康躍虎， 裴 亮， 李曉彬

(1.中國科學院 地理科學與資源研究所 陸地水循環及地表過程重點實驗室, 北京 100101； 2.中國科學院大學, 北京 100049)

1 研究背景

隨著水資源的日益緊缺和現代科學技術的發展，精準灌溉是未來農業灌溉的發展趨勢。滴灌作為現代灌溉技術，可根據作物需水需肥規律進行適時適量施肥灌溉，可以最大限度地減少蒸發、深層滲漏和地表徑流等帶來的水分和養分的損失，是實現精準灌溉的主要技術之一[1-2]。Wan等[3]、Kang等[4]和Wang等[5]在對滴灌農田水分養分循環過程、作物耗水規律、作物生長與耗水關系等進行了長期、系統、深入研究的基礎上，提出了監測滴頭正下方20 cm深度土壤基質勢下限的滴灌農田水鹽調控方法，為精準滴灌的實現提供了簡單、可靠的方法。

根據土壤基質勢精準滴灌農田條件下具有少量高頻灌溉的特點，土壤水分得到高度控制，農田土壤整體處于水分非飽和狀態，并且在時間和空間上均保持一個較為穩定的關系[6-8]。只在關鍵生育期灌溉的傳統地面灌溉農田條件下，每次灌溉時土壤表層為飽和狀態，長時間的灌溉間隔期內受降雨、作物耗水、蒸發等影響，土壤水分在時間和空間上表現出不穩定、變化劇烈的關系[9]。因此，與傳統地面灌溉農田相比，精準滴灌農田的土壤水分運動與再分布過程、灌溉或降雨產生的深層滲漏量、灌溉產生的棄水、降雨產生的地表徑流等均不相同，造成滴灌農田的排水過程和排水量也與傳統地面灌溉農田有所不同[10]。隨著滴灌在灌區農業中的大面積應用，以及灌區管理的現代化、信息化發展，亟需對精準滴灌農田的排水過程進行研究，為農田排水系統的規劃設計與改造、灌區灌排信息化管理提供參考。

滴灌農田排水過程分為土壤水分入滲至淺層地下水和隨地下水運動流入排水系統兩個過程。目前國內外學者對滴灌條件土壤水分在非飽和帶的運動、區域淺層地下水運動以及土壤水與地下水運動耦合方面做了大量的研究[11-13]，構建起了單點源和多點源的土壤水分在非飽和帶的運動模型，即LINKFLOW、MODFLOW、FEFLOW等地下水運動模型以及SWMS-3D等耦合模型[14-17]。綜觀現有研究成果，對于精準滴灌農田條件下土壤水分隨灌溉水、降水入滲至淺層地下水再運動到排水系統的整個排水過程的研究還未見到報道。本文在監測土壤基質勢指導滴灌灌溉條件下，根據農田土壤水分分布特征，綜合考慮淺層地下水運動的尺度問題，通過構建飽和帶-非飽和帶水分運動耦合模型，建立精準滴灌農田排水過程系統模型，對精準滴灌農田排水過程進行準確高效地數值模擬，對現代灌區排水系統規劃及信息化管理具有重要的應用價值和意義。

2 精準滴灌農田排水過程的數學模型

隨著灌溉和降雨，土壤水分垂向運動入滲至淺層地下水中，在水力梯度作用下隨淺層地下水流入到排水溝中。理論上，農田中離排水溝越遠的區域，地下水埋深越淺。這就造成了在排水過程中，非飽和帶厚度在空間上的變化，不同位置的深層滲漏水量也不一樣，排水過程也不相同。在非飽和帶中，土壤水分除了存在垂向交換之外，還存在水平交換，水分在非飽和帶的運動可采用三維Richards方程。在農田中，飽和帶含水層厚度較大，淺層地下水向排水溝的運動主要以水平運動為主。水分在飽和帶的運動可采用基于三維Richards方程沿飽和含水層垂向積分平均得到的水平二維方程[18]。在農田排水過程中，非飽和帶通過深層滲漏水量補給飽和帶。通過深層滲漏水量將非飽和帶水分運動方程與飽和帶水分運動方程耦合起來，構建精準滴灌農田排水過程系統模型。

2.1 模型的控制方程

2.1.1 非飽和帶水分運動控制方程 在滴灌條件下，當滴灌帶和滴頭間距較小時，滴灌帶間濕潤鋒會發生重疊，形成濕潤帶，土壤水分運動為三維運動。假設土壤為均質、各向同性的剛性多孔介質，不考慮溫度、生物或化學反應對水分運動的影響，也不考慮土壤水分運動的滯后效應，土壤非飽和帶水分運動可用三維Richards方程描述:

(1)

式中：K(h)為非飽和導水率，cm/min；θ為土壤體積含水率，cm3/cm3；h為土壤負壓水頭，cm；t為時間，min；x、y、z為空間坐標，x、y為水平方向坐標，z為垂直方向坐標，cm；規定z向下方向為正。

2.1.2 飽和帶水分運動控制方程 在飽和帶中，地下水流運動為三維運動，地下水流運動可用以 為因變量的Richards方程描述:

(2)

式中：Ks為土壤飽和導水率，cm/min；μ為比儲水系數;I為源匯項。

在農田排水過程中，飽和帶含水層的厚度較大，地下水流運動主要以水平方向運動為主，垂向運動較小。假設飽和帶含水層隔水底板為水平，將公式(2)方程沿著飽和帶含水層厚度垂向積分即可平均得到水平二維方程。

公式(2)左邊的積分項為：

(3)

(4)

公式(2)右邊的積分項為：

(5)

綜合以上各式，飽和帶地下水運動方程可采用以下水平二維方程形式：

(6)

2.1.3 飽和帶和非飽和帶的耦合項 在農田排水過程中，非飽和帶通過深層滲漏水量補給飽和帶，通過深層滲漏水量將非飽和帶水分運動方程與飽和帶水分運動方程進行耦合。公式(6)中的V|Z1便為非飽和帶補給飽和帶的深層滲漏水量，其值可用公式(1)從地表Z0沿著垂向積分至潛水水面Z1，并利用Darcy定理求得：

(7)

將公式(7)代入公式(6)中，可得農田排水過程系統模型：

(8)

2.2 定解條件

2.2.1 初始條件

h(x,y,z,t)=h0(x,y,z)

(9)

(Z0

h(x,y,t)=hs(x,y)

(10)

(Z1

2.2.2 邊界條件

(11)

(z=Z0,t≥0)

h(z,t)=0 (z=Z1,t>0)

(12)

(13)

式中：E為作物蒸散發強度，cm/min；q=q0/A，q0為滴頭流量，mL/min,A為土壤濕潤面積，cm2；Ω為計算區域； ?Ω為水平四周邊界；n為(x,y)∈Ω的外法向量；g(x,y)為水平向排水溝流動的水流通量。

2.3 模型的求解

2.3.1 離散方程構建 利用Du Fort-Frankel有限差分格式對非飽和帶水分運動方程和飽和帶水分運動方程進行空間和時間離散，建立矩陣求解方程。采用矩形網格對空間進行劃分，非飽和帶含水率在垂向上急劇變化，需要較細的空間劃分網格；飽和帶空間尺度較大，采用較大的空間劃分網格。

解 法 四 ：由f（x）=sin2x+a cos2x= （2x+φ）的周期為π，并于直線x=-π對稱，可知f8，則，即a=-1.

2.3.2 程序編寫和計算 計算程序編寫采用的是MATLAB軟件。首先，對非飽和帶進行了網格劃分和離散，形成非飽和帶有限差分矩陣方程；再通過對飽和帶進行網格劃分和離散，形成飽和帶有限差分矩陣方程，根據耦合關系形成飽和帶-非飽和帶統一矩陣方程，最后進行求解。具體計算流程見圖1。

圖1 程序設計流程圖

首先輸入節點、步長、定解條件和參數等信息，形成飽和帶和非飽和帶水流矩陣方程。假設淺層地下水初始水頭為h1，確定非飽和帶垂向上的節點數，利用非飽和帶求解矩陣求得非飽和帶土壤含水率分布，再利用飽和-非飽和統一矩陣方程求得當前時刻淺層地下水水頭h2，若h2與h1收斂，說明h2為當前時刻水頭值；若h2與h1不收斂，則將h2重新賦值給h1，重復上一步驟，不斷迭代，直到h2與h1收斂為止，然后進行下一時刻的運算，直到運行結束。

3 模型驗證

3.1 實測數據獲取

在天津靜海建立了精準滴灌農田排水過程淺層地下水位的田間觀測系統。所選區域為45 m×95 m大小的四方形地塊，地塊的三周邊界為排水溝渠，溝渠寬2 m；另一邊界為水泥道路。所選地種植冬小麥，在地塊上布置滴灌帶，滴灌帶間距0.6 m，滴頭間距0.3 m，滴頭流量1.38 L/h。所選地塊及滴灌帶布置見圖2。

3.1.1 土壤基本特性 田間觀測區域土壤類型為黏壤土，土壤的非飽和導水率K(h)與壓力水頭h為高度非線性函數，可采用Van-Genuchten[19]模型所給出的表達式。土壤的飽和導水率Ks采用環刀法進行測定，比儲水系數 采用田間試坑法測定。具體參數見表1。

圖2 所選地塊及滴灌帶布置示意圖(單位：m)

土壤類型θr/(cm-3·cm-3)θs/(cm-3·cm-3)α/cm-1nKs/(cm-3·min-1)μ黏壤土0.0580.6360.00891.5770.00470.042

注：θr為土壤殘余含水率；θs為土壤飽和含水率；α為土壤進氣值的倒數；n為形狀系數；Ks為土壤飽和導水率；μ為土壤比儲水系數。

3.1.2 觀測孔布置 在精準滴灌灌溉條件下，灌溉水分在整個地塊上均勻分布，因此在相同土壤地質條件下，整個地塊土壤的地下水分運動是軸對稱的，可選地塊的一半作為觀測區域進行觀測點布置。橫向、縱向分別均勻布置10行(列)觀測點，共布置100個觀測點。

所選地塊地下水位觀測孔布置方式及編號如圖3所示。

圖3 所選地塊地下水位觀測孔布置示意圖(單位：m)

觀測孔編號：從下到上為1～10，從左到右為1～10，即左邊第1列的最下方的點編號為NO 101，從下往上為NO 101～NO 110；右邊第1列的最下方的點編號為NO 1001，從下往上為NO 1001～NO 1010；中間各點編號以此類推。

3.1.3 滴灌灌溉頻率和灌溉水量 根據Kang等[6～7]提出的滴灌農田精準水鹽調控方法，將滴頭正下方20 cm深度處的土壤基質勢閾值下限定為-20 kPa進行灌溉時，隨灌溉水進入土壤中的鹽分得到有效淋洗，土壤鹽分平衡，作物獲得高產。因此在整個地塊內選擇3個有代表性的點分別在距離滴頭正下方20 cm深度處埋設負壓計，當負壓計顯示的土壤基質勢低于-20 kPa時，對地塊進行灌水，每次灌水量為10 mm。

3.1.4 觀測內容 每天早上7:00記錄觀測孔潛水水位、排水溝水位和冠層上20 cm蒸發皿水量蒸發量。在非飽和帶水分運動過程中，需要考慮作物的蒸散發對水分運動的影響。目前測定作物蒸散發量的方法主要有水量平衡法、梯度法、蒸滲儀法、波文比法和渦度相關法等，這些測定方法所需參數較多，對參數要求精度高[20]。在田間實際觀測過程中，由于實驗條件有限，難以應用以上方法獲取日尺度上的作物實際蒸散發量。Liu等[21]通過在華北地區對冬小麥蒸散發量的研究發現作物冠層頂部蒸發皿的蒸發量與作物實際蒸散發量有良好的經驗關系，兩者的回歸系數趨近于1，根據冠層頂部蒸發皿的蒸發量能夠較為準確地估算冬小麥的日蒸散發量。在已有的實驗條件下，本文采用該方法對作物實際蒸散發量進行估算。

在拔節到灌漿期，冬小麥蒸散發量與冠層上20 cm蒸發皿水量蒸發量的關系式為[21]:

ETc=1.124ET20

(14)

式中：ETc為冬小麥蒸騰蒸發量,mm;ET20為冠層上20 cm蒸發皿蒸發量,mm。

圖4為冠層上20 cm蒸發皿蒸發量在觀測期間內隨時間的變化過程。

圖4 冠層上20 cm蒸發皿蒸發量隨時間的變化

3.2 模擬結果與實測結果的對比與分析

將精準滴灌農田排水過程系統模型潛水水位模擬值與實測值進行比較，初步識別和驗證模型。田間地下水位觀測一共持續了30 d，系統模型模擬運行時間分為兩個時間段：第3 d到第20 d和第21 d到30 d。因為在田間觀測的第3 d和第21 d有強降雨，降雨量分別為62和48 mm。

3.2.1 灌溉周期內淺層地下水位的變化和模擬結果 任取田間觀測1個灌水周期內的潛水水位與模擬水位值進行對比分析。圖5為灌溉周期內淺層地下水位實測值變化過程，圖6為灌溉周期內淺層地下水位模擬值變化過程，圖7為灌溉周期內淺層地下水位模擬值與實測值之差的變化過程。其中田間觀測的第7、8、9 d正好為1個灌溉周期，且沒有降雨。從圖5、6中可以發現，實驗區地塊靠公路端中部區域潛水水位最高，并沿著向排水溝方向逐漸降低，即離排水溝越近，潛水水位越低。潛水水位降低梯度沿著農田中間區域至排水溝方向逐漸減小，農田中間區域水位降低速度快，而排水溝渠附近區域水位降低速度慢。灌溉周期內的第1 d水位最高，并且隨著時間的推移，水位逐漸降低且水位降低速率逐漸減小，即灌溉后隨時間推移農田各點水位逐漸趨于穩定。結合圖5和6，系統模型潛水水位模擬值與實測值在時間和空間上的變化趨勢一致。

圖7中，系統模型淺層地下水位模擬值與實測值之差很小，隨時間推移，兩者的誤差由正值逐漸向負值轉化，誤差絕對值有增大趨勢。中間區域水位誤差大于兩側區域水位誤差，受邊界條件約束，離排水溝渠越近，水位模擬值越趨近于實測值。表2、3、4為灌溉周期內部分觀測孔潛水水位模擬值與實測值的具體數值。表2、3、4中潛水水位模擬值與實測值的絕對誤差在-13～21 mm范圍內，潛水水位模擬值與實測值吻合效果較好。

3.2.2 降雨后淺層地下水位的變化和模擬結果 取NO 101、NO 501、NO 505這3個地下水位觀察孔的潛水水位模擬值與實測值進行擬合分析。這3個觀測孔分別位于地塊公路端中間、中部中間和中部側邊，能一定程度上反映整個農田的模擬情況。圖8為NO 101、NO 501、NO 505觀測孔淺層地下水位模擬值與實測值的擬合結果。從圖8中發現，第2次降雨后的潛水水位下降速度明顯高于第1次降雨后的下降速度，這是因為第2次降雨后的作物蒸散發和溝渠水面蒸發速度高于第1次降雨后的蒸發速度。在模擬初始階段，淺層地下水位模擬值與實測值擬合效果。

圖5 淺層地下水位實測值(圖中展示的是整個實驗區域的一半，左邊為整個研究區域的中軸線，下圖類同)

圖6 淺層地下水位模擬值

圖7 淺層地下水位模擬值與實測值之差

很好，隨時間推移模擬值與實測值的誤差有增大趨勢。部分原因是模型在處理作物蒸散發量時采用的是一個階段的平均值，所以在第1個時間段先模擬值大于實測值，而后模擬值小于實測值；而第2個時間段為模擬值一直大于實測值。在模型運行第1個時間段，NO 101觀測孔潛水水位模擬值與實測值的絕對誤差在-27～43 mm之間；NO 501觀測孔潛水水位模擬值與實測值的絕對誤差在-22 ～45 mm之間；NO 505觀測孔潛水水位模擬值與實測值的絕對誤差在-27 ～47 mm之間。在模型運行第2個時間段，NO 101觀測孔潛水水位模擬值與實測值的絕對誤差在-34 ～-77 mm之間；NO 501觀測孔潛水水位模擬值與實測值的絕對誤差在-27 ～-76 mm之間；NO 505觀測孔潛水水位模擬值與實測值的絕對誤差在-27 ～-79 mm之間?？傮w上淺層地下水位模擬值與實測值擬合效果較好。雖然隨著模型運行，潛水水位模擬值與實測值的誤差有增大趨勢，但是對于強降雨后短期內的農田排水過程，系統模型能較好地模擬精準滴灌農田的潛水水位變化過程。

表2 第7 d部分觀察孔水位實測值與模擬值 m

表3 第8 d部分觀察孔水位實測值與模擬值 m

表4 第9 d部分觀察孔水位實測值與模擬值 m

注：表2～4中括號外的值為潛水水位模擬值，括號內的值為潛水水位實測值。

圖8 NO 101、NO 501、NO 505觀測孔地下水位模擬值與實測值擬合結果圖

4 結 論

(1)本文根據精準滴灌農田土壤水分的分布特征，利用Richards方程建立了基于有限差分法的精準滴灌農田排水過程系統模型，模擬了精準滴灌農田排水過程。

(2)潛水水位降低梯度沿著農田靠公路端中間區域至排水溝方向逐漸減小，農田中間區域水位降低速度快，而排水溝渠附近區域水位降低慢。灌溉周期內的第1 d水位最高，并且隨著時間的推移，水位逐漸降低且水位降低速率逐漸減小，即灌溉后隨時間推移農田各點水位逐漸趨于穩定。

(3)在降雨后農田排水過程的模擬中，潛水水位模擬值與實測值在時間和空間上變化趨勢一致，兩者差異較小，吻合較好。雖然隨著模擬時間變長，潛水水位模擬值與實測值的差異有增大趨勢，但總體上，系統模型能較好地真實反映強降雨后短期內精準滴灌農田的排水過程，可為大面積應用滴灌技術的現代灌區排水系統規劃設計和改造以及灌排信息化管理提供參考依據。