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

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

1 研究背景

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

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

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

2 精準滴灌農(nóng)田排水過程的數(shù)學(xué)模型

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

2.1 模型的控制方程

2.1.1 非飽和帶水分運動控制方程 在滴灌條件下，當(dāng)?shù)喂鄮Ш偷晤^間距較小時，滴灌帶間濕潤鋒會發(fā)生重疊，形成濕潤帶，土壤水分運動為三維運動。假設(shè)土壤為均質(zhì)、各向同性的剛性多孔介質(zhì)，不考慮溫度、生物或化學(xué)反應(yīng)對水分運動的影響，也不考慮土壤水分運動的滯后效應(yīng)，土壤非飽和帶水分運動可用三維Richards方程描述:

(1)

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

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

(2)

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

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

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

(3)

(4)

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

(5)

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

(6)

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

(7)

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

(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為作物蒸散發(fā)強度，cm/min；q=q0/A，q0為滴頭流量，mL/min,A為土壤濕潤面積，cm2；Ω為計算區(qū)域； ?Ω為水平四周邊界；n為(x,y)∈Ω的外法向量；g(x,y)為水平向排水溝流動的水流通量。

2.3 模型的求解

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

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

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

圖1 程序設(shè)計流程圖

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

3 模型驗證

3.1 實測數(shù)據(jù)獲取

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

3.1.1 土壤基本特性 田間觀測區(qū)域土壤類型為黏壤土，土壤的非飽和導(dǎo)水率K(h)與壓力水頭h為高度非線性函數(shù)，可采用Van-Genuchten[19]模型所給出的表達式。土壤的飽和導(dǎo)水率Ks采用環(huán)刀法進行測定，比儲水系數(shù) 采用田間試坑法測定。具體參數(shù)見表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為土壤飽和含水率；α為土壤進氣值的倒數(shù)；n為形狀系數(shù)；Ks為土壤飽和導(dǎo)水率；μ為土壤比儲水系數(shù)。

3.1.2 觀測孔布置 在精準滴灌灌溉條件下，灌溉水分在整個地塊上均勻分布，因此在相同土壤地質(zhì)條件下，整個地塊土壤的地下水分運動是軸對稱的，可選地塊的一半作為觀測區(qū)域進行觀測點布置。橫向、縱向分別均勻布置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 滴灌灌溉頻率和灌溉水量 根據(jù)Kang等[6～7]提出的滴灌農(nóng)田精準水鹽調(diào)控方法，將滴頭正下方20 cm深度處的土壤基質(zhì)勢閾值下限定為-20 kPa進行灌溉時，隨灌溉水進入土壤中的鹽分得到有效淋洗，土壤鹽分平衡，作物獲得高產(chǎn)。因此在整個地塊內(nèi)選擇3個有代表性的點分別在距離滴頭正下方20 cm深度處埋設(shè)負壓計，當(dāng)負壓計顯示的土壤基質(zhì)勢低于-20 kPa時，對地塊進行灌水，每次灌水量為10 mm。

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

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

ETc=1.124ET20

(14)

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

圖4為冠層上20 cm蒸發(fā)皿蒸發(fā)量在觀測期間內(nèi)隨時間的變化過程。

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

3.2 模擬結(jié)果與實測結(jié)果的對比與分析

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

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

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

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

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

圖6 淺層地下水位模擬值

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

很好，隨時間推移模擬值與實測值的誤差有增大趨勢。部分原因是模型在處理作物蒸散發(fā)量時采用的是一個階段的平均值，所以在第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上淺層地下水位模擬值與實測值擬合效果較好。雖然隨著模型運行，潛水水位模擬值與實測值的誤差有增大趨勢，但是對于強降雨后短期內(nèi)的農(nóng)田排水過程，系統(tǒng)模型能較好地模擬精準滴灌農(nóng)田的潛水水位變化過程。

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

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

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

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

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

4 結(jié) 論

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

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

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