負水頭環境土壤水分濕潤鋒運移三維模擬

朱 君，李 婷，謝 添，張艾明

（中國輻射防護研究院 核環境模擬與評價技術重點實驗室，太原 030006）

0 引 言

負水頭供水可以精準控制土壤含水率，維持土壤的非飽和狀態，避免因干濕交替使植物遭受旱澇脅迫和土壤養分淋溶損失，是開發新型節水灌溉技術的有效措施[1-2]。因此，對負水頭條件下土壤水分運動特征的研究，具有十分重要的理論和實踐意義。目前，大量的實驗結果表明負水頭條件下，累積入滲量與時間、最大水平和垂直濕潤距離與時間都符合冪函數關系，且濕潤體的寬深比隨入滲時間和水頭的增加而增加[3-6]。Rodriguez-Sinobas等[7]也認為濕潤體體積和寬深比隨著供水壓力的增加而增加，且入滲速度與供水壓力呈冪函數關系。以上研究不是建立在土壤水動力學的基礎上，是直接擬合試驗數據得到的經驗公式。如何結合試驗現象，確定負水頭環境下土壤水分運移模型，并預測不同供水水頭、不同土壤質地等條件下灌水器的間距，實現土壤含水率的持續精確控制則更加重要。但是現有模型都只適用于描述負水頭條件下土壤濕潤鋒一維垂直、水平運動規律[8-10]。鄒朝望等[11-12]基于Gardner和Russo公式表征了土壤水分特征曲線與非飽和土壤的導水率，推導了負水頭下濕潤峰、入滲率、累積入滲量與入滲時間的關系，并用試驗數據驗證了其合理性。范軍亮等[13-14]發現Philip模型和Kostiakov公式均可以較好描述負水頭條件下，4種典型黃土高原土壤的累積入滲量、濕潤鋒運移距離、入滲率和入滲時間之間的關系。上官玉鐸等[15-16]對比了負水頭條件下的土壤水分入滲規律，發現累積入滲量與時間呈冪函數關系，濕潤鋒與時間的平方根呈線性關系，Kostiakov公式比Philip模型、Green-Ampt模型能更準確描述入滲率和時間的關系。辛琛等[17-18]通過實驗獲得了不同負水頭下累積入滲量、濕潤鋒和入滲時間的關系，發現Philip模型、Green-Ampt模型和Kostiakov公式都可以描述負水頭條件下的土壤濕潤鋒運移規律。王佳佳等[1]將對土壤濕潤鋒運移規律的研究拓展到二維模型，利用室內實驗觀測了不同負壓水頭條件下土壤水鹽的運移特征，通過實驗數據驗證了Hydrus-2D模型的適用性，對不同質地土壤中水分分布、累積入滲量、水分入滲速率、鹽分分布和蒸發量進行了模擬。冀榮華等[19]結合土壤水動力學建立了負壓灌溉條件下的Hydrus-2D水分入滲模型，模擬結果與實測結果的相對誤差為2%～4%。Yao等[20]通過二維濕潤鋒實驗發現，濕潤體近似同心橢球面，且橢球中心含水率最高。現階段負水頭條件下的土壤濕潤鋒運移模型以一、二維為主[21-26]，三維模型鮮有報道。本文通過室內實驗數據驗證負水頭環境下Hydrus-3D土壤水分模型的準確性和有效性，對完善負水頭條件下的土壤水分運動理論十分必要。

1 材料與方法

1.1 儀器與試劑

實驗裝置由馬氏瓶、有機玻璃土箱、多孔板灌水器、橡皮輸水管組成。有機玻璃箱為扇形柱體，高40.0 cm，徑向長度35.0 cm，夾角30o；多孔板灌水器的出水孔為 3～4 μm（15 cm×10 cm×2.5 cm）；馬氏瓶高75.0 cm，內徑9.0 cm；橡皮輸水管內徑6.0 mm。試驗采用超純水（Milli-Q Element超純水制備機，日本Millipore公司）。

1.2 試驗樣品

選擇山西省榆次地區的砂土、壤土介質為試驗對象，用搖擺式篩析機（AS200型，德國Retsch公司）以及激光粒度分析儀（Mastersizer 3000E，英國Malvern公司）測定粒徑分布。其中，砂土的砂粒量85.50%、粉粒量 5.47%、黏粒量 9.03%，有機質量1.7 g/kg，土壤pH值8.34。壤土的砂粒量54.21%、粉粒量29.65%、黏粒量16.14%，有機質量9.3 g/kg，土壤pH值8.08。

1.3 土壤水力參數測定

測定砂土、壤土的水分特征曲線與飽和導水率。水分特征曲線采用離心機法測定，由低到高依次設置12個不同的轉速對飽和試樣進行脫水，得到不同體積含水率下的土壤吸力值。飽和導水率采用定水頭法測定，首先從土柱底端通水，使整個土柱體飽和，上部安裝定水頭裝置，調整入水口供水量，使試驗期間水頭差H保持不變，記錄時間Δt和底部出口流量Q。

1.4 試驗方法

砂土、壤土均按比例1.25 g/cm3分層裝入有機玻璃箱。將馬氏瓶、灌水器用橡皮管連接起來，埋入有機玻璃箱的土壤介質中，并保證整套裝置的密封性。以多孔板灌水器頂部為參考平面，設計3種不同的作用水頭H，如下：

1）H=0 m，作用水頭為 0，馬氏瓶進氣口與灌水器平。

2）H=-0.5 m，作用水頭為負壓，馬氏瓶進氣口比灌水器低0.5 m。

3）H=-1.0 m，作用水頭為負壓，馬氏瓶進氣口比灌水器低1.0 m。

記錄累積入滲量和濕潤鋒運移曲線。在有機玻璃箱的表面粘貼透明膠片記錄水分入滲過程，累積入滲量通過馬氏瓶上的刻度讀取。

2 結果與分析

2.1 土壤水分濕潤鋒運移結果

在有機玻璃箱表面粘貼的透明膠片上記錄相應時刻砂土、壤土的濕潤鋒運移過程，并在試驗結束后將膠片上的水分運移曲線轉化為Excel數據，見圖2、圖3。

圖2 砂土水分濕潤鋒運移曲線Fig.2 Wetting front moving curve of sand

圖3 壤土水分濕潤鋒運移曲線Fig.3 Wetting front moving curve of loam

在重力勢和基質勢作用下，濕潤鋒向垂直和水平方向運動，濕潤體近似為1/4橢圓狀，濕潤鋒包絡面積隨時間逐漸增加。

1）H=0 m時，砂土的濕潤鋒包絡面積在30 min到達904.1 cm2，壤土的濕潤鋒包絡面積在240 min到達為791.7 cm2。

2）H=-0.5 m時，3 480 min砂土的濕潤鋒包絡面積到達742.3 cm2，2 820 min壤土的濕潤鋒包絡面積到達745.7 cm2。

3）H=-1.0 m時，4 200 min砂土的濕潤鋒包絡面積到達241.2 cm2，4 200 min壤土的濕潤鋒包絡面積到達629.5 cm2。

2.2 入滲量與時間的關系

3種不同負水頭條件下，砂土、壤土的累積入滲量隨著時間的增加逐漸增大，結果見圖4。可以用冪函數I=m×tn描述累積入滲量與時間的關系，隨著負水頭高度的增加，m、n值均減小。砂土的m值由66.51減少至1.08，n值由0.96減少至0.42；壤土的m值由38.49減少至6.59，n值由0.69減少至0.56，擬合值見表1。與文獻[3-7]的試驗結果一致。

表1 累積入滲量與時間擬合值Table 1 Fitting parameters of cumulative infiltration and time

H=0 m，砂土m值是壤土的172.80%，原因是壓力勢占主導作用，且砂土的滲透系數大于壤土。H=-0.5 m，砂土m值是壤土的13.52%；H=-1.0 m，m值砂土是壤土的 16.34%，壤土的累積入滲量曲線高于砂土，基質勢逐漸占主導作用。

圖4 累積入滲量隨時間變化曲線Fig.4 Curve of cumulative infiltration versus time

2.3 濕潤峰與時間的關系

3種不同負水頭條件下，砂土、壤土濕潤鋒隨時間的增加逐漸向水平、垂直方向擴大，結果見圖5。水平和垂直最大濕潤距離與時間的關系，用下式描述：

式中：Zf為最大濕潤距離（mm）；t為時間（min）；p與q為常數，擬合關系見表2和表3。砂土和壤土的水平、垂直最大濕潤距離與時間的平方根呈良好的線性關系[15-16]，相關系數R2都大于0.96，p隨著負壓水頭的增加而減小，砂土p值由36.20減少至1.17，壤土p值由13.24減少至2.43。

H=0 m，砂土p值是壤土的273.41%，單位時間砂土濕潤鋒遷移距離大于壤土，原因是壓力勢占主導作用。H=-0.5 m，砂土p值是壤土的78.90%；H=-1.0 m，砂土p值是壤土的48.15%，單位時間壤土濕潤鋒遷移距離大于砂土，基質勢逐漸占主導作用。

表2 水平最大濕潤距離與時間擬合值Table 2 Fitting parameters of maximum horizontal wetting distance and time

表3 垂直最大濕潤距離與時間擬合值Table 3 Fitting parameters of maximum vertical wetting distance and time

圖5 濕潤鋒距離與時間變化曲線Fig.5 Curve of wetting distance versus time

圖6 濕潤鋒入滲速度與時間變化曲線Fig.6 Curve of infiltration velocity of wetting front versus time

2.4 入滲速度與時間的關系

負水頭高度從0 m增加至-1.0 m時，濕潤鋒入滲速度逐漸減小，見圖6。砂土、壤土的水平濕潤鋒入滲速度（VH）和垂直濕潤鋒入滲速度（VV）與時間的關系用下式描述：

式中：V為濕潤鋒入滲速度（mm/min）；t為時間（min）；a、b為常數，擬合關系見表4和表5。砂土和壤土的水平、垂直濕潤鋒入滲速度與時間均呈良好的冪函數關系[3-7]，相關系數R2都大于0.99。

H=0 m，曲線上同一時間點對應的濕潤鋒入滲速度砂土>壤土；H=-0.5 m和-1.0 m，曲線上同一時間點對應的濕潤鋒入滲速度壤土>砂土。

表4 水平濕潤鋒入滲速度與時間擬合值Table 4 Fitting parameters of infiltration velocity of Horizontal wetting front and time

表5 垂直濕潤鋒入滲速度與時間擬合值Table 5 Fitting parameters of infiltration velocity of vertical wetting front and time

3 模型與驗證

應用 Hydrus-3D建立三維土壤水分運移數值模型[27-29]，通過試驗數據驗證模型的準確性和有效性。

試驗扇形柱體高40.0 cm，徑向長度35.0 cm，夾角300。模型在z方向離散為5 mm，剖分為20層，共生成節點數111 383個，三角網格298 711個。計算時間步長為1 s。

3.1 參數選取

根據試驗結果，采用Van Genuchten模型對砂土、壤土的飽和含水率θs、殘余含水率θr、進氣值倒數α和擬合參數n進行求解，同時計算飽和滲透系數。參數結果見表6。

表6 土壤水力參數Table 6 Soil hydraulic parameters

3.2 模型計算結果

將累積入滲量與時間的冪函數關系作為流量邊界輸入模型，其他邊界處理為0流量邊界。按照物理模型試驗數據對應的時間點輸出計算結果，限于篇幅，3種負水頭條件下只給出了砂土、壤土最后一個時間點的濕潤鋒分布。

1）對于砂土，H=0 m時，30 min灌水器周圍達到飽和含水率0.47 cm3/cm3，最大的水平、垂直濕潤鋒距離為290、366 cm；H=-0.5 m時，3 480 min含水率降低至0.202 cm3/cm3，最大的水平、垂直濕潤鋒距離為278、314 cm；H=-1.0 m時，4 200 min含水率降低至0.129 cm3/cm3，最大的水平、垂直濕潤鋒距離為154、182 cm。見圖7。

圖7 砂土水分濕潤鋒運移計算結果Fig.7 Calculation results of wetting front moving in sand

2）對于壤土，H=0 m時，240 min灌水器周圍達到飽和含水率0.503 cm3/cm3，最大的水平、垂直濕潤鋒距離為289、330 cm；H=-0.5 m時，2 820 min含水率降低至0.487 cm3/cm3，最大的水平、垂直濕潤鋒距離為289、323 cm；H=-1.0 m時，4 200 min含水率降低至0.443 cm3/cm3，最大的水平、垂直濕潤鋒距離為261、292 cm。對基質勢較大的壤土，隨著負水頭高度的增加，灌水器周圍含水率只降低了11.9 %。結果見圖8。

圖8 壤土水分濕潤鋒運移計算結果Fig.8 Calculation results of wetting front moving in loam

3.3 模型驗證

3種負水頭條件，以負水頭高度0 m的實測濕潤鋒包絡面積作為模型率定數據，確定水分運移模式和參數等；然后，只改變流量邊界，即累積入滲量與時間的冪函數關系，以負水頭高度-0.5 m和-1.0 m的實測濕潤鋒包絡面積作為模型驗證數據。見圖9、圖10。

由表7可知3種負水頭模型，砂土計算濕潤鋒包絡面積與實測濕潤鋒包絡面積的偏差 0.51%～7.21%。由表8可知壤土計算濕潤鋒包絡面積與實測濕潤鋒包絡面積的偏差0.22%～16.03%。所建三維模型可以用于描述負水頭環境下土壤水分濕潤鋒的運移特征和規律。

表7 砂土模型驗證結果Table 7 Results of sand model validation

表8 壤土模型驗證結果Table 8 Results of loam model validation

圖9 砂土水分濕潤鋒運移模型驗證Fig.9 Model verification of wetting front moving in sand

圖10 壤土水分濕潤鋒運移模型驗證Fig.10 Model verification of wetting front moving in loam

3.4 濕潤深度和濕潤體體積

負水頭灌溉的優勢是可以控制土壤含水率，維持非飽和狀態，通過率定和驗證后的模型，可以得到不同負水頭、不同灌水時長條件下，各含水率的濕潤深度及濕潤體體積的動態變化（表9、表10）。

1）砂土負水頭高度H=0 m時，6～30 min灌水器周圍的濕潤深度由20.5 cm增加至36.7 cm，濕潤體體積由1 227.95 cm3增加至5 909.18 cm3。含水率0.376～0.344的濕潤深度由18.0 cm增加至30.6 cm，濕潤體體積由763.02 cm3增加至2 829.91 cm3，占總體積的62%減少至48%。

負水頭高度H=-0.5 m時，3 480 min后濕潤深度為31.5 cm，濕潤體體積為4 478.77 cm3，最大含水率為0.162。負水頭高度H=-1.0 m時，4 200 min后濕潤深度為31.5 cm，濕潤體體積為820.68 cm3，最大含水率為0.103。

2）壤土負水頭高度H=0 m時，10～240 min灌水器周圍的濕潤深度由16.5 cm增加至33.0 cm，濕潤體體積由345.51 cm3增加至5114.94 cm3。含水率0.402～0.390的濕潤深度由15.05 cm增加至30.6 cm，濕潤體體積由189.07 cm3增加至2 262.82 cm3，由占總體積的55%減少至44%。

負水頭高度H=-0.5 m時，50～2 820 min灌水器周圍的濕潤深度由21.0 cm增加至32.8 cm，濕潤體體積由484.66 cm3增加至5 015.60 cm3。含水率0.354～0.218的濕潤深度由15.45 cm增加至29.95 cm，濕潤體體積由256.09 cm3增加至3 521.29 cm3，由占總體積的53%增加至70%。

負水頭高度H=-1.0 m時，240～4 200 min灌水器周圍的濕潤深度由17.5 cm增加至28.75 cm，濕潤體體積由654.44 cm3增加至3 572.29 cm3。含水率0.320～0.183的濕潤深度由16.25 cm增加至26.35 cm，濕潤體體積由510.16 cm3增加至2 681.76 cm3，由占總體積的75%～78%。

表9 砂土濕潤體體積隨時間的動態變化情況Table 9 Dynamic change of wetting volume of sand with time cm3

表10 壤土濕潤體體積隨時間的動態變化情況Table 10 Dynamic change of wetting volume of loam with time cm3

由試驗和模擬結果可知，砂土負水頭高度H=0 m時，濕潤體48%～62%接近飽和含水率；H=-0.5 m時，濕潤體的最大含水率為16.2%。砂土的負水頭灌溉高度應介于-0.5~0 m之間較為合理。

壤土負水頭高度H=-0.5 m時，2 820 min后含水率0.286～0.218的濕潤深度為29.95 cm，濕潤體的體積為1 891.96 cm3，占總體積的38%；H=-1.0 m時，2 760 min后含水率0.286～0.218的濕潤深度為23.5 cm，濕潤體體積為1 434.91 cm3，占總體積的40%。壤土的負水頭灌溉高度可以介于-1.0～-0.5 m之間。

4 結 論

1）砂土、壤土的累積入滲量隨著時間的增加而逐漸增大，呈良好的冪函數關系。H=0 m，砂土的累積入滲量曲線高于壤土；H=-0.5 m或者-1.0 m，壤土的累積入滲量曲線高于砂土。

2）砂土、壤土的濕潤鋒隨著時間的增加逐漸向水平、垂直方向擴大。水平、垂直最大濕潤距離與時間的平方根呈良好的線性關系。H=0 m，單位時間內砂土濕潤鋒遷移距離大于壤土；H=-0.5 m或者-1.0 m，單位時間內壤土濕潤鋒遷移距離大于砂土。

3）砂土、壤土的濕潤鋒入滲速度隨著負水頭高度的增加逐漸減小，與時間呈良好的冪函數關系。H=0 m，曲線上同一時間點對應的濕潤鋒入滲速度砂土>壤土。H=-0.5 m或者-1.0 m，曲線上同一時間點對應的濕潤鋒入滲速度壤土>砂土。

4）應用Hydrus-3D軟件，建立三維土壤水分運移數值模型。率定和驗證后的模型計算濕潤鋒包絡面積與實測的偏差，砂土為 0.51%～7.21%，壤土為0.22%～16.03%，所建三維模型可以計算不同負水頭、不同灌水時長條件下，各含水率的濕潤深度及濕潤體體積的動態變化。