基于水平吸滲確定土壤水分擴散率的解析-修正法

帥嘉偉，龍 平,※，胡世麗，王觀石,，楊耀杰，羅嗣海

（1. 江西理工大學土木與測繪工程學院，贛州 341000；2. 江西理工大學江西省環(huán)境巖土與工程災害控制重點實驗室，贛州 341000）

0 引 言

研究土壤水分運動對農(nóng)業(yè)灌溉、水利工程和涉及水分入滲等的實際工程問題具有重要意義。土壤水分擴散率是分析土壤水分運動的關鍵參數(shù)之一，準確測試土壤水分擴散率是高精度模擬土壤水分運動的前提，因而，建立準確測試土壤水分擴散率的方法具有重要意義。

現(xiàn)有測試土壤水分擴散率的方法主要分為兩種。第一種是反分析法，Simu??nek等在Richards方程的基礎上，提出了估算水力參數(shù)的數(shù)值反分析法，進而確定土壤水分擴散率。大量研究表明，數(shù)值反分析法估算水力參數(shù)容易出現(xiàn)適定性問題（由于試驗誤差，使得數(shù)值反分析法存在多解性的問題，即適定性），且需要事先假設水力性質(zhì)的數(shù)學模型，容易引入模型誤差。第二種是水平吸滲法，Bruce等結合一維水平運動的Richards方程和Boltzmann變換，提出了測試土壤水分擴散率的水平吸滲法。水平吸滲法由于試驗簡單、耗時短而被廣泛使用。鄒小陽等采用水平吸滲法分析了殘膜對土壤水分擴散率的影響。姚淑霞等運用水平吸滲法研究了科爾沁沙地不同生境條件下的土壤水分擴散率。然而測試土柱瞬時剖面的含水率比較繁瑣，容易出現(xiàn)較大誤差。針對難以準確測試土柱剖面上含水率的問題，Whisler等提出在土柱上布置含水率測試儀，測試所布置位置處含水率隨時間的變化，再采用水平吸滲法分析含水率隨時間變化的數(shù)據(jù)，雖然能簡單、準確地獲得試驗數(shù)據(jù)，但對于高含水率，土壤水分擴散率易受邊界效應影響而出現(xiàn)較大誤差。

為了解決改進算法易受邊界效應影響的問題，本文分析一維水平入滲Richards的解析解，在此基礎上，對受邊界效應影響的數(shù)據(jù)進行修正，再采用水平吸滲法確定土壤水分擴散率，綜合水平吸滲法的優(yōu)勢，提出推導土壤水分擴散率的解析-修正法，并通過數(shù)值模擬和室內(nèi)試驗驗證方法的合理性，以期為農(nóng)業(yè)灌溉等土壤水運動的模擬提供精準參數(shù)。

1 解析-修正法的提出

1.1 一維水平吸滲Richards方程

第一類邊界條件下，一維無限長水平吸滲定解問題可以采用一維水平吸滲Richards方程式（1）和式（2）描述。

式中為時間，d；表示以土柱最左端為原點、水平向右為正方向的水平坐標，m；為土柱在= 0處的含水率，m/m；為初始含水率，m/m；為土壤水分擴散率，m/d；為體積含水率，m/m。

1.2 Richards方程求解

式中erfc為互補誤差函數(shù)；和為待定參數(shù)，為無量綱量；為互補誤差函數(shù)中自變量值，為無量綱量；為Boltzmann變換參數(shù)，=/，m/d。

采用Boltzmann變換處理式（1），可以得到式（8）。

將式（8）寫成差分形式，結果如式（9）所示，結合=/，基于～或～原始數(shù)據(jù)，采用式（9）確定土壤水分擴散率的方法稱為水平吸滲法。

式中為離散節(jié)點編號（= 1, 2, 3, …,-1，+1為總結點數(shù)），為無量綱量；λ為離散節(jié)點處的Boltzmann變換參數(shù)（λ=+(-1)?，?為的步長），m/d；D為離散節(jié)點處的土壤水分擴散率，m/d；θ為離散節(jié)點處的體積含水率，m/m；為離散節(jié)點編號（= 1, 2, 3, …,），為無量綱量；θ為離散節(jié)點處的體積含水率，m/m；λ為離散節(jié)點處的Boltzmann變換參數(shù)（λ=+(-1)?），m/d。

根據(jù)=/，將測點處含水率隨時間變化的數(shù)據(jù)（～數(shù)據(jù)）轉(zhuǎn)化為～數(shù)據(jù)，～數(shù)據(jù)分為兩部分，分別為受邊界和不受邊界影響的數(shù)據(jù)，如圖1所示。計算～曲線的斜率，在未受邊界影響的區(qū)域，斜率d/d隨著的減小而增加，當含水率受到邊界效應影響時，斜率d/d快速減小，隨著的減小，d/d由增加轉(zhuǎn)變?yōu)闇p小的轉(zhuǎn)折點對應的記為，根據(jù)～曲線確定對應的含水率，即為受邊界和不受邊界影響的界限含水率。先采用式（9）分析未受邊界影響的數(shù)據(jù)，得到[,]區(qū)間對應的土壤水分擴散率，然后采用式（7）擬合得到的～數(shù)據(jù)，確定參數(shù)和。將式（7）所示的解析解寫成差分形式，結果如式（10）所示。將式（10）代入式（9）中，整理可得式（11），采用式（11）修正[,]區(qū)間對應的含水率，再采用式（9）計算修正部分含水率對應的土壤水分擴散率，這種計算土壤水分擴散率的方法稱為解析-修正法，其具體流程如圖2所示。

圖1 受邊界影響和不受邊界影響數(shù)據(jù)的示意圖 Fig.1 Schematic diagram of data affected and unaffected by boundaries

圖2 解析-修正法流程圖 Fig.2 Flow chart of analytical-modified method

1.3 Hydrus-1D模擬含水率數(shù)據(jù)

采用Hydrus-1D軟件模擬水分在壤砂土（Loamy Sand，LS）、砂壤土（Sandy Loam，SL）、砂質(zhì)黏壤土（Sandy Clay Loam，SCL）、壤土（Loam，LO）、粉砂黏土（Silt Loam，SIL）和粉砂（Silt，SI）六種一維水平土柱中的入滲過程。水平土柱的長度設置為0.75 m，土柱的水力性質(zhì)采用van Genuchten 模型（簡稱VG模型）描述，VG模型所描述土壤水分擴散率如式（12）所示，六種土柱的VG模型參數(shù)的取值如表1所示。左邊界（= 0）設置為定含水率邊界（=，為土柱在= 0處的含水率，m/m），右邊界（=，為土柱長度，m）設置為自由排水邊界，初始含水率均設置為0.150 m/m。設置時間節(jié)點數(shù)量為50，記錄不同時刻下，體積含水率隨水平坐標的變化情況（～數(shù)據(jù)）。在= 0.30 m、= 0.45 m和= 0.60 m處布置3個體積含水率測點，監(jiān)測測點處體積含水率隨時間的變化（～數(shù)據(jù)）。

表1 基本參數(shù)[24] Table 1 Basic parameters

式中為土壤水分擴散率，m/d；=/[(-)]，和分別為飽和含水率和殘余含水率（m/m），、和為VG模型參數(shù)（= 1-1/），為飽和滲透系數(shù)（m/d）；為標準化含水率，= (-)/(-)，m/m；為體積含水率，m/m。

2 案例分析

2.1 供試土壤

試驗選擇了三種質(zhì)地的土壤，分別取自江西省信豐縣、江西省尋烏縣和福建省屏南縣的離子型稀土礦的表土，深度為2.3 m，采用比重計法測得各土壤的機械組成，結果如表2所示。根據(jù)《美國制土壤質(zhì)地分類標準》，信豐（Xinfeng）、尋烏（Xunwu）和屏南（Pingnan）的土壤質(zhì)地分別為砂質(zhì)壤土、粉砂質(zhì)壤土和壤土。

表2 土壤的機械組成 Table 2 Mechanical composition of soil

2.2 室內(nèi)試驗

取10根內(nèi)徑為0.14 m、長度為0.30 m的有機玻璃管（一端為開口端、一端為閉口端）進行水分傳感器標定試驗，距開口端0.15 m處開設0.04 m×0.02 m的矩形孔，并貼上硅膠墊。土柱質(zhì)量含水率分別設置為10%、11%、12%、13%、15%、16%、17%、18%、19%和20%，往土樣中分別加入對應質(zhì)量去離子水，充分攪拌均勻后裝入密封箱靜置48 h。設置土柱孔隙比為1.0，土壤容重為1.35×10N/m，分5層裝入土樣，將水分傳感器垂直插入硅膠墊處，設置采樣時間間隔為5 min，待傳感器讀數(shù)穩(wěn)定后停止采集數(shù)據(jù)，取最后5個數(shù)據(jù)的平均值作為傳感器最終讀數(shù)，采用烘干法測試探針周圍土樣的含水率。將質(zhì)量含水率換算為體積含水率（質(zhì)量含水率與體積含水率的換算關系為=/γ，為土壤容重，N/m；γ為水容重，N/m），即可得到水分傳感器讀數(shù)與體積含水率的關系數(shù)據(jù)，采用合適的數(shù)學表達式擬合試驗數(shù)據(jù)，確定水分傳感器的標定函數(shù)。

以信豐土壤為例，闡明室內(nèi)一維水平吸滲試驗過程如圖3所示。取內(nèi)徑0.14 m、長度0.80 m的有機玻璃管（一端為開口端，一端為閉口端，兩端裝有法蘭盤），有機玻璃管閉口端部開設一定數(shù)目直徑約為2 mm的小孔，在距離有機玻璃管開口端0.45 m和0.60 m處分別開設一個0.04 m×0.02 m的矩形孔，在矩形孔處的內(nèi)壁貼上硅膠墊。將土壤配制成質(zhì)量含水率約為7.0%的濕土壤，采用烘干法測試土壤準確的初始含水率，設置土柱柱長為0.73 m、孔隙比為1.0。首先往管內(nèi)裝入0.02 m厚度的粗砂，然后分5層裝入土樣，第一層為0.13 m，其他層為0.15 m，裝入下一層前刮毛土層上表面，待填樣完成后，裝入0.04 m厚度的粗砂，再放置一塊直徑為0.13 m、厚度為0.01 m且?guī)в幸欢〝?shù)目小孔的有機玻璃擋板，開口端連接直角彎管。在與土柱的上端平齊的直角彎管上開設一個溢流孔，在溢流孔處裝上一根溢流管。在設定位置處分別插入一個FDS-100水分傳感器（邯鄲開發(fā)區(qū)清易電子科技有限公司，河北邯鄲，測試精度：±3%），傳感器從左往右依次記為和，設置水分傳感器采樣時間間隔為5 min。

圖3 室內(nèi)一維水平吸滲試驗示意圖 Fig.3 Schematic diagram of indoor one-dimensional horizontal imbibition test

直角彎管內(nèi)快速加水至溢流孔處，之后采用蠕動泵全速注水。待兩個傳感器的讀數(shù)穩(wěn)定后試驗結束。結合水分傳感器標定函數(shù)，可以得到兩個測點處體積含水率隨時間的變化情況。取樣測定試驗結束時土柱的含水率。尋烏和屏南土壤的一維水平吸滲試驗與信豐的試驗步驟相同。

2.3 數(shù)據(jù)分析與統(tǒng)計方法

3 結果與分析

3.1 基于Hydrus-1D模擬的含水率的水分擴散率計算

采用Hydrus-1D軟件模擬水分在LS、SL、SCL、LO、SIL和SI六種土柱中的入滲過程，得到不同時刻的體積含水率在剖面的分布情況如圖4所示。以和為待定參數(shù)，采用式（6）擬合圖4所示的～數(shù)據(jù)，結果如表 3所示。對于同一種土，三個時刻的參數(shù)的擬合結果的離散程度較小，其中SI土的參數(shù)擬合結果的離散程度最大，變異系數(shù)也僅為2.31%，由此可見，對于同一種土，參數(shù)和變化較??；對于不同質(zhì)地的土，參數(shù)和取值具有明顯差異，參數(shù)和的取值與土壤質(zhì)地相關。六種土三個時刻的決定系數(shù)均大于0.990。由此可見，式 （4）或式（6）可以較好地描述第一類邊界條件下，水分在半無限長土柱中入滲的含水率隨時間和坐標的變化情況，即在第一類邊界條件下，式（4）或式（6）可以作為式（1）的近似解析解。

表3 參數(shù)擬合結果 Table 3 Parameter fitting results

圖4 體積含水率隨水平土柱上x軸坐標長度的變化 Fig.4 Variation of volumetric water content with the length ofx-axis coordinate on horizontal soil column

= 0.30 m、= 0.45 m和= 0.60 m處體積含水率隨時間變化的模擬結果如圖5所示。利用方程=/將圖 5的～數(shù)據(jù)轉(zhuǎn)化為～數(shù)據(jù)，再對～數(shù)據(jù)進行等距節(jié)點插值（節(jié)點數(shù)量為30），結果如圖6所示。采用解析-修正法分析～數(shù)據(jù)，得到～數(shù)據(jù)的修正結果同樣繪制于圖6中，土壤水分擴散率的計算結果如圖7所示。對于同一種土，半無限長土柱內(nèi)的～曲線是唯一的，即相同的，對應唯一的，不隨時間或坐標的影響，當時間較小時，水分還未入滲至土柱最右端，其含水率隨坐標的變化情況與半無限長土柱內(nèi)的相同，因此可以采用較小時～數(shù)據(jù)確定唯一的～曲線（本文取圖4中第一個時刻數(shù)據(jù)，該該時刻記為）。由圖6可知，修正后的～數(shù)據(jù)與處的～數(shù)據(jù)基本完全重合，說明修正方法是合理的。越靠近入水端的測點，得到的修正后的含水率范圍越大（如對于LS土柱，= 0.30 m和= 0.60 m修正后的最大含水率分別為0.373和0.346 m/m），使得土壤水分擴散率計算結果的范圍也更廣（如圖7所示）。這是由于對于同一種土，3個測點的總時長相同，越靠近入水端，得到的越小，與呈負相關關系，因而含水率范圍越大。

圖5 體積含水率隨時間的變化情況 Fig.5 Change of volumetric water content with time

采用水平吸滲法分析圖6中的數(shù)值模擬～數(shù)據(jù)，得到土壤水分擴散率的結果如圖7所示。由圖7可知，采用水平吸滲法分析受邊界效應影響的數(shù)據(jù)，得到擴散率的計算結果明顯偏離真實值（擴散率的真實值為Hydrus-1D模擬時所采用的土壤水分擴散率，即把表1的參數(shù)代入式（12）得到的結果），出現(xiàn)較大的誤差。由此可見，采用水平吸滲法分析～數(shù)據(jù)得到的～數(shù)據(jù)具有較大的局限性，僅能較為準確地得到不受邊界效應部分含水率對應的土壤水分擴散率。

圖6 體積含水率隨λ參數(shù)的變化情況 Fig.6 Variation of volumetric water content withλ parameter

解析-修正法的計算結果也繪制于圖7中。由圖7可知，相比水平吸滲法，在受邊界效應影響的含水率區(qū)域，采用解析-修正法計算得到的土壤水分擴散率也基本與1∶1線重合，計算結果精度較高；整個含水率區(qū)域，采用解析-修正法得到的土壤水分擴散率與真實值的決定系數(shù)均在0.900以上，而水平吸滲法計算結果的決定系數(shù)基本在0.600以下，甚至出現(xiàn)負相關的情況。由此可見，解析-修正法能利用～數(shù)據(jù)準確確定土壤水分擴散率。

圖7 土壤水分擴散率的解析-修正法和水平吸滲法的計算結果 Fig.7 Results of the analytical-method and horizontal imbibition method of soil water diffusivity

3.2 解析-修正法驗證結果

基于室內(nèi)試驗，對兩個水分傳感器進行標定，如圖8所示，線性回歸方程擬合精度較高（>0.97），能夠用于土壤含水率計算?；趫D8中方程得到兩個測點的含水率動態(tài)變化數(shù)據(jù)（～）如圖9所示。利用方程=/將～數(shù)據(jù)轉(zhuǎn)化為～數(shù)據(jù)（等距差值數(shù)量為15），結果如圖10所示，根據(jù)～數(shù)據(jù)可以確定和，結果也列于圖10中，采用式（9）分析圖10中[,]區(qū)間的數(shù)據(jù)，確定相應的土壤水分擴散率D，再采用式（10）擬合未受邊界效應影響的～數(shù)據(jù)，得到近似解析解參數(shù)和的值如表4所示。再采用式（11）修正圖10中[,]區(qū)間對應的含水率，修正后的～曲線繪于圖10中，采用式（9）分析修正后的～曲線，得到土壤水分擴散率如圖11所示，水平吸滲法的結果也繪于圖11中。由圖11可知，當>時，水平吸滲法得到的土壤水分擴散率隨著含水率的增加而減小，土壤水分擴散率隨含水率變化的趨勢不相符，說明不能直接采用水平吸滲法計算受邊界效應影響的土壤水分擴散率。

圖8 土壤水分傳感器標定曲線 Fig.8 Calibration curve of soil moisture sensor

表4 近似解析解參數(shù)的擬合結果 Table 4 Fitting results of approximate analytic solution parameters

圖9 室內(nèi)一維水平土柱吸滲試驗θ～t數(shù)據(jù) Fig.9θ-t data of one-dimensional horizontal soil column imbibition test

圖10 體積含水率θ隨Boltzmann參數(shù)λ變化的試驗結果與修正結果 Fig.10 Experimental and modified results of volumetric water contentθ varying with Boltzmann parameterλ

圖11 解析-修正法和水平吸滲法的土壤水分擴散率計算結果 Fig.11 Results of soil water diffusivity calculated by analytic - modified method and horizontal imbibition method

對于解析-修正法的結果，當>時，土壤水分擴散率隨著含水率的增加而快速增加，滿足土壤水分擴散率隨含水率變化的基本趨勢；對于同一土柱，根據(jù)兩個測點計算得到的土壤水分擴散率吻合度很好，= 0.60 m的含水率范圍，土壤水分擴散率計算結果的決定系數(shù)分別為0.958、0.962和0.951，說明不同測點的計算結果具有很好的一致性。

4 結 論

1）結合線性化水平入滲Richards方程的解析解和常數(shù)變易法，推導了一維水平入滲Richards方程的近似解析解，與Hydrus-1D模擬結果的決定系數(shù)均在0.990以上，驗證了近似解析解是合理的。近似解析解中的參數(shù)與土壤性質(zhì)相關，對于同一種土，參數(shù)變化較小。

2）先采用近似解析解修正受邊界效應影響的數(shù)據(jù)，再采用水平吸滲法確定土壤水分擴散率，建立了解析-修正法，結果發(fā)現(xiàn)，解析-修正法能精確地修正受邊界效應影響的數(shù)據(jù)，且越靠近入水端，修正的含水率范圍越大。

3）室內(nèi)試驗和Hydrus-1D模擬試驗都表明，含水率隨時間變化的數(shù)據(jù)易受邊界效應的影響，直接采用水平吸滲法確定土壤水分擴散率，在高含水率區(qū)的精度較差，整體決定系數(shù)均在0.600以下，而采用解析-修正法得到的土壤水分擴散率的決定系數(shù)均在0.900以上，彌補了水平吸滲法的不足。

4）本文方法沒有事先對試驗數(shù)據(jù)進行處理，減少測試誤差的影響，導致得到的土壤水分擴散率的計算易受測試誤差的影響。進一步研究中需引入試驗數(shù)據(jù)預處理算法（如移動平均法），提高測試結果的精度。