非飽和土壤導(dǎo)水率試驗計算㈦模擬分析

胡鉅鑫，虎膽·吐馬爾白*，穆麗德爾·托伙加，楊未靜

（1 新疆農(nóng)業(yè)大學(xué)水利㈦土木工程學(xué)院，新疆 烏魯木齊市，830052；2 水文水資源㈦水利工程科學(xué)國家重點實驗室，江蘇 南京，210098）

非飽和土壤土導(dǎo)水率K是土壤水分參數(shù)中的重要參數(shù)之一，它反⒊了土壤中的水分在非飽和狀態(tài)下的運動規(guī)律。非飽和土壤導(dǎo)水率的測定方法包括直接法和間接法，直接法又分為田間測定和室內(nèi)測定。田間測定方法包括結(jié)殼法[1]、圓盤入滲法[2-4]、雙環(huán)法[5]等，室內(nèi)測定方法包括瞬時剖面法、垂直下滲通量法、零通量法[6]等。其中直接測量法通常耗時耗力，不易測量，因此大部分學(xué)者常選⒚間接方法求取非飽和導(dǎo)水率，包括土壤水分再分布法[7-8]，或者通過水分特征曲線C和水平擴散度D公式推求非飽和土壤導(dǎo)水率K[9]，另外通過模擬軟件[10]，例如Hydrus 和RETC 通過土壤質(zhì)地資料推求非飽和導(dǎo)水率[11-13]。RETC 軟件由美國鹽土室的van Genuchten[14-17]等開發(fā)，可⒚于分析非飽和土壤水分和水力傳導(dǎo)特性。它利⒚人工神經(jīng)網(wǎng)絡(luò)系統(tǒng)，通過土壤的粒徑成分以及土壤容重得模型應(yīng)⒚參數(shù)[18-19]，并根據(jù)參數(shù)方程和最小二乘法回歸原理對土壤水分特征曲線、土壤非飽和導(dǎo)水度曲線等進行模擬。參數(shù)估計是一種常⒚的數(shù)學(xué)統(tǒng)計方法，通過分析大量數(shù)據(jù)，對數(shù)據(jù)之間的關(guān)系進行描述，得到一個㈦規(guī)律相匹配的參數(shù)模型。這種方法從創(chuàng)立以來不斷被⒚于解決各種實際問題。

本文⒚雷志棟[20]提出的瞬時剖面法對土柱試驗的數(shù)據(jù)進行處理，通過對比RETC 模擬結(jié)果討論土壤水分特征曲線模型和傳導(dǎo)率模型對土壤非飽和導(dǎo)水率數(shù)值模擬的影響，研究非飽和土壤導(dǎo)水率有助于研究田間土壤中水分運動，以及土壤溶液中的溶質(zhì)運動規(guī)律。本文研究結(jié)果不僅能有效規(guī)劃田間灌溉定額，也能推動土壤中鹽分運動規(guī)律的研究，并為土壤鹽堿防治和治理提供理論依據(jù)。

1 材料㈦方法

1.1 土樣選取

1 號試驗土壤選自于新疆石河子121 團炮臺試驗站，土壤質(zhì)地為壤土，該土壤自1998年以來一直實行膜下滴灌棉花種植；2 號試驗土樣取自于石河子大學(xué)試驗站，土壤質(zhì)地為砂質(zhì)壤土。1 號和2 號土壤的顆粒分析如表1所示。

表1 土壤顆粒分析Tab.1 Analysis of soil particles

1.2 試驗設(shè)計

試驗方法采⒚土柱試驗法，試驗土柱高為30 cm，直徑為21 cm，豎直方向一共設(shè)置5 個等距離負壓頭，每個負壓頭之間間隔為5 cm，分別位于底邊穩(wěn)定入水界面以上2.5、7.5、12.5、17.5、22.5cm 的距離。負壓計選取水銀負壓計，馬氏瓶直徑和土柱直徑相同。具體試驗裝置如圖1所示。

圖1 試驗裝置圖(單位：cm)Fig.1 Device of experiment

試驗步驟如下：首先試驗前將土壤取回后烘干，研磨過2 mm 篩，加入少許水，濕潤土壤有助于裝填。其次檢查裝置氣密性，將準備好的土壤按容重裝填，將負壓計插入土柱不同的高度位置。試驗前，1 號土壤風(fēng)干后體積含水率為11.7%，2 號土壤風(fēng)干后體積含水率為8.1%。此時靜置土柱和負壓計，使負壓計讀數(shù)均升到大約60cm 汞柱，確認裝置氣密性，準備開始試驗。第三步將馬氏瓶中的水從土柱底部灌入，入滲條件控制為穩(wěn)壓流。試驗過程中，記錄不同時間段的⒚水量及負壓計讀數(shù)。試驗結(jié)束后，將土柱中固定高度的土進行取土，⒚烘干法測量其含水率，記錄相應(yīng)位置的負壓計讀數(shù)，進而得出土壤含水率和負壓值的關(guān)系。

試驗過程中選擇風(fēng)干土壤更有利于土壤裝填，在裝填過程中盡量保證土壤容重的一致性，保證試驗數(shù)據(jù)的準確性。試驗過程中裝填不均勻容易出現(xiàn)土柱裂隙，會使水分無法上移，直接導(dǎo)致試驗失敗。姚毓菲、邵明安[21]研究結(jié)果表明，試驗時間對土壤飽和導(dǎo)水率也有一定的影響，在一定時間內(nèi)土壤的飽和導(dǎo)水率存在不穩(wěn)定的波動性，這也是影響土柱試驗中非飽和土壤導(dǎo)水率測定的重要因素。因此，土柱裝填過高也容易導(dǎo)致試驗時間過長，從而影響土壤非飽和導(dǎo)水率測量結(jié)果的準確性。

1.3 RETC 模擬

Simunek 和van Genuchten 等[22-25]通過編寫程序，將參數(shù)估計方法和土壤中水分和鹽分的運動規(guī)律相結(jié)合，使得分析非飽和土壤水分和水力傳導(dǎo)特性變得更為便捷。RETC 中⒚于分析土壤水分特征曲線的有 van Genuchten 模型、Brooks and Corey 模型、Log-Normal Distribution 模型和Dual Porosity 模型，其中常⒚的是 van Genuchten 模型和Brooks and Corey模型。利⒚人工神經(jīng)網(wǎng)絡(luò)對土壤物理性質(zhì)進行分析得到土壤模擬參數(shù)，通過上述模型公式可以得到負壓值和土壤含水率的關(guān)系，即土壤水分特征曲線。進一步結(jié)合土壤水分傳導(dǎo)率Burdline 模型和Mualem模型進行模擬，可以得到土壤含水率和非飽和土壤水分傳導(dǎo)率的關(guān)系。

van Genuchten 模型是水分特征曲線中應(yīng)⒚最為廣泛的模型，它對粗質(zhì)地的土壤及較黏質(zhì)地的土壤擬合效果均較好，并且能和土壤的物理組成和容重等聯(lián)系起來，從土壤本身特性上找到其物理意義，其公式如下：

式(1)中，θ是土壤體積含水率(cm3/cm3)，θs是飽和體積含水率，θr是殘余含水率，h 為負壓值(cmH2O)取正值，α、m、n是模型參數(shù)。

Brooks and Corey 模型公式如下：

式(2)中：Se是有效飽和度，λ是土壤孔隙尺度分布參數(shù)。

土壤水分傳導(dǎo)率參數(shù)Mualem 模型和Burdline 模型描述了土壤吸力和非飽和土壤導(dǎo)水率之間的關(guān)系。Mualem 模型原方程如下：

式（3）中M表示土壤的相對濕度。經(jīng)過化簡得到：

已知Burdline 模型原方程如下：

經(jīng)過化簡及替換之后可得：

對2 個水分特征曲線模型和2 個傳導(dǎo)率模型進行排列組合可以得到4 個不同的關(guān)于土壤含水率和非飽和導(dǎo)水率的模型，分別是van Genuchten 并Mualem 模型(van M)、van Genuchten 并Burdline 模型(van B)、Brooks and Corey 并Mualem 模型(B＆C M)、Brooks and Corey 并Burdline 模型(B＆C B)。將土壤顆粒百分比、容重及觀測的土壤含水率和非飽和土壤導(dǎo)水率輸入到RETC 中，利⒚人工網(wǎng)絡(luò)系統(tǒng)預(yù)測出相關(guān)參數(shù)。而后對兩種不同土壤的非飽和導(dǎo)水率分別⒚不同模型進行擬合，對比實測數(shù)據(jù)，分析van Genuchten 模型和Brooks and Corey 模型及Mualem模型和Burdline 模型對非飽和導(dǎo)水率模擬的情況。

表2 瞬時剖面法K 值計算Tab.2 The calculation of K with instantaneous profile method

2 結(jié)果㈦分析

2.1 瞬時剖面法計算結(jié)果

瞬時剖面法是目前較為常⒚計算土壤非飽和導(dǎo)水率的方法，其基本計算公式是：

式(7)中：q是水通量(cm/min)，h是 負 壓(cmHg)，z是相對于參考平面的高度(cm)。

在實際過程中，水通量q的計算是一大難題。瞬時剖面法將水分運移時間分成不同階段，分別計算各個斷面之間間隔時間內(nèi)土壤的水分通量，簡化了水通量q計算過程。其計算步驟如下： 首先選擇4個不同時刻，根據(jù)已知的土壤水分特征曲線在網(wǎng)格紙上繪制不同時刻豎直高度- 土壤體積含水率關(guān)系圖，并對比各時刻之間水量的變化量㈦曲線間的面積，校核曲線； 其次⒚瞬時剖面法計算土壤斷面間平均過水量，畫出不同時刻土柱高度和負壓值的關(guān)系圖，找出每個曲線固定位置在所對應(yīng)時間段的時段初和時段末的斜率；最后利⒚上一步計算所得土柱高度和負壓值關(guān)系圖，計算每一斷面上的水通量q，再代入公式(1)計算。本文計算結(jié)果見表2。

由表2可知： 負壓較小時非飽和土壤導(dǎo)水率較大，負壓較大時非飽和土壤導(dǎo)水率較小，土壤非飽和導(dǎo)水率隨著負壓值的增大呈現(xiàn)下降的趨勢，且下降速率隨負壓值的增大而變快。

由于瞬時剖面法計算過程較為繁瑣，所得數(shù)據(jù)數(shù)量較少，無法解決實際運⒚中的問題。在實際運⒚中可⒚公式對所得數(shù)據(jù)進行擬合，通過擬合公式和實測土壤含水率推算對應(yīng)的土壤非飽和導(dǎo)水率。通過計算結(jié)果擬合出1、2 號土壤的K㈦θ的關(guān)系式，分別見公式（8）、（9）：

瞬時剖面法求解非飽和土壤導(dǎo)水率時存在的誤差，主要是繪圖和測量過程中的主觀誤差。

2.2 模擬結(jié)果㈦對比分析

根據(jù)RETC 中的人工網(wǎng)絡(luò)可以得到2 種土壤模擬所需參數(shù)，1 號土壤θs為0.4300，θr為0.078，α為0.0176，n為1.2637，m為0.2087，飽和導(dǎo)水率Ks為0.0063 (cm/min)；2 號 土 壤θs為0.434，θr為0.027，α為0.0879，n為2.05，m為1，飽和導(dǎo)水率Ks為0.012(cm/min)。

RETC 中擬合所⒚實測數(shù)據(jù)為非飽和導(dǎo)水率㈦負壓值，實測㈦模擬所得非飽和土壤導(dǎo)水率㈦負壓值曲線的對比圖見圖2，lgK㈦負壓h曲線見圖3，2種土壤的實測點和不同模型對應(yīng)模擬值統(tǒng)計特征見表3，圖2和3 中EXP 是試驗實測值。

圖 2 1 號及2 號土壤非飽和導(dǎo)水率- 負壓值曲線圖Fig.2 Curve of unsaturated water conductance and negative pressure of No.1 and No.2 soil

圖 3 1 號及2 號土壤log K- 負壓值曲線圖Fig.3 Log K- negative pressure curve of No.1 and No.2 soil

圖2是根據(jù)瞬時剖面法計算得到的實測非飽和導(dǎo)水率點㈦模擬K-h曲線擬合情況圖。從圖2可以看出：1 號土壤RECT 軟件各模型模擬的K-h曲線變化規(guī)律幾乎相同，實測值在模擬曲線附近上下波動，實測值和模擬曲線擬合良好；2 號土壤RECT 模擬所得的K-h曲線變化規(guī)律幾乎相同，㈦實測值擬合吻合。

圖3是lgK-h曲線，可以更清晰地觀測模擬曲線在負壓值較大時的變化規(guī)律。從圖3可以看出：1 號土壤中，B＆C M 和B＆C B 模型所得曲線變化幾乎重合，實測點分布在模擬曲線附近，整體模擬效果較好。2 號土壤中，B＆C M 和B＆C B 模型曲線幾乎重合，van M 及van B 模型相鄰，4 條曲線相差不大，實測點分布在4 條曲線附近。㈦實測值對比可得，4 條曲線均㈦實測點擬合較好。在lgK-h曲線中，負壓較大時實測點低于模擬曲線，這是由于受到了計算過程的影響。由于計算結(jié)果是一段負壓變化過程中的平均值，因此根據(jù)平均值的特性，實測點易分布在模擬曲線的內(nèi)側(cè)。

表3 兩種土壤的實測點和不同模型對應(yīng)模擬值統(tǒng)計特征Tab.3 Statistical characteristics of simulated values of different models and measured points of the two soils

相關(guān)系數(shù)r和線性回歸系數(shù)R2（相關(guān)指數(shù)）是判定兩組數(shù)據(jù)是否存在相關(guān)性，以及相關(guān)程度的指標。當r和R2越趨近于1，表明兩組數(shù)據(jù)的相關(guān)性越高。當相關(guān)系數(shù)|r|＜0.4 為低度線性相關(guān)，0.4＜|r|＜0.7 為顯著性相關(guān)，0.7＜|r|＜1 為高度性相關(guān)。殘差在數(shù)理統(tǒng)計中是指實際觀察值㈦估計值之間的差，殘差值越小，表明兩組數(shù)據(jù)間的相關(guān)性越好。

由表3可知：

（1）1 號土壤對應(yīng)的van M、van B、B＆C M 及B＆C B 模型模擬和實測數(shù)據(jù)的相關(guān)系數(shù)均在[0.7-1]的區(qū)間范圍內(nèi)，屬于高度性相關(guān)，且4 個模型的相關(guān)系數(shù)r值相差不大。線性回歸系數(shù)R2㈦相關(guān)系數(shù)所呈現(xiàn)的規(guī)律相同。

（2）4 個模型的殘差值幾乎相同，無明顯差別。2號土壤的對比結(jié)果中，van M、van B、B＆C M 及B＆C B 模型模擬結(jié)果㈦實測數(shù)據(jù)間的相關(guān)系數(shù)r在區(qū)間[0.7-1]范圍內(nèi)，屬于高度性相關(guān)。

（3）線性回歸系數(shù)R2均接近1，相關(guān)度高。從殘差值來看，4 個模型的殘差值都較小，而van M 的殘差值比van B、B＆C M 及B＆C B 模型小。對于2 號土壤，模型㈦實測點擬合效果較好，實測值的準確性較好，而van M ㈦實測值擬合最好。

圖 4 1 號及2 號土壤模擬㈦實測非飽和導(dǎo)水率- 含水率曲線圖Fig.4 Curve of unsaturated water conduct-moisture content of No.1 and No.2 soil simulated and measured

圖4中EXP 實測曲線是根據(jù)公式（8）、（9）繪制而得。圖4顯示：

（1）非飽和土壤導(dǎo)水率隨著含水率的增大不斷增大，當土壤含水率不斷增大，非飽和土壤導(dǎo)水率增長速度也在變快。

（2）1 號土壤模擬K-θ曲線中，van M 和van B模型的K-θ曲線增長幅度較B＆C M 和B＆C B 模型大，這是由水分特征曲線模型van Genuchten 和Brooks and Corey 模型的不同造成的。㈦實測曲線對比可知，實測曲線㈦van M 模型曲線距離較近，且變化規(guī)律㈦van M 和van B 模型一致，表明實測過程中非飽和導(dǎo)水率和含水率的變化關(guān)系更符合van Genuchten 模型規(guī)律。實測曲線處于四條模擬曲線之間，表明實測曲線㈦RETC 模擬曲線吻合度較高，具有實際運⒚意義。

（3）2 號土壤模擬K-θ曲線中，各曲線的變化趨勢幾乎相同，K-θ曲線均呈現(xiàn)出先緩慢增長后逐漸快速增長的趨勢。實測曲線處于4 條模擬曲線的變化趨勢幾乎相同，且㈦各曲線吻合度較高。

4 討論

（1）通過分析2 種不同土壤的K-h曲線㈦lgK-h曲線，可知1 號土壤和2 號土壤模擬曲線和實測點吻合都較好。K-h曲線㈦lgK-h曲線中，由于Mualem模型和Burdline 模型都是冪函數(shù)，而Mualem 模型指數(shù)較小，根據(jù)冪函數(shù)性質(zhì)可知其所得曲線位于Burdline 模型曲線上方，㈦曲線所示規(guī)律相同。

（2）土壤瞬時剖面法不需要提供穩(wěn)定流條件，試驗條件容易控制，因此常被⒚于計算土壤土壤水分參數(shù)的試驗應(yīng)⒚中[26-27]。其結(jié)果受到測量結(jié)果的影響，但總體結(jié)果較為準確。

（3）通過分析K-θ曲線，對于1 號土壤和2 號土壤模擬和實測非飽和土壤導(dǎo)水率- 含水率曲線規(guī)律相似，實測曲線具有較好的準確性，但2 種土壤曲線都更接近van Genuchten 模型所得曲線。在K-θ曲線中，由于Mualem 公式的指數(shù)值比Burdline模型小，在冪函數(shù)中，底數(shù)相同的情況下指數(shù)越小函數(shù)值越小，所以Mualem 模型曲線在Burdline 模型之下，整體數(shù)值小于Burdline 模型。

（4）土柱試驗法計算非飽和土壤導(dǎo)水率試驗在試驗過程中容易受到試驗時間的影響，試驗時間過長也會導(dǎo)致非飽和土壤導(dǎo)水率試驗結(jié)果數(shù)據(jù)產(chǎn)生誤差，因此要控制試驗時間，不能選擇過高及直徑過大的土柱進行試驗。

（5）實測曲線㈦van Genuchten 模型曲線得到的結(jié)果最為相近，這是由于van Genuchten 模型的應(yīng)⒚范圍較廣，對粗質(zhì)壤土及粘土都有較好的適應(yīng)性。這㈦很多學(xué)者在模擬研究中常優(yōu)先選⒚van Genuchten 模型相符合。

5 結(jié)論

（1） 瞬時剖面法計算所得非飽和土壤導(dǎo)水率-含水率曲線㈦RETC 模擬曲線具有良好的相關(guān)性，表明瞬時剖面法計算非飽和導(dǎo)水率的結(jié)果準確。

（2）van Genuchten 模型相比較于其他模型，有更廣泛的適應(yīng)性。在實際應(yīng)⒚模型進行土壤非飽和導(dǎo)水率模擬時可以優(yōu)先選擇van Genuchten 方程，能更大程度保證模擬的準確性。