基于分層水量平衡模型的蓄水坑灌果園土壤水分動態模擬

張甜甜，郭向紅，2，雷濤，孫西歡，馬娟娟，鄭利劍

（1.太原理工大學水利科學與工程學院，太原 030024；2.中國水利水電科學研究院流域水循環模擬與調控國家重點實驗室，北京 100038）

0 引 言

土壤水分與作物生長息息相關，農田土壤水分的動態變化直接影響作物的生長發育，農田土壤水分含量長期處于過高或過低的水平都會對作物的生長產生不利的影響，因此對田間土壤水分含量進行動態監測與預測對保障農作物的產量具有重要的意義[1]。

對于土壤水分動態的獲取方法主要有儀器監測和模型模擬。儀器監測采用專有土壤水分測量儀器長期對田間水分進行測量，比較費時費力，監測成本較高。而近些年隨著計算技術的不斷發展，模型模擬逐漸成為土壤水分動態研究的重要手段。常見的土壤水分動態模型有三類，即機理模型、神經網絡模型、水量平衡模型。機理模型[2-4]是根據土壤水分運動方程并結合實際邊界條件，建立不同條件下的土壤水分運動模型，可以模擬土壤水分的空間分布和動態變化，但由于模型為非線性偏微分方程，數值模型的收斂性較差，需要精確水分運動參數，計算較為復雜，比較費機時。神經網絡模型[5-7]是目前較為流行的基于數據驅動的預測方法，可以實現對土壤水分的實時動態預測，但該類模型屬于黑箱模型，不能解釋模型背后的機理，且需要大量的田間數據支撐，才能保證預測精度。水量平衡模型[8]是根據根區土壤水量平衡原理，建立的反映根區水量動態變化的模擬模型，原理簡單，所需數據較少，雖然不能精準模擬土壤水分空間分布，但模擬的根區土壤水分動態變化可以滿足農田水分管理的需求，已經廣泛應用于農田水分動態變化的研究。

目前土壤水量平衡模型結構已由單層水量平衡模型發展到多層水量平衡模型[9-11]，主要涉及的作物有冬小麥、夏玉米、牧草等[12-14]，模型模擬的灌水方法有地面灌溉包括漫灌、溝灌、噴灌等全面灌水方法，但對于如何模擬局部灌溉對土壤水分的影響尚未研究。蓄水坑灌是一種典型的適用于果林的局部灌溉方式，水分通過坑壁側向滲入根區土壤，入滲界面并非普通灌溉方式的垂向入滲，而是側向入滲[15]，因此現有的分層水量平衡模型無法直接應用于蓄水坑灌土壤水分動態模擬，需要對現有的模型進行改進。基于此，本文將考慮蓄水坑灌的水分分布特征，對現有分層水量平衡模型進行改進，建立蓄水坑灌條件下的分層水量平衡模型，以期為蓄水坑灌果園土壤水分動態預測提供新途徑。

1 材料和方法

1.1 試驗區概況

試驗于山西省晉中市太谷縣農業科學院果樹研究所（112°32′E，37°23′N）進行，該地平均海拔約781.9 m，無霜期175 d，平均氣溫一般在5～10 ℃之間。全年降雨量分布不均勻，全年降雨量最多的時期為6-8月，降雨量最少的時期為12-2月，多年平均降雨量約為460 mm，氣候類型屬于典型的溫帶性大陸氣候。

依據試驗區土壤質地對土體進行分層，將地表以下200 cm的土體分為3層，第一層厚度為70 cm，第二層厚度為50 cm，第三層厚度為80 cm。土壤的主要物理參數如表1所示。

表1 土壤物理參數表Tab.1 Table of physical parameters of soils

1.2 試驗設計

試驗所用果樹品種為紅富士，南北方向種植，行距4 m，株距2 m。試驗用樹選擇樹齡與樹勢基本一致、無病蟲害的蘋果樹，并在樹下布置蓄水坑，蓄水坑的直徑為30 cm，深度為40 cm。本試驗共設2個處理，分別為處理一（灌水上下限為田間持水量的50%～80%）和處理二（灌水上下限為田間持水量的60%～90%），每個處理3個重復。試驗于2015年4-10月進行，如表2所示為試驗期間灌水安排的具體情況。試驗期間每隔7～10 d對各層土壤含水率進行監測，并在降雨、灌溉前后進行加測。

表2 灌水時間及灌水定額Tab.2 Irrigation time and irrigation quota

1.3 測試項目

1.3.1 土壤水分的測定

各測點土壤體積含水率采用TRIME-PICO IPH土壤水分測量系統進行監測，監測范圍為果樹的樹干周圍地表以下0～200 cm深度。

土壤水分測點分過坑、不過坑、處于過坑與不過坑的中間位置3種情況。過坑測點布置于樹干與蓄水坑中心連線上，距離樹干分別為30 cm（g30）、50 cm（g50）、100 cm（g100）、120 cm（g120）；不過坑測點布置于相鄰兩蓄水坑分界線與樹干連線上，距離樹干分別為30 cm（b30）、60 cm（b60）、90 cm（b90）、120 cm（b120）；中間位置測點距離樹干分別為90 cm（z90）、120 cm（z120）。在垂直方向上每20 cm取一個水分測點，深度取到地表之下200 cm處，水平方向的水分測點布置如圖1所示。

圖1 含水率監測點位置示意圖（單位：cm）Fig.1 Schematic diagram of water content monitoring points

1.3.2 葉面積指數的測定

在試驗期間每隔一周采用LAI-2200冠層分析儀對葉面積指數進行測量。

1.3.3 氣象數據的測定

本次試驗采用精度較高的小型自動氣象監測站對試驗地的參考作物需水量ET0以及降雨量進行了實時監測，圖2為試驗期間ET0及降雨量變化圖。

圖2 試驗期間ET0及降雨量變化圖Fig.2 ET0 and rainfall variation during the experiment

2 模型建立

2.1 蓄水坑灌果園土壤分層水量平衡模型的建立

土壤水量平衡模型根據一定時段內水分的輸入和輸出來確定土壤水分的動態變化，反映了作物根系層水量變化與水分收、支之間的關系。根據蓄水坑灌的布置以及土壤質地的分層情況，將0～200 cm土層分為上、中、下3層，如圖3所示為蓄水坑灌果園土壤分層水量平衡模型概化圖。第一層土壤的水分收入有降雨量P去除植被截留量C之后的凈雨量(PC)、灌溉量I和第二層土壤分配至第一層的水量q1,2，水分支出有地表蒸發量Es、蓄水坑的坑壁蒸發量Ep、下滲至第二層的水量Q1、根系吸水量S1；第二層土壤的水分收入有第一層的下滲水量Q1和第三層土壤分配至第二層的水量q2,3，水分支出有下滲至第三層的水量Q2、分配至第一層的水量q1,2和根系吸水量S2；第三層土壤的水分收入為第二層的下滲水量Q2，水分支出有第三層的下滲水量Q3、分配至第二層的水量q2,3和根系吸水量S3。不考慮地下水補給，結合土壤水量平衡模型可得蓄水坑灌果園的各層土壤水量平衡方程如下：

圖3 蓄水坑果園土壤水量平衡模型概化圖Fig.3 Soil water balance model generalization of water storage pit orchard

由各層土壤的儲水量變化可得出每一層的土壤含水率計算公式：

2.2 下滲水量

蓄水坑灌是一種典型的局部灌溉技術，在布設蓄水坑并進行灌水之后，水分將通過坑的側壁滲入土壤，土壤水分呈現出非均勻特征[16]，且表層土壤可能沒有濕潤或者僅很小部分濕潤。因此，蓄水坑灌條件下降雨和灌溉之后的水分運移情形是不同的，降雨之后水分從地表開始滲入土壤，逐漸濕潤下層土壤；灌溉時水分側向入滲，并且水分并未完全濕潤上層土壤就向下運動。因此在本模型中對于第一層土壤的下滲水量而言，降雨和灌溉需要分開考慮。

2.2.1 降雨之后的下滲水量

雨水經過植被攔截之后由地表向下層土壤入滲，根據滿溢滲透的假定[17]，當土層中的水分達到田間持水量之后再向下滲透，直到降雨量和下滲水分滲透完為止。降雨之后超出飽和含水量的部分全部下滲，由于實際情形中即使沒有達到田間持水量，也會發生土壤水分的下滲，所以本文認為未超出飽和含水量的部分當其達到臨界含水量[18]時就發生下滲。

在實際計算下滲水量時，分為兩種情況，一種是該層土壤的水分增加量超過該層的容水量即該層所能保持的水分（時段初土壤儲水量與飽和含水率對應的土壤儲水量的差值），此時各層的下滲水量公式計算如下[18,19]：

式中：Pi為第i天的降雨量，cm；為第i天第j-1層土壤下滲水量，cm；為第i天第j層土壤容水量，cm；n為土壤剖面的排水常數[18]，取值范圍一般為0.30～0.35；θsj為第j層土壤的飽和含水率，cm3/cm3；θj為入滲之前第j層土壤的含水率，cm3/cm3；θfj為第j層土壤的田間持水率，cm3/cm3；LAI為作物葉面積指數；LAImax為某一時段內作物葉面積指數的最大值；θkj為第j層土壤的臨界含水率（當土壤水分低于某一臨界點時，導致作物的生長狀況發生顯著變化時所對應的土壤含水率[20]；θkj與作物葉面積指數及各層土壤的田間持水率呈線性關系[18]），cm3/cm3。

另一種是該層土壤的水分增加量并未超過該層的容水量，若此時該層的土壤含水率小于該層土壤的臨界含水率時則下滲水量為0，當該層的土壤含水率大于臨界含水率時按下式計算下滲水量：

2.2.2 灌溉之后的下滲水量

灌水之后，蓄水坑中的水并未完全濕潤表層就向下運動，不滿足滿溢滲透假定，為了使模型反映此種現象，假定第一層土壤所能保持的水分即容水量減少，即第一層入滲水量超過容水量就向下層土壤入滲。因此，本文引入局部灌溉系數ξ對灌溉之后第一層土壤的容水量進行了修正，修正之后第一層土壤的容水量為ξM1(ξ＜1)。在計算灌溉之后土壤的下滲水量時只用修正第一層土壤0～70 cm土層的容水量，如公式（7）所示，剩余兩層土壤的下滲水量仍按照降雨時的下滲水量公式進行計算。

2.3 土壤重分配水量

由于各層的土壤含水率差異較大，需考慮在相鄰土壤層間擴散的水量，引入下述土壤水重分配公式[19]對重分配的水量進行計算：

式中：qj,j+1為第(j+1)層的水分再分配量（以水分向上擴散為正），cm；θwj為第j層土壤的凋萎含水率，cm3/cm3；ηj為土壤水分擴散系數，與含水量有關，按下式計算：

植物不能從土壤中吸水時的土壤含水率稱為凋萎含水率[21]，本文蓄水坑灌果園土壤凋萎含水率對應的負壓水頭為-150 m[22]。

2.4 蒸發量

蓄水坑灌法是將灌溉水注入樹冠下的若干個蓄水坑內，比起傳統灌溉多了坑壁蒸發，蓄水坑灌的棵間蒸發包括地表蒸發和坑壁蒸發。在蓄水坑灌條件下，蒸發強度主要與土壤含水率因素有關，本文采用谷瓊瓊[23]擬合得到的蓄水坑灌果園地表蒸發強度和坑壁蒸發強度計算公式：

式中：es和ep分別為蓄水坑灌地表蒸發強度和坑壁蒸發強度，mm/d，通過已知初始含水率計算得到初始的地表蒸發強度和坑壁蒸發強度。

2.5 植被截留量

當發生降雨時，部分雨水經過植物的葉片時會被截留，這部分雨水不會滲入土體，超過截留能力的雨水才會到達土壤表面。植被的最大截留降水量指的是單位面積上植被所能截留的最大降水量，即聚積在葉面所在水平面上的水層深度。植被截留量[24]的大小主要與植被的葉面積指數LAI和降雨強度有關。葉面積指數為單位面積上的植物葉面積總和，即葉面積總和與占地面積的比值。本文采取馬生軍等人的蘋果樹冠層降雨截留模型[25]：

式中：C為截留量，mm；P為計算時段內的降雨量，mm。

2.6 根系吸水量

本文采用的根系吸水模型為[26]：

式中：S(x,y,z,t)為根系實際吸水量，L；Vmax(x,y,z,t)為根系最大吸水速率，L/min；γ(h)為水分脅迫系數。

γ(h)可表示為：

式中：h0、h1、h2、h3為θs、80%θf、60%θf、θw對應的基質勢水頭。Vmax(x,y,z,t)采用Vrugt[27]的表述。

式中：Tpot為潛在蒸騰強度，cm/min；Kc為作物系數，采用FAO[28]發布的作物蒸散量、作物需水量計算指南中的數據；ET0為參考作物需水量，cm/min，采用氣象監測站采集的數據。

依據上述根系吸水模型和氣象站收集到的數據，采用MATLAB進行編程，計算得出模擬區域各層各個節點處的根系吸水速率，根據求得的節點處的平均根系吸水速率進而求得各層的根系吸水速率以及根系吸水量。

2.7 模型精度評價指標

精度評價的指標為均方根誤差RMSE、平均相對誤差MAPE和平均絕對誤差MAE，計算公式如下：

3 結果分析

3.1 模型率定與驗證

將2015年4月13日測得的土壤含水率作為計算的初始含水率，采用上述的模型對蓄水坑灌條件下4-10月份的田間土壤含水率進行模擬。采用處理一的數據率定局部灌溉系數ξ，并采用處理二的數據對模型進行驗證。系數ξ的率定過程為，采用分層水量平衡模型，分別計算ξ為0.6、0.7、0.8、0.9、1.0時的土壤含水率，通過精度評價比較得出在本模型下參數ξ的最優值，表3為不同ξ取值下處理一模擬計算精度評價表。

表3 不同ξ取值下處理一模擬計算精度評價表Tab.3 Evaluation of the accuracy of treatment 1 simulation calculation under different values of ξ

從表3的計算精度評價表來看，隨著ξ的取值從0.6變化到1.0，模型的計算誤差先減小后增大，當ξ=0.8時的模擬效果最好，所以ξ的率定值為0.8。同時由表3可知，模型各個土層深度的計算誤差隨深度增加基本呈先減小后增大的趨勢，70～120 cm土層模擬精度最高。這可能是由于0～70 cm土層直接受到降雨、灌溉、地表蒸發等氣象條件和蓄水坑邊界等多因素影響，造成該層模型模擬的精度比70～120 cm低；120～200 cm土層在垂向上處于模型模擬的最深層，在模型中只考慮了70～120 cm 土層和120～200 cm土層之間的水分交換，忽略200 cm以下區域與120～200 cm土層之間的水分交換，造成該層模擬精度比70～120 cm土層低，但從表3的誤差指標來看，忽略這一部分的水量雖然模擬精度降低但在可接受范圍之內。進一步將ξ=0.8代入模型，進行處理二模擬計算，采用處理二的實測數據對模型進行驗證，表4為處理二模型計算精度評價表。

表4 處理二模擬計算精度評價表Tab.4 Evaluation of the accuracy of treatment 2 simulation calculation

從表4來看，當設定局部灌溉系數ξ為0.8時，處理二各土層土壤含水率模擬值與實測值的3種誤差都較小，模型模擬的精度較高，可以用于蓄水坑灌土壤水分動態模擬。

3.2 土壤水分動態變化對比分析

圖4和圖5是2015年4-10月蓄水坑灌果園的分層水量平衡模型模擬的土壤含水率隨時間變化的動態。表5和表6是灌水和降雨之后各層土壤含水率模擬值與實測值的增加幅度情況。表7和表8是兩種處理下各層土壤含水率的極差值以及平均含水率值。

圖4 處理一各土層土壤含水率模擬值與實測值比較Fig.4 Comparison of simulated and measured soil water content values for each soil layer in Treatment 1

圖5 處理二各土層土壤含水率模擬值與實測值比較Fig.5 Comparison of simulated and measured soil water content values for each soil layer in Treatment 2

表5 灌水之后各層土壤含水率模擬值與實測值的增加幅度Tab.5 Increase in simulated and measured soil water content of each layer after irrigation

表6 降雨之后各層土壤含水率模擬值與實測值的增加幅度Tab.6 Increase in simulated and measured soil water content of each layer after rainfall

表7 各層土壤含水率極差值Tab.7 Value of extreme difference of soil water content in each layer

表8 各層土壤平均含水率 cm3/cm3Tab.8 Average water content of soil in each layer

從圖4和圖5土壤含水率的動態變化來看，兩種處理條件下3層土壤含水率模擬值與實測值的變化趨勢較為一致，結合表3和表4的計算精度評價指標表明實測值與模擬值吻合較好。兩個水分處理的水分動態模擬結果表明，由于蒸發、水分下滲、果樹的根系吸水等耗水因素，隨著時間的推移各層土壤含水率呈下降趨勢，當有灌溉和降雨時呈上升趨勢，總的來說各層土壤含水率在一定范圍內不斷波動。在灌水日期附近，兩種處理下各層土壤含水率均出現了較大幅度的增加，符合實際情況。處理一在模擬期間于5月19日和7月12日進行了兩次灌水，處理二在模擬期間于5月4日、6月20日和8月26日進行了3次灌水，灌水之后各土層土壤含水率的模擬值與實測值分別較上個測量時間的增加幅度如表5所示。綜合表5的數據來看，灌水之后70～120 cm和120～200 cm土層土壤含水率模擬值與實測值的增加幅度明顯大于0～70 cm土層，且70～120 cm土層的增加幅度最大，這是由于蓄水坑灌是中深層立體的局部灌溉方式，灌溉水更集中分布在整個土體的中層，這也是果樹根系分布相對密集的區域。說明采用蓄水坑灌的灌溉方式可以有效減少地表蒸發，給果樹生長提供更充足的水分。在整個模擬期間降雨量較大的時段出現在9月8日-10日，累計降雨量為57.4 mm，這次降雨之后兩種處理下各層土壤含水率模擬值與實測值也出現了較大幅度的增加，具體增加幅度如表6所示。從表6的數據來看，降雨之后0～70 cm土層土壤含水率模擬值與實測值的增加幅度明顯大于70～120 cm和120～200 cm土層，這是由于降雨增加的水量是直接從地表開始入滲的，因此上層土壤含水率的增加幅度大于中下層土壤。在蓄水坑灌條件下灌溉和降雨之后0～70 cm土層的下滲機理不同，因此灌溉和降雨發生之后從淺到深各土層土壤含水率的增加幅度呈現了不同的趨勢，具體表現為灌溉之后從淺到深各層土壤含水率增加幅度呈先增大后減小的趨勢，降雨之后從淺到深各層土壤含水率的增加幅度呈逐漸減小的趨勢。從表5及表6的數據來看，灌水和降雨之后各土層土壤含水率實測值的增加幅度通常低于模擬值，這是由于實際水分的變化具有一定的滯后性，沒有模型模擬對灌水和降雨的響應迅速。

從3個土壤層水分的動態模擬情況和表7的數據來看，最接近地表的0～70 cm土層由于受到外界氣象因素的影響較大，比如降雨、灌溉以及地表蒸發和坑壁蒸發，相比較其他土層，0～70 cm土層的含水率變化趨勢波動幅度較大，70～120 cm土層以及120～200 cm土層含水率波動趨勢相較上層土壤更為平緩。從表8的數據來看，模擬期間兩個處理的平均含水率均表現為0～70 cm土層最低，這是因為表層直接受到太陽輻射，地表蒸發和坑壁蒸發較為強烈。處理二的灌水上限高于處理一，且處理二的灌水次數多于處理一，處理二各土層的土壤平均含水率均高于處理一。

4 討 論

土壤分層水量平衡模型已廣泛應用于農田各層土壤水分動態的模擬，陳皓銳[19]、張雪飛[29]等人根據水量平衡原理建立了土壤分層水量平衡模型，模擬了各土層土壤含水率的動態變化。但現有的分層水量平衡模型無法直接應用于蓄水坑灌土壤水分動態模擬，本文基于土壤的水量平衡原理，在現有模型的基礎上引入局部灌溉系數ξ，建立了蓄水坑灌條件下0～200 cm土層的分層水量平衡模型，通過精度評價，表明模型的精度較高，證實本文的分層水量平衡模型可適用于蓄水坑灌果園土壤水分動態模擬。局部灌溉系數ξ的合理取值十分關鍵，ξ的大小會影響到模型的下滲水量，在模擬計算前需要對此參數進行率定。

由于0～70 cm土壤與外界的水分交換較為頻繁，本模型模擬的該層土壤含水率變化對外界的響應更敏感，波動幅度更明顯，這與劉雪靜[30]等人的研究一致。由于測量儀器在淺層更容易產生誤差以及灌溉和降雨會引起含水率的突變，導致該層土壤的模擬誤差較大。表層土壤直接受到日光照射和溫度的影響，水分流失較快，模型模擬易高估0～70 cm土層土壤含水率，這與劉昭[10]等人的研究結果一致。下層土壤最接近地下水，在實際中即使超過了田間持水率也有可能并不會立即發生水分下滲，因此在用本模型模擬時會出現低估120～200 cm土層土壤含水率的現象，這與郭其樂[9]等人的研究結果一致。當長時間無降雨或者灌溉，無明顯的水分收入之后，模型的模擬能力會下降[9]，在降雨以及灌溉之后，模型的水分收支情況更新，會使各層土壤水分不斷達到新的平衡，使模型的模擬效果更好。

5 結 論

本文在現有分層水量平衡模型的基礎上引入局部灌溉系數ξ，建立了蓄水坑灌果園的分層水量平衡模型，得到以下結論：

（1）采用處理一土壤含水率模擬值與實測值的數據率定了局部灌溉系數ξ的值為0.8，采用處理二的數據對模型進行了驗證，精度評價表明模型的模擬精度較高，可以用于蓄水坑灌土壤水分動態模擬。

（2）在蓄水坑灌條件下，灌溉之后從淺到深各土層土壤含水率的增加幅度呈先增大后減小的趨勢，水分主要集中在中深層土壤深度范圍內；降雨之后從淺到深各土層土壤含水率的增加幅度呈逐漸減小的趨勢，水分主要集中在上層土壤。

（3）兩個水分處理的水分動態模擬結果表明，各土層土壤含水率在一定范圍內不斷波動，0～70 cm土層的土壤含水率波動范圍最大。