康守旋 ，費良軍 ※，鐘 韻 ，趙彭輝 ，楊 震 ，樊倩雯

（1. 西安理工大學省部共建西北旱區(qū)生態(tài)水利國家重點實驗室,西安 710048；2. 三峽大學三峽庫區(qū)生態(tài)環(huán)境教育部工程研究中心,宜昌 443002）

0 引言

中國黃河流域大部分區(qū)域降水稀少[1-2]，農業(yè)灌溉用水嚴重不足[3-4]，制約了當?shù)剞r業(yè)發(fā)展。該區(qū)域降雨分布集中，水土流失嚴重[5]，使得黃河成為世界上含沙率最高的河流[6]。為充分利用水資源來滿足農業(yè)灌溉用水的需求，該地區(qū)的一些灌區(qū)大力開展引渾灌溉。渾水灌溉可以在提高水資源利用率的同時改良土壤、增加土壤肥力。進入農田的渾水入滲時，部分泥沙顆粒隨水分運移進入土壤中，滯留在土層內，形成滯留層，而渾水中大部分泥沙顆粒逐漸沉積在土壤表面，形成沉積層。滯留層和沉積層共同形成了結構不同于原土的致密層[7-9]，從而改變了上邊界入滲條件，對入滲過程產生影響。因此，深入研究渾水入滲過程對緩解引渾灌區(qū)水資源緊缺、擴灌增產和提高灌溉質量等具有重要的理論價值和現(xiàn)實意義。

渾水灌溉是引用含沙河水作為水源進行灌溉。渾水在進入田間后，其所含的泥沙隨水流推進、沉降，在土體表面產生致密層，形成了不同于清水的土壤水分入滲規(guī)律[10-11]。一些研究采用渾水試驗方法探索渾水入滲規(guī)律及影響因素[12-13]，發(fā)現(xiàn)：累積入滲量與入滲時間呈冪函數(shù)變化關系[12]；渾水含沙率、泥沙顆粒級配、土壤容重及土壤初始含水率等均能影響渾水入滲[14-16]；入滲率和濕潤鋒運移距離隨渾水含沙率的增大而減小[16]；泥沙中小于0.01 mm 顆粒含量越大，阻滲作用越顯著[17]；渾水入滲能力和濕潤鋒運移距離均隨土壤容重的增大而減小[18]。也有研究分析了膜孔灌溉下的渾水入滲規(guī)律[19-21]：LIU 等[22-23]研究了不同渾水肥液濃度和初始含水率下渾水膜孔灌的累積入滲量變化規(guī)律、水分分布和再分布及氮素分布規(guī)律；姜瑞瑞等[24]研究了不同膜孔直徑下多向交匯渾水入滲濕潤體特征及灌水均勻度。另外的研究關注了渾水入滲模型，如引入遺傳算法擬合帶有經驗值的Kostiakov 模型[17]、將渾水入滲形成致密層的阻滲作用歸結為濕潤鋒平均吸力值的Green-Ampt 模型[25]、利用渾水波涌灌大田試驗數(shù)據并結合Kostiakov 模型推導了3 種經驗模型[26]。除農業(yè)領域外，渾水入滲在環(huán)境、工程方面也得到廣泛關注[27-31]。

渾水入滲過程中，所形成致密層阻礙了水分入滲，減小了入滲量。隨著渾水入滲量的減小，渾水中形成致密層的泥沙量也將變少，使得致密層厚度增加速率變緩，土壤入滲能力不至于急速減小。因此，渾水入滲過程是入滲量變化與致密層阻滲能力相互影響的過程[18-19]。除入滲量外，致密層導水能力與其內部孔隙大小、結構有關[32-33]，即形成致密層的渾水泥沙數(shù)量及粒徑（含沙率及泥沙顆粒級配）對致密層導水能力造成影響。通常采用土柱來確定土壤的導水能力[34]，為了便于研究入滲過程，減少不可控因素，已有研究采用飽和土柱進行渾水滲流試驗，來探索該狀況下導水能力的變化規(guī)律[8]。利用數(shù)值模擬方法來探索土壤入滲特性已成為研究熱點，如介飛龍等[35]采用橢圓方程來模擬膜孔灌濕潤鋒形狀，研究了初始含水率對膜孔灌濕潤體特征的影響；馮正江等[36]在使用小波分析和通徑分析法分析Kostiakov 模型標定因子和土壤特性參數(shù)的基礎上，利用多元線性回歸、BP 神經網絡和支持向量機建立估算標定因子的土壤傳遞函數(shù)；王曉彤等[37]利用Hydrus-1D 軟件來模擬分析黃河泥沙充填復墾土壤的入滲和蒸發(fā)特性，優(yōu)化了復墾土壤的夾層結構，為研究設計黃河泥沙夾層式土壤剖面提供了一種可靠方法。

目前，現(xiàn)有研究多集中在渾水入滲特性、影響因素及入滲模型等方面，缺乏渾水入滲致密層形成過程對導水能力影響的研究。在致密層形成過程中，水分入滲上邊界條件不斷變化，土壤的導水能力也隨之改變。研究致密層在土體表面形成后，不同因素下導水能力變化規(guī)律對揭示渾水入滲機理具有十分重要的意義。因此，為了探究渾水入滲對土壤導水能力的影響，本文進行飽和土柱渾水滲流試驗，深入分析多因素（含沙率、泥沙顆粒級配、入滲時間）對致密層形成條件下土體導水能力的影響，建立多因素導水率動態(tài)模型，為進一步揭示渾水入滲規(guī)律提供理論依據。

1 材料與方法

1.1 供試材料

供試土樣采自引黃灌區(qū)0～30 cm 深度范圍的農田土壤，采樣地點分別為西安市灞橋區(qū)、西安北郊及寧夏吳忠市。土壤樣品自然粉干碾碎后，過2 mm 篩以備使用。土壤粒徑采用Mastersizers-2000 型激光粒度分析儀(英國馬爾文儀器公司，測量范圍為0.02～2 000 μm)測定，結果如表1 所示。按照國際制土壤質地分類標準，采自西安灞橋、西安北郊、寧夏吳忠的土樣分別為粉壤土、砂壤土及砂土。

表1 試驗土壤和泥沙粒徑組成Table 1 Particle size composition of tested soils and sediments

入滲試驗渾水中所含泥沙取自涇惠渠灌區(qū)干渠，取回的泥沙經風干過1 mm 篩，人工配置出不同顆粒組成的5 種泥沙，其中3 種不同級配用于多因素試驗，按照黏粒含量大小分別記為J1（低黏粒含量）、J2（中黏粒含量）及J3（高黏粒含量），剩余2 種（Y1 和Y2）用于驗證試驗。

1.2 試驗設計

開展多因素試驗研究渾水含沙率（3%、6%、9%）和渾水泥沙種類（J1、J2、J3）對砂土飽和土柱入滲的影響并建立入滲模型。各組試驗各進行3 次重復。為了驗證多因素渾水入滲模型的可行性，另外設置了8 組試驗。共計17 組試驗，具體方案如表2 所示。

表2 渾水入滲試驗方案Table 2 Experimental scheme of muddy water infiltration

1.3 試驗過程

渾水配置：按照試驗設計選取泥沙種類，并依據含沙率稱量相應質量的泥沙和去離子水，將其倒入馬氏瓶后充分攪拌以配置渾水。

試驗裝置準備：2022 年5 月在西安理工大學農水試驗大廳進行試驗。渾水飽和土柱入滲裝置由土柱和渾水馬氏瓶兩部分組成，如圖1 所示。土柱材質為有機玻璃，內徑為5 cm、高為13 cm，為了通氣和更好收集土柱滲出的水流，底部裝有開孔的墊片，并在土柱底蓋加有管嘴。試驗前，按預定容重分層填裝土壤，層間刮毛，裝土高度為8 cm。為防止土壤從底部小孔中損失，裝土前在底部墊入濾紙。將裝好的土柱在水中浸泡12 h，使其飽和后進行飽和土柱渾水入滲試驗。

圖1 試驗裝置圖Fig.1 Schematic diagram of experimental apparatus

試驗過程：為保持穩(wěn)定的渾水含沙率，將馬氏瓶與磁力攪拌器結合組成渾水馬氏瓶進行入滲。將帶有磁性的攪拌子放入裝有渾水的馬氏瓶中，再將馬氏瓶放置于磁力攪拌器上，利用磁場力使渾水馬氏瓶內的攪拌子進行旋轉來不斷攪拌渾水，保持渾水含沙率的穩(wěn)定。渾水入滲前，為排除初始滲流速度的影響，先進行清水入滲，待入滲穩(wěn)定后，再進行渾水滲流。在管嘴下放置量筒，按照先密后疏的時間間隔量取滲出水的體積并記錄，以計算出不同時間的導水率。根據Darcy 定律[38]，流量與水力梯度成正比，即：

式中q為滲透流量，cm3/min；Kh為導水率，cm/min；i為水力梯度；A為土樣的橫截面積，cm2。

當渾水入滲時，飽和土柱上界面會發(fā)生變化：部分泥沙顆粒隨水分運移進入土壤中，滯留在土層內，而渾水中大部分泥沙顆粒逐漸沉積在土壤表面改變了土壤導水性能，使得導水率不斷發(fā)生變化而無法準確獲得，這里使用一段時間間隔的平均出流量來代替滲透流量，再利用式（1）計算得出不同入滲時刻的土柱導水率值。平均出流量與導水率存在以下關系為[8]

式中Q(t+Δt)為一段時間間隔后出流量，cm3；Q(t)為時間間隔前出流量；Δt為時間間隔，min；二者比值為一段時間間隔平均出流量，近似為滲透流量。

1.4 評價指標

為評價模型準確性，利用統(tǒng)計學中決定系數(shù)（R2），均方根誤差（SRMSE）和相對誤差絕對值均值（SMARE）對模型中計算值與試驗實測值之間的符合度進行評價分析。通常R2越接近于1，SRMSE和SMARE越接近0，表明模型計算精度越高，即實測值與計算值越接近。

2 結果與分析

2.1 多因素影響下導水率動態(tài)變化

根據不同處理渾水飽和土柱入滲試驗，得到每組試驗的導水率值隨入滲時間的變化過程（圖2）。由圖2可以看出，導水率隨入滲的進行不斷減小，在同一入滲時刻，各組試驗累積入滲量存在差異，表明各因素對導水率影響程度不同。

圖2 不同試驗處理導水率隨時間變化Fig.2 Hydraulic conductivity versus time for different treatments

為進一步分析各因素對飽和土柱入滲下導水率的影響，采用多因素方差分析方法分析渾水含沙率與黏粒含量（體積分數(shù)）對導水率的影響，結果見表3。由表3可知，渾水含沙率、各顆粒含量和入滲時間對導水率影響極顯著(P＜0.01)。

表3 不同因素不同水平對導水率的影響Table 3 Effects of different factors and different levels of hydraulic conductivity

采用多元回歸法[39-40]構造渾水含沙率和＜0.002 mm顆粒（黏粒）含量影響下導水率的動態(tài)經驗計算式：

式中S為渾水含沙率，%；N為黏粒含量%；a為導水能力系數(shù)；b、c、d分別為各項因素的指數(shù)。

多元回歸分析結果如式（4）所示，R2為0.853，SRMSE為0.004 cm/min，模型擬合效果良好。

所求經驗模型（式（4））中包含3 個影響因素，不同因素之間單位和數(shù)量級存在差異，為分析各因素對導水率影響的重要程度，對數(shù)據進行標準化處理，處理后b、c、d的標準系數(shù)分別為-0.422、-0.295 和-0.789，表明渾水含沙率、黏粒含量和入滲時間對單位膜孔面積累積入滲量均有影響，其中受入滲時間影響最大，渾水含沙率次之，黏粒含量影響最小；標準系數(shù)均小于0，表明導水率隨入滲時間、渾水含沙率和黏粒含量間增大而減小。

為檢驗式（4）的可靠性，利用試驗處理10、11 對其進行驗證，將試驗實測值和模型計算值進行對比分析，結果見圖3。圖3 中模型（式（4））和實測值間SMARE分別為6.99%和5.93%（小于7%），SRMSE分別為0.008和0.006 cm/min（處理10 和11），總體誤差較?。ㄐ∮?.01 cm/min），說明所建的動態(tài)模型能有效地描述導水率與各因素及入滲時間的量化關系。

圖3 砂土導水率實測值和計算值Fig.3 Measured and calculated values of sand soil hydraulic conductivity

2.2 多因素影響下導水率變化率動態(tài)變化

將導水率對入滲時間求導數(shù)，即得渾水飽和土柱入滲條件下土壤導水率隨入滲時間變化的函數(shù)關系：

式中k為導水率變化率，cm/min2。從式（5）可以看出，導水率變化率隨入滲時間的延長逐漸減小，當時間足夠長時，導水率變化率逐漸趨于0，導水率趨于穩(wěn)定。

由于不同因素對導水率變化率的影響不同，分別取其對含沙率、黏粒含量求偏導數(shù)的絕對值，分析導水率變化率受各因素影響的敏感程度[39]。

通過式（6）和式（7）可分別定量計算出渾水含沙率、黏粒含量對導水率變化率的敏感性指標，敏感性指標越大，相應因素變化對導水率變化率的影響越大。各試驗處理下導水率變化率均隨著各因素增大而減小，故以處理1 為例，計算出各因素對導水率變化率的敏感性指標，點繪其隨各因素變化的曲線，如圖4 所示。

圖4 處理1 導水率變化率敏感性指標與各因素關系Fig.4 Relationship between sensitivity of hydraulic conductivity change rate and each factor in experimental treatment 1

由圖4 可知，各敏感性指標隨著相應因素增大而明顯減小，表明含沙率和黏粒含量變化對導水率變化率均有顯著影響；含沙率敏感性指標為1.64×10-5～2.60×10-3，黏粒含量敏感性指標為4.95×10-5～2.10×10-3，表明黏粒含量和含沙率對導水率變化率影響程度十分接近。入滲時間對各敏感性指標的影響較大，在入滲時間為10 min時，隨著各因素的增大各敏感性指標減小幅度較大，而在入滲時間為30 和90 min 時，各敏感性指標變化幅度明顯減小，在10 min 含沙率敏感性指標隨著含沙率的增大而減小了2.2×10-3，而入滲時間為30 和90 min 僅分別減小了5×10-4和1×10-4。這也說明隨著入滲的進行，導水率變化率對含沙率和黏粒含量的敏感程度逐漸降低，各因素影響下導水率變化趨勢相一致。

2.3 不同土壤質地下動態(tài)導水率模型的適用性

土壤質地對土壤孔隙大小、形態(tài)及分布產生很大的影響，因此不同質地的土壤其飽和導水率一般不相同。

式（3）中包含了導水能力系數(shù)項（a）、含沙率影響項（Sb）、泥沙顆粒級配影響項（Nc）及時間影響項（td）。渾水飽和土柱入滲條件下，部分泥沙進入土體中，而泥沙沉積在土柱表面形成沉積層時會受到入滲速率的影響。為了便于研究，若忽略此影響并不考慮滯留層的作用，則沉積層內部結構是連續(xù)均勻的。當含沙率和泥沙相同時，形成的沉積層結構相同，其厚度隨時間不斷變化，因此導水能力只受到入滲時間和土柱中土壤土質的影響。

對于不同土質的渾水飽和土柱導水率動態(tài)模型，由砂土得出的含沙率影響項（Sb）和泥沙顆粒級配影響項（Nc）中的指數(shù)仍然有效，只有導水能力系數(shù)項（a）和時間影響項（td）不同。即不同土質渾水飽和土柱入滲下，導水率模型為

分別選取處理12 和15（土質分別為砂壤土和粉壤土，含沙率S均為6%，黏粒體積分數(shù)N分別為4.22%和1.01%）來推求式（8）中的導水能力系數(shù)a、時間指數(shù)d。此時式（8）中S、N值均已確定，式（8）可視為導水率隨入滲時間變化的冪函數(shù)：

式中α為擬合系數(shù)，α與導水能力系數(shù)a、含沙率影響項Sb、泥沙顆粒級配影響項Nc存在以下的關系：

圖5 為處理12 和15 導水率隨入滲時間變化圖，對導水率和入滲時間進行冪函數(shù)擬合（見圖5），擬合結果如下：

圖5 處理12 和15 導水率隨入滲時間變化Fig.5 Hydraulic conductivity versus time for treatments 12 and 15

圖5 中擬合結果的決定系數(shù)分別為0.912 和0.930（處理12 和15），SRMSE分別為2×10-3和5×10-5cm/min（處理12 和15），說明式（8）對處理12 和15 擬合結果良好。由此得出砂壤土和粉壤土的時間項指數(shù)d分別為-0.081 和-0.062，分別將處理12 和15 中的S和N值代入式（10）后得到砂壤土和粉壤土的導水能力系數(shù)a分別為0.059 和0.011。

得出由處理12 和15 推求的a、d值后，將其代入式（8），即可得到砂壤土和粉壤土導水率動態(tài)模型，分別如下：

為檢驗式（13）和式（14）的可靠性，分別利用試驗處理13、14、16 及17（處理13、14 土質為砂壤土，處理16、17 土質為粉壤土）對其進行驗證，將試驗實測值和模型計算值進行對比分析，結果見圖6。圖6 中試驗實測值和動態(tài)模型（式（13）和式（14））計算值的SRMSE分別為2.1×10-3、1.2×10-3、6.0×10-4及3.0×10-4cm/min（＜0.01 cm/min），SMASE分 別 為16.58%、14.90%、15.50%及10.46%（處理13、14、16 及17），均小于17%，所建立的砂壤土和粉壤土動態(tài)模型仍能有效地描述導水率與各因素及入滲時間的量化關系。土壤質地為粉壤土（處理16、17）的導水率動態(tài)模型SRMSE、SMASE均小于砂壤土（處理13、14），模型模擬效果更好，這是由于粉壤土較砂壤土為細質土，土體內部結構更為致密，減弱了入滲過程中渾水泥沙顆粒進入土體內部的滯留作用；另外粉壤土飽和導水率小于砂壤土，因此渾水粉壤土飽和土柱入滲速率更小，泥沙沉積過程受入滲速率影響更小，形成的沉積層也更均勻，綜合來看更近于假設條件，因此粉壤土導水率動態(tài)模型模擬效果更好。

圖6 不同處理導水率實測值和計算值Fig.6 Measured and calculated values of hydraulic conductivity for different treatments

3 結論

以土壤質地、含沙率及泥沙種類為影響因素，共進行17 組渾水飽和土柱入滲試驗，利用多因素分析法建立和驗證了渾水砂土飽和土柱入滲下導水率動態(tài)模型，并將此模型推廣至砂壤土和粉壤土，得出如下結論：

1）渾水含沙率、黏粒含量和入滲時間對導水率影響極顯著（P＜0.01），影響程度由大到小依次為：入滲時間、渾水含沙率、黏粒含量，導水率與渾水含沙率、黏粒含量和入滲時間均為負相關；建立了導水率與各影響因素之間的動態(tài)模型，決定系數(shù)（R2）為0.853，均方根誤差（SRMSE）為0.004 cm/min，模型驗證試驗結果中模型計算值與實測值的一致性較好，兩者間的SRMSE小于0.01 cm/min，相對誤差絕對值均值（SMARE）小于7%，說明基于渾水砂土飽和土柱多因素分析法得到的導水率動態(tài)模型可靠性較高。

2）導水率變化率隨入滲時間的延長逐漸減小，渾水含沙率和黏粒含量對導水率變化率影響顯著且影響程度相近，各敏感性指標受入滲時間的影響較大。

3）基于砂土導水率動態(tài)模型和假設條件，建立了適用于砂壤土和粉壤土導水率動態(tài)模型，砂壤土和粉壤土導水率模型R2分別為0.912 和0.930，SRMSE分別為2×10-3和5×10-5cm/min；模型驗證試驗結果中模型計算值與實測值的一致性較好，兩者間的SRMSE小于0.01 cm/min，SMARE小于17%，表明模型能夠較好地反應砂壤土和粉壤土導水率與各影響因素的量化關系。