胡文同，栗現(xiàn)文，江思珉

（1.武漢大學水資源與水電工程科學國家重點實驗室，武漢430072；2.西北農林科技大學水利與建筑工程學院，陜西楊凌712100；3.同濟大學土木工程學院，上海200092）

0 引言

我國是一個農業(yè)大國，又是一個水資源短缺且時空分布極不平衡的國家，淡水資源短缺已成為影響農業(yè)可持續(xù)發(fā)展的主要問題。同時，我國微咸水資源儲量多、分布廣、開采條件較好[1,2]，高效利用微咸水資源已成為緩解淡水供需矛盾的重要途徑。膜下滴灌技術被認為是目前最合理有效的微咸水開發(fā)利用技術[3,4]，但使用不當會造成土壤積鹽，不利于農業(yè)的可持續(xù)發(fā)展[5-7]。有學者認為微咸水膜下滴灌農田是否積鹽的臨界灌溉水礦化度為3.0 g/L[8]。

目前，對層狀土壤水鹽運移規(guī)律的研究，主要集中在黏土夾層和夾沙層[9,10]，而對犁底層的研究相對不足。犁底層表現(xiàn)為土壤密實程度變化[11]，改變土壤水力性質，影響土壤水分及溶質運移和分布。研究發(fā)現(xiàn)犁底層改變層狀土壤水分入滲過程、濕潤鋒、累積入滲量[12]。張建兵等[13]發(fā)現(xiàn)犁底層能夠阻礙鹽分運移而匯聚于該層，而樊慧敏等[14]通過對渭北農田取樣發(fā)現(xiàn)，土壤容重對含鹽量有明顯的異位效應。但在膜下滴灌條件下，不同容重犁底層對土壤水鹽運移分布的影響尚不明確。

隨著計算機技術的快速發(fā)展，數(shù)學模擬模型HYDRUS 已廣泛應用于滴灌及層狀土壤水鹽運移分布模擬研究，為探究土壤水鹽變化提供了技術支持[12,15,16]。如Roberts 等[17]采用HYDRUS-2D/3D 模擬地下滴灌土壤水鹽運動，分析不同滴灌管埋深和灌溉水礦化度的地表積鹽問題；Xu 等[18]通過田間灌溉試驗與HYDRUS-2D數(shù)值模擬發(fā)現(xiàn)灌溉對土壤脫鹽有明顯影響，模擬值與實測值的偏差小于0.5%；楊霞等[12]通過東北寒區(qū)分層土壤、均質土壤室內入滲試驗與HYDRUS-2D數(shù)值模擬得出，犁底層抑制土壤水分下滲，模擬值與實測值吻合；李遠[19]利用HYDRUS 進行三維入滲試驗發(fā)現(xiàn)，土壤容重越大，垂直濕潤鋒越小，水平濕潤鋒越大、脫鹽系數(shù)越小。但很少有學者利用HYDRUS-2D系統(tǒng)研究犁底層容重對微咸水膜下滴灌土壤水鹽運移分布的影響。

本研究利用HYDRUS-2D軟件建立相同土質不同犁底層容重的二維層狀土壤水鹽運移模型，分析犁底層容重對微咸水膜下滴灌土壤水鹽運移分布的影響，并通過模擬結果的溶質示蹤量化不同模擬情景的排水來源，可為田間合理利用微咸水及灌溉排水提供參考。

1 模型與方法

1.1 概念模型

在HYDRUS-2D 的Geometry 模塊構建水平方向為150 cm(田間一膜寬度，灌水模式為“一膜雙管”)，垂直深度為100 cm 的土壤水流、溶質二維運移模擬域，如圖1所示。上邊界不考慮蒸發(fā)和降水，滴頭處設為變流量邊界；下邊界不考慮地下水的影響，設為自由排水邊界；左右邊界為零流量邊界。溶質運移模型邊界條件均為第3類邊界條件，在滴頭處，滴灌時溶質通量為灌溉水溶液濃度，無滴灌時為零；零流量邊界處溶質通量也為零。

圖1 土壤水流、溶質運移概念模型Fig.1 Conceptual model of soil water and solute transport

1.2 水分運動方程

模型采用修改后的Richards方程描述土壤水分運動[20]：

式中：θ為土壤體積含水率；h為土壤壓力水頭，cm；t為入滲的時間，h；x為橫向坐標，cm；z為垂向坐標，cm；K(h)為非飽和導水率，cm/h。

土壤水分運動初始條件為：

式中：θ(x，z，0)、θi(x，z)為0 時刻i節(jié)點土壤體積含水率，假設各模擬情景初始土壤體積含水率為0.15。

土壤水分運動邊界條件：

式中：qd為滴頭流量，cm/h，有灌水時為100 cm/h，無灌水時為0。

土壤水力特性根據(jù)van Genuchten-Mualem 模型[21]表示為：

式中：θs為土壤飽和體積含水率；θr為土壤剩余體積含水率；Ks為土壤飽和水力傳導度，cm/h；Se為土壤有效飽和度，%；α、n、m、l均為擬合參數(shù)，l通常取0.5。

1.3 溶質運移方程

在溶質運移模擬中考慮對流和彌散2種運動的影響。溶質運移方程如下[20]：

式中：c為土壤溶液中溶質濃度，g/L；D為彌散系數(shù)，cm2/h；q為體積通量，cm/h。

土壤溶質運移初始條件為：

式中：c(x,z,0)、ci(x,z)為0 時刻i節(jié)點土壤水礦化度，g/L，假設各模擬情景土壤水初始礦化度為1.0 g/L。

土壤溶質運移邊界條件：

式中：cd為滴灌用水礦化度，g/L，為3.0 g/L。

1.4 參數(shù)設定

1.4.1 土壤性質

模型土壤質地選自文獻[19]，其中沙粒(0.05～2 mm)、粉粒(0.05～0.002 mm)和黏粒(＜0.002 mm)含量分別為51.4%、26.6%和22%。

1.4.2 土壤水力參數(shù)

本研究利用HYDRUS-2D 模型中的Rosseta 模塊，根據(jù)土壤質地(沙粒、粉粒、黏粒含量)和土壤容重，預測不同容重土層土壤水力參數(shù)。

1.4.3 溶質運移參數(shù)

研究表明，溶質運移彌散度與土壤容重呈線性相關[22]。因各層位土壤為各向同性，故認為縱、橫向彌散度相同。彌散度(α，cm)與土壤容重(ρ，g/cm3)的關系參照文獻[19]：

1.5 情景設定

為系統(tǒng)探討犁底層容重對土壤水鹽運移分布影響，設置4種情景(見表1)，其中S1 為均質情景，S2～S4 情景20～40 cm 處犁底層容重不斷增大。由HYDRUS-2D 的Rosseta 模塊和式(11)，確定不同容重土層的土壤水力和彌散度參數(shù)(見表2)。

表1 不同犁底層容重模擬情景土壤性質Tab.1 Soil properties of different plow pan bulk density simulation scenarios

表2 模型水力及溶質運移參數(shù)Tab.2 Hydraulic and solute transport parameters of the model

模擬時間設定為63 h，其中，灌水時間15 h，重分布時間48 h。在進行迭代計算時，離散時間單位為h，初始時間步長0.000 1，最短時間步長0.000 1，最大時間步長為120，θ迭代精度為0.01，最大迭代次數(shù)取10。設置網(wǎng)格尺寸為2 cm，在滴頭和和土壤容重突變處加密單元格。在左側滴頭及兩滴頭中間10、30、50、70、90 cm 深度處設置觀測點(見圖1)，用于觀測各土層土壤水鹽動態(tài)變化。因膜下滴灌土壤鹽分主要集中在0～60 cm 土層[23，24]，在模型10、30、50 cm 深度每隔10 cm 設置一個觀測點，每個深度設置16 個觀測點，分別用于對耕作層（0～20 cm）、犁底層（20～40 cm）和底土層（40～60 cm）水鹽狀況進行統(tǒng)計分析。

1.6 敏感度指數(shù)

土壤水力參數(shù)對土壤水鹽運移會產生重要影響。前人研究表明，土壤水鹽運移對θs、n、Ks、αL的敏感性較強，因此，選擇上述參數(shù)進行局部敏感性分析。采用擾動分析法，將各參數(shù)在基準參數(shù)的基礎上分別上下擾動5%和10%。定義敏感度指數(shù)L為：

式中：n為擾動次數(shù)；Δi為擾動值，分別取-10%、-5%、+5%、+10%；m為觀測點數(shù)（分別為N1、N2、N3、N6、N7、N8），m=6；p為觀測時間點，分別為灌水結束、灌后24 h 及48 h；θ0為基準參數(shù)情景對應的土壤含水量。

1.7 數(shù)據(jù)統(tǒng)計與情景

利用Excel 對HYDRUS-2D 輸出文件的土壤水鹽數(shù)據(jù)進行匯總和計算；采用SPSS 進行統(tǒng)計特征值分析、單因素方差分析和配對樣本T檢驗。

2 結果與討論

2.1 犁底層容重對土壤水鹽運移規(guī)律的影響

2.1.1 犁底層容重對土壤水分運移規(guī)律的影響

模擬結果顯示飽和體積含水量(θs)隨土壤容重(ρ)增加而降低(θs=-0.252 6ρ+0.773 9，R2=0.999 5)。圖2 為模擬過程各情景滴頭及2 滴頭中間不同深度θ變化過程。對比發(fā)現(xiàn)，處于相同層位的觀測點有相似的變化規(guī)律：飽和含水率相同且模擬結束后θ變化規(guī)律基本一致。各情景N6～N8 處θ開始變化所需時間明顯長于相同層位觀測點N1～N3，但N9、N10處θ開始變化所需時間與相同層位觀測點N4和N5基本一致，說明滴灌條件下，土層水分經過犁底層后逐漸呈層狀運動。灌水結束后，犁底層容重越大，耕作層N1 和N6 處θ越大，犁底層N2和N7 處Se越高，說明犁底層對該層及其上層土壤涵養(yǎng)水分具有積極作用，且容重越大效果越明顯。其原因為犁底層容重越大，土壤穩(wěn)定入滲率越小[25]，等量水分下滲需要的時間越長。底土層N3 處θ開始變化所需時間：S1(6.5 h)＜S2(6.7 h)＜S3(7.1 h)＜S4(7.2 h)，但開始變化與達到飽和的時間間隔為：S1(5.6 h)＞S2(5.0 h)＞S3(4.3 h)＞S4(3.9 h)。這是因為犁底層容重越大，該土層對濕潤鋒的吸力越大，濕潤鋒會暫停在容重突變界面，隨著水分的不斷下滲，吸力減弱，土壤水分才進入下層土壤中；同時，犁底層容重越大，該土層θs越低，土壤達到飽和所需灌水量越少，犁底層土壤入滲率達到穩(wěn)定(最大)所需時間越短，故下層土壤達到飽和所需時間越短。

圖2 各情景滴頭位置N1～N5和2滴頭中間N6～N10處θ變化曲線Fig.2 Dynamic changes of θ at N1～N5 at emitter position and N6～N10 between the two emitters in each scenario

2.1.2 濕潤鋒動態(tài)變化

圖3顯示各情景左側滴頭濕潤鋒隨時間動態(tài)變化。隨灌水時間增加，各情景濕潤鋒水平和垂直推進距離均逐漸增加。滴灌后期，模型邊界土壤水分垂直運移速度明顯增加，這是滴灌和不透水邊界共同作用的結果；隨滴灌進行，2滴頭濕潤鋒發(fā)生交匯，圖3中右側交匯區(qū)土壤水分運移速度較左側明顯增加。濕潤鋒垂直推進距離主要表現(xiàn)出3 個階段：各情景前3 h 濕潤鋒垂直推進距離基本相同；第4～7 h 隨犁底層容重增加而降低；9 h 后隨犁底層容重增加而增加；濕潤鋒水平推進距離表現(xiàn)出2 個階段：前3 h 基本相同，之后濕潤鋒水平推進距離或模型邊界垂直運移深度均隨犁底層容重增加明顯增加。有研究發(fā)現(xiàn)微潤灌溉土壤容重越大，濕潤鋒向下運移速度越慢[26]，沙土中的黏土夾層可以減小濕潤鋒的垂直推進速度[27]，這均與本研究有相似之處。

圖3 濕潤鋒動態(tài)變化Fig.3 Dynamic changes of wet front

2.1.3 犁底層容重對土壤鹽分運移規(guī)律的影響

圖4 為各情景N1～N5、N6～N10 土壤水礦化度隨時間變化規(guī)律。隨土層深度增加，滴頭位置剖面土壤水礦化度開始變化時間不斷變長，且最大值不斷降低。犁底層容重越大，滴頭下方N2～N4 土壤水礦化度越低，2 滴頭中間N6～N10 土壤水礦化度越高。其原因為犁底層阻礙土壤鹽分隨水分垂直運動，促進其水平運動，容重越大效果越明顯，有學者通過溝灌研究發(fā)現(xiàn)，土壤容重越大，鹽分運動受到的阻礙越強[28]。各情景耕作層N1和均質條件犁底層N2處土壤水礦化度均在灌水結束前達到入滲水礦化度，之后一直維持不變，S2～S4模擬時段內犁底層N2 處土壤水礦化度最大值分別是S1(3.0 g/L)的97.67%(2.93 g/L)、95.33%(2.86 g/L)和93.67%(2.81 g/L)；底土層N3 和N4 在灌水結束前后均有增長，S1 底土層N3 處土壤水礦化度最大值(2.86 g/L)分別是S2～S4 的1.05(2.73 g/L)、1.08(2.64 g/L)和1.16 倍(2.59 g/L)。對比2 滴頭中間剖面：各情景耕作層N6 處土壤水礦化度最先開始變化，但觀測后期卻低于犁底層N7 處。其原因為犁底層容重越大，灌溉水到達犁底層后在耕作層及犁底層交界面的橫向運動速度越快，由于“鹽隨水動”的特性，N7處土壤水礦化度增加越快。

圖4 各情景滴頭位置N1～N5和2滴頭中間N6～N10處土壤水礦化度變化曲線Fig.4 Dynamic changes of soil water salinity at N1～N5 at emitter position and at N6～N10 between the two emitters in each scenario

2.1.4 土壤水礦化度空間分布

模擬結束時土壤水礦化度等值線如圖5所示?？梢钥闯?，受微咸水膜下滴灌影響，滴頭下方土壤水礦化度最高，因此在含鹽量較低土壤中用微咸水膜下滴灌時，土壤鹽分會優(yōu)先在滴頭下方聚集。犁底層顯著影響土壤水礦化度等值線，土壤水礦化度高值區(qū)垂直分布深度隨犁底層容重增加明顯減小，水平方向則隨犁底層容重增加而增加，可見犁底層具有和黏土夾層一樣的隔水滯鹽作用[29]。有研究發(fā)現(xiàn)土層之間的滲透性差異越大，黏土夾層的滯鹽作用越強[30]。本研究中犁底層容重越大，鹽分往深層運移越慢。滴頭下方土壤水礦化度變化界面深度均約90 cm，但2 滴頭中間變化界面存在差別，表現(xiàn)為隨犁底層容重增加而增加，但均小于滴頭下方土壤水礦化度。不過也有學者研究發(fā)現(xiàn)犁底層土壤鹽分含量較耕作層顯著增大[13,31]，這是長期往復灌溉和蒸發(fā)多種因素綜合作用的結果。

圖5 模擬結束時土壤水礦化度空間分布Fig.5 Spatial distribution of soil water salinity at the end of the simulation

2.2 敏感性分析

根據(jù)上文定義的敏感度指數(shù)，計算土壤水鹽對各參數(shù)的敏感度指數(shù)L（見表3）。從表3 中可以看出，θ對各參數(shù)的敏感度指數(shù)大于土壤水礦化度。θ對各參數(shù)的敏感度指數(shù)為：θs＞n＞Ks，土壤水礦化度對各參數(shù)的敏感度指數(shù)為：θs＞n＞αL＞Ks。敏感度指數(shù)L大于1，說明參數(shù)對輸出變量的影響大；L小于1，說明參數(shù)對輸出變量的影響小[32]。土壤水鹽對4 個參數(shù)的敏感度指數(shù)均小于1，說明其影響較小。

表3 土壤水鹽敏感度指數(shù)Tab.3 Soil water and salt sensitivity index

2.3 土壤水鹽空間變異分析

2.3.1 土壤含水量空間變異分析

由于灌水結束時0～60 cm 深度θ處于飽和狀態(tài)，故僅對灌后重分布48 h 土壤水分進行統(tǒng)計分析(見表4)。重分布48 h，耕作層θ均值隨犁底層容重增加而增加(p＜0.05)；犁底層θ均值雖然逐漸降低，但其Se分別為79.61%、80.65%、83.05%和87.02%，隨犁底層容重增加而增加(p＜0.05)；S2～S4 底土層θ均值隨犁底層容重增加而降低(p＜0.05)，故犁底層不僅對本土層土壤水分重分布產生重要影響，而且會產生異位影響。各情景θ變異系數(shù)隨土壤深度增加而增大；S4各土層變異系數(shù)最低(S3 和S4 耕作層均為0)，S1 各土層變異系數(shù)最大，故犁底層容重的增加可以降低本層及相鄰土層θ的空間變異程度。變異系數(shù)可以反映樣本的空間變異性，低于10%為弱變異，10%～100%為中度變異，大于100%為強變異[33]，各情景θ變異系數(shù)遠小于10%，均為弱變異。

表4 模擬結束時土壤水分統(tǒng)計特征值 %Tab.4 Statistical characteristic value of soil moisture at the end of simulation

2.3.2 土壤水礦化度空間變異分析

表5 為灌水結束及重分布48 h 后土壤水礦化度統(tǒng)計特征值。隨犁底層容重增大，各土層土壤水礦化度最小值逐漸增大，而犁底層和底土層最大值逐漸減小。各情景土壤水礦化度均值表現(xiàn)為：耕作層＞犁底層＞底土層，耕作層土壤水礦化度均值隨犁底層容重增加而增加，S2～S4 和S1 的差異均達顯著水平(p＜0.05)。有研究發(fā)現(xiàn)，上層土壤鹽分與下層土壤容重正相關[14]，這與本研究中犁底層對耕作層土壤水礦化度的影響結果一致；但層狀土壤中的黏土夾層在抑制表土積鹽的同時會增加本土層含鹽量[9]，而本研究中犁底層的土壤水礦化度卻低于均質條件相同位置，可能原因是土質不同及本研究未考慮蒸發(fā)和地下水的影響。灌水結束時，S1 土壤水礦化度變異系數(shù)隨土層深度增加而增加，但其余情景犁底層位置最低，且各土層土壤水礦化度變異系數(shù)隨犁底層容重增加而降低；S4 耕作層、犁底層和底土層土壤水礦化度變異系數(shù)分別是S1 相同土層的61.92%、35.05%和35.51%，因此，犁底層能夠顯著降低土壤鹽分的空間異質性。重分布48 h，各情景各土層土壤水礦化度最小值、最大值、均值有不同程度的增長或不變，且增長量隨土層深度增加而增加；S2～S4 和S1 耕作層土壤水礦化度均值差異仍達到顯著水平(p＜0.05)；各情景各土層土壤水礦化度變異系數(shù)較灌水結束時均有所降低，但降低幅度隨犁底層容重增大而減小。

表5 灌水結束時(模擬灌水時長15 h)及重分布48 h后(模擬時長63 h)土壤水礦化度統(tǒng)計特征值 g/LTab.5 Statistical characteristic values of soil water salinity after irrigation and 48 hours'redistribution

2.4 犁底層容重對排水排鹽的影響

4 種情景模擬結束時累積排水排鹽量見表6。S4 累積排水排鹽量最大，其次是S3，而S1 最??；灌水結束時，犁底層容重越大，透過犁底層到達底土層的水量越大（見圖3），底土層排出的水量越多，但犁底層上部土層的下滲水量越少。在本研究的特定情境下，受犁底層的影響，多透過犁底層的水量大于犁底層導致的上部土層減滲的水量，且其影響隨容重增加而變大，故犁底層容重越大，模型底部累積排水排鹽量越大(p＜0.05)。S1 和S4 累積排水量差為32.157 cm3/cm，遠小于模擬結束時犁底層容重對本土層產生的水量差異（124.2 cm3/cm，犁底層面積×模擬結束時兩情景犁底層位置θ之差），各情景排水礦化度隨犁底層容重增加而增大，但遠小于滴灌用水礦化度。

表6 各情景模擬結束時排水排鹽量Tab.6 Drainage and salt discharge at the end of simulation

通過將礦化度作為溶質示蹤劑，發(fā)現(xiàn)排水中新水所占比例很小(見表7)，且該比例隨犁底層容重增大而略有增加。有學者通過同位素技術和溶質示蹤劑揭示降雨前儲存在土壤中的老水是降雨徑流的主要來源，而雨水主要起驅遣老水的作用[34]?；诒狙芯康哪M域及灌溉情景，灌溉水會驅使背景礦化度較低的土壤水排出，而礦化度較高的灌溉水會滯留在土壤中；故對于不受鹽堿化影響或鹽堿化程度較輕的農田，使用礦化度較低的微咸水進行灌溉仍會加重土壤鹽堿化程度。

表7 排水來源百分比 %Tab.7 Percentage of drainage water

3 結論

本文通過Hydrus-2D軟件模擬犁底層容重對微咸水膜下滴灌土壤水鹽變化、統(tǒng)計特征值及排水排鹽的影響，得出以下主要結論。

（1）在相同土質和滴灌條件下，犁底層容重越大，耕作層土壤水流橫向滲流擴展能力越強，土壤水分下滲速度越慢，底土層含水量變化所需時間越長。犁底層對濕潤鋒垂直推進距離的影響表現(xiàn)為先抑制后促進，而水平方向一直為促進。犁底層顯著影響土壤水礦化度空間分布，該層容重越大，滴頭下方土壤水礦化度越小，2滴頭中間越大。模擬結束時累積排水排鹽量隨犁底層容重增大而增大(p＜0.05)。通過將礦化度作為溶質示蹤劑發(fā)現(xiàn)，所排水量中灌溉水所占比例很小，且隨犁底層容重和累積排水量的增加而增加。

（2）灌水結束時，各情景土壤水礦化度均值隨土層深度增加而降低；均質條件土壤水礦化度變異系數(shù)隨土層深度增加而增加，但其余情景犁底層最低，且各土層土壤水礦化度變異系數(shù)均隨犁底層容重增加而降低。重分布48 h，各情景土壤水礦化度變異系數(shù)較灌水結束時均有所降低，但降低幅度隨犁底層容重增大而減小。

（3）本研究僅考慮了犁底層容重對微咸水膜下滴灌土壤水鹽運移分布的影響，未考慮蒸發(fā)、地下水、植被與犁底層的耦合影響規(guī)律。應用本文研究成果時，還需根據(jù)田間具體條件開展相關試驗對數(shù)值模型參數(shù)進行進一步率定，并在必要時完善模型結構。