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

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

0 引 言

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

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

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

1 解析-修正法的提出

1.1 一維水平吸滲Richards方程

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

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

1.2 Richards方程求解

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

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

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

式中為離散節(jié)點編號（= 1, 2, 3, …,-1，+1為總結(jié)點數(shù)），為無量綱量；λ為離散節(jié)點處的Boltzmann變換參數(shù)（λ=+(-1)?，?為的步長），m/d；D為離散節(jié)點處的土壤水分?jǐn)U散率，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āng)含水率受到邊界效應(yīng)影響時，斜率d/d快速減小，隨著的減小，d/d由增加轉(zhuǎn)變?yōu)闇p小的轉(zhuǎn)折點對應(yīng)的記為，根據(jù)～曲線確定對應(yīng)的含水率，即為受邊界和不受邊界影響的界限含水率。先采用式（9）分析未受邊界影響的數(shù)據(jù)，得到[,]區(qū)間對應(yīng)的土壤水分?jǐn)U散率，然后采用式（7）擬合得到的～數(shù)據(jù)，確定參數(shù)和。將式（7）所示的解析解寫成差分形式，結(jié)果如式（10）所示。將式（10）代入式（9）中，整理可得式（11），采用式（11）修正[,]區(qū)間對應(yīng)的含水率，再采用式（9）計算修正部分含水率對應(yīng)的土壤水分?jǐn)U散率，這種計算土壤水分?jǐn)U散率的方法稱為解析-修正法，其具體流程如圖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）六種一維水平土柱中的入滲過程。水平土柱的長度設(shè)置為0.75 m，土柱的水力性質(zhì)采用van Genuchten 模型（簡稱VG模型）描述，VG模型所描述土壤水分?jǐn)U散率如式（12）所示，六種土柱的VG模型參數(shù)的取值如表1所示。左邊界（= 0）設(shè)置為定含水率邊界（=，為土柱在= 0處的含水率，m/m），右邊界（=，為土柱長度，m）設(shè)置為自由排水邊界，初始含水率均設(shè)置為0.150 m/m。設(shè)置時間節(jié)點數(shù)量為50，記錄不同時刻下，體積含水率隨水平坐標(biāo)的變化情況（～數(shù)據(jù)）。在= 0.30 m、= 0.45 m和= 0.60 m處布置3個體積含水率測點，監(jiān)測測點處體積含水率隨時間的變化（～數(shù)據(jù)）。

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

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

2 案例分析

2.1 供試土壤

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

表2 土壤的機(jī)械組成 Table 2 Mechanical composition of soil

2.2 室內(nèi)試驗

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

以信豐土壤為例，闡明室內(nèi)一維水平吸滲試驗過程如圖3所示。取內(nèi)徑0.14 m、長度0.80 m的有機(jī)玻璃管（一端為開口端，一端為閉口端，兩端裝有法蘭盤），有機(jī)玻璃管閉口端部開設(shè)一定數(shù)目直徑約為2 mm的小孔，在距離有機(jī)玻璃管開口端0.45 m和0.60 m處分別開設(shè)一個0.04 m×0.02 m的矩形孔，在矩形孔處的內(nèi)壁貼上硅膠墊。將土壤配制成質(zhì)量含水率約為7.0%的濕土壤，采用烘干法測試土壤準(zhǔn)確的初始含水率，設(shè)置土柱柱長為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ù)目小孔的有機(jī)玻璃擋板，開口端連接直角彎管。在與土柱的上端平齊的直角彎管上開設(shè)一個溢流孔，在溢流孔處裝上一根溢流管。在設(shè)定位置處分別插入一個FDS-100水分傳感器（邯鄲開發(fā)區(qū)清易電子科技有限公司，河北邯鄲，測試精度：±3%），傳感器從左往右依次記為和，設(shè)置水分傳感器采樣時間間隔為5 min。

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

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

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

3 結(jié)果與分析

3.1 基于Hydrus-1D模擬的含水率的水分?jǐn)U散率計算

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

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

圖4 體積含水率隨水平土柱上x軸坐標(biāo)長度的變化 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處體積含水率隨時間變化的模擬結(jié)果如圖5所示。利用方程=/將圖 5的～數(shù)據(jù)轉(zhuǎn)化為～數(shù)據(jù)，再對～數(shù)據(jù)進(jìn)行等距節(jié)點插值（節(jié)點數(shù)量為30），結(jié)果如圖6所示。采用解析-修正法分析～數(shù)據(jù)，得到～數(shù)據(jù)的修正結(jié)果同樣繪制于圖6中，土壤水分?jǐn)U散率的計算結(jié)果如圖7所示。對于同一種土，半無限長土柱內(nèi)的～曲線是唯一的，即相同的，對應(yīng)唯一的，不隨時間或坐標(biāo)的影響，當(dāng)時間較小時，水分還未入滲至土柱最右端，其含水率隨坐標(biāo)的變化情況與半無限長土柱內(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），使得土壤水分?jǐn)U散率計算結(jié)果的范圍也更廣（如圖7所示）。這是由于對于同一種土，3個測點的總時長相同，越靠近入水端，得到的越小，與呈負(fù)相關(guān)關(guān)系，因而含水率范圍越大。

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

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

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

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

圖7 土壤水分?jǐn)U散率的解析-修正法和水平吸滲法的計算結(jié)果 Fig.7 Results of the analytical-method and horizontal imbibition method of soil water diffusivity

3.2 解析-修正法驗證結(jié)果

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

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

表4 近似解析解參數(shù)的擬合結(jié)果 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ù)λ變化的試驗結(jié)果與修正結(jié)果 Fig.10 Experimental and modified results of volumetric water contentθ varying with Boltzmann parameterλ

圖11 解析-修正法和水平吸滲法的土壤水分?jǐn)U散率計算結(jié)果 Fig.11 Results of soil water diffusivity calculated by analytic - modified method and horizontal imbibition method

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

4 結(jié) 論

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

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

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

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