基于HYDRUS-1D模型的灌排聯(lián)合下的水鹽運移模擬

賈 浩，李寶珠，李文昊，3

（1.石河子大學水利建筑工程學院，新疆石河子832000；2.新疆天業(yè)節(jié)水灌溉股份有限公司，新疆石河子832000；3.現(xiàn)代節(jié)水灌溉兵團重點實驗室，新疆石河子832000）

大面積實施節(jié)水灌溉，誘發(fā)土壤次生鹽漬化、林地退化和草地荒漠化等生態(tài)環(huán)境問題[1]。據(jù)統(tǒng)計，全球鹽堿土面積達到95 億hm2，占耕地面積的18%，且每年以150 萬hm2的速度增長[2]，而中國鹽漬土地占全國可利用土地4.8%，其中，新疆鹽堿土面積為2.2×1012hm2，占全國鹽堿土總面積的22%。近20年來，膜下滴灌技術(shù)[3]在新疆廣泛使用，淺灌和勤灌無法將鹽分排除土體導致次生鹽漬化問題出現(xiàn)，農(nóng)田“有灌無排”的現(xiàn)狀，威脅著綠洲糧食安全和生態(tài)系統(tǒng)的穩(wěn)定[4]。

為了解決此類問題，新疆采取了冬春灌、水利措施等方案。但根據(jù)“鹽隨水來，鹽隨水去”的原理，只有灌排系統(tǒng)協(xié)同發(fā)展，土壤就會發(fā)揮最大效益，明溝和暗管被認為是改良鹽堿地的有效水利措施。國內(nèi)外專家對其進行了研究，Hou等[5]開展灌水量和排水溝深度兩因素三水平試驗，指出3 種排水方式對土壤鹽分的總?cè)コ蕿?.7%～13.2%；Tao 等[6]開展了基于土柱的改進式暗管排水試驗，暗管排水流量相比傳統(tǒng)排水流量增加2～3 倍；王振華等[7]開展了膜下滴灌條件下的暗管排水大田試驗，得出間距15m、埋深0.6 m、管徑90 mm 時，土壤脫鹽率最高；劉洪光等[8]開展了明溝排水和暗管排水對比大田試驗，得出暗管排水的脫鹽效率優(yōu)于明溝。國內(nèi)外學者還采用HYDRUS 模型來深入模擬水鹽動態(tài)分布[9-11]，左強等[12]對所建立的數(shù)值模型和水鹽運動參數(shù)進行了校驗，數(shù)據(jù)結(jié)果能很好反應其水鹽動態(tài)分布。

上述研究均表明暗管排水對土壤脫鹽效率較高，且HYDRUS 數(shù)值模型可以較好地模擬水鹽動態(tài)分布，但灌排聯(lián)合下的水鹽運移問題，至今尚不清楚。因此，本文為探索灌排聯(lián)合下的水鹽運移規(guī)律，利用HYDRUS-1D模型對灌排聯(lián)合下的水鹽運移規(guī)律進行模擬驗證，探索灌排聯(lián)合下水鹽運移的機理，揭示灌排聯(lián)合下排水排鹽效果，為干旱區(qū)灌排聯(lián)合模式下的水鹽運移理論研究提供幫助。

1 材料與方法

1.1 試驗材料

本試驗于2018年7-10月在現(xiàn)代節(jié)水灌溉兵團重點實驗室內(nèi)進行。試驗用土取自八師炮臺鎮(zhèn)鹽堿地0～30 cm土壤，初始含鹽量為為12.46 g/kg，通過風干和研磨，之后過篩（5 mm）使其鹽分均勻，分層裝填，每隔10 cm 進行夯實。根據(jù)國際制土壤質(zhì)地分級標準，土壤質(zhì)地為砂壤土，土壤顆分及基本物理性質(zhì)如表1 所示。試驗用水為地下水，其礦化度為0.25 g/L。滴頭采用貼片式滴頭，滴頭流量為1.8 L/h。試驗模型裝置見圖1，利用馬氏瓶（高65 cm，容積5 102.5 cm3）作為水源，5 mm 厚有機玻璃管制作而成的土柱（半徑20 cm，裝置高130 cm），土柱表層覆蓋新疆天業(yè)公司生產(chǎn)的聚乙烯塑料地膜，水表用來記錄灌水量。底部開3～5 mm 的孔，開孔率為5%，底部鋪設高度為5 cm 的砂石反濾層。排水管材質(zhì)為PVC 管、管徑分為d=5、7 cm，排水管上開直徑4或2 mm的圓孔，開孔率為6%，排水管外包一層透水無紡布，無紡布外覆3～5 cm 細砂和粗砂作為濾層，細砂平均粒徑為0.35～0.5 mm，粗砂的平均粒徑為5～10 mm[22]。排水管距土柱表層80 cm。排水管通過圓孔與圓柱體相連，用玻璃膠進行密封防止漏水，排水管坡度設置為5‰，使排出的水能夠順利地沿排水管坡度方向流入收集器內(nèi)，土柱放置于盛放有水塑料盆內(nèi)，用來模擬地下水位。

表1 供試土壤物理參數(shù)Tab.1 Test soil physical parameters

1.2 試驗設計及方法

試驗采用灌水量、排水模式二因素三水平完全處理，灌水設置3個水平：灌水定額分別為15、25、35 L（分別標記為W15、W25、W35），排水模式設置3 種：排水管d=7 cm、排水管d=5 cm、無排水管（分別標記為P7、P5、P0）。為保證水分和鹽分在土壤內(nèi)重新分布且時空變化較小，每次灌水時間間隔基本相同，共進行3次灌水淋洗。排水結(jié)束后水分經(jīng)過重分布（24 h）后立即進行取樣，取樣點分別在土柱中心處、沿暗管水平方向距離滴頭±15 cm 處，共3 個取樣點，取樣深度為80 cm，每隔10 cm取樣一個，重復三次，分析灌排聯(lián)合對土壤水鹽運移及土壤脫鹽效果的影響。將所取土樣在105 ℃下烘干8 h，用烘干法測定其含水量，將烘干土樣研磨、過0.5 mm 土壤篩，按土水比為1：5 配置成溶液，用玻璃棒攪拌、靜置使上層液澄清，用電導率儀測定其電導率，用干燥殘渣法確定土壤含鹽量與電導率之間的標定關(guān)系式，如公式(1)所示。

式中：S為土壤含鹽量，%；EC為電導率，μs/cm。

為了進一步分析土壤脫鹽狀況，將每次灌水后土樣含鹽量與初始含鹽量進行比較，計算含鹽量相對變化率，根據(jù)每次取樣測得的土壤鹽分質(zhì)量分數(shù)，得出土壤脫鹽率公式如下：

式中：η 為脫鹽率%；Wh為不同深度的含鹽量，%；W0為初始含鹽量，%。

1.3 模型建立

HYDRUS 模型是一種有效的數(shù)值模型工具，是利用土壤物理參數(shù)模擬水、熱及溶質(zhì)在飽和非飽和土壤中的一維（二維）運動的有限元計算機模型[13]。

土壤水分運動方程：

式中：θ 為土壤體積含水率，cm3/cm3；K(θ)為非飽和土壤導水率，cm/d；t 為時間，d；z 為土壤深度，cm。

土壤水分特征特征曲線和非飽和導水率用van Genuchten方程擬合：

式中：Ks為土壤飽和導水率，cm/d；θe土壤相對飽和度；θr為土壤剩余體積含水率；θs為土壤飽和體積含水率；θ(h)為土壤體積含水率，cm3/cm3；h 為負壓水頭，cm；K(θ)為土壤非飽和導水率，cm/h；n、m、α 均為經(jīng)驗參數(shù)；其中m=1-1/n，α 是與土壤物理性質(zhì)有關(guān)的參數(shù)；l 為經(jīng)驗擬合參數(shù)，通常取平均值0.5。

鹽分運移基本方程：

式中：C 為溶質(zhì)濃度，g/cm3；qi為水通量，cm/h；Dij為擴散度，cm2/h；xi為空間坐標（i=1，2），x1=x，x2=z，D11=Dxx，D12=Dxz。

1.4 模型初始條件和邊界條件

土柱上邊界是地膜下面有貼片滴頭，采用第一類邊界條件，確定上、下邊界的初始含水率分別為0.028 cm3/cm3、0.39 cm3/cm3，上下邊界初始含鹽量均為12.46 g/ kg。在整個灌水過程中，上邊界為常壓，下邊界的排水管設置為滲透邊界，灌水結(jié)束，上邊界為零通量。將土壤含水率和土壤含鹽量的實測值及模擬值分別進行對比，相應調(diào)整土壤水力特性參數(shù)，使兩者充分接近為止在模型中輸入供試土壤物理參數(shù)，修正的模型參數(shù)如表2所示。根據(jù)土柱試驗計算土壤鹽分的縱向彌散系數(shù)，土壤鹽分彌散系數(shù)表現(xiàn)了土壤中溶質(zhì)隨著水分運動的變化規(guī)律，修正的模型參數(shù)如表3所示。

表2 土壤物理特性參數(shù)Tab.2 Soil physical property parameters

表3 土壤鹽分橫縱彌散系數(shù)Tab.3 Horizontal and vertical dispersion coefficients of soil salinity

1.5 模型模擬

將土柱簡化為計算區(qū)域為高(垂直方向)120 cm 的一維模型，基于本模型來探索干旱區(qū)灌排聯(lián)合條件下土壤水鹽的動態(tài)變化規(guī)律。由于排水管埋深為80 cm，所以本模型模擬0～80 cm 深度范圍土壤水分和鹽分變化特征，模擬時間共計45 d，采用變時間步長剖分方式，根據(jù)收斂迭代次數(shù)調(diào)整時間步長。設定初始時間步長為0.00l d，最小步長為0.000 l d，最大步長為5 d；土壤含水量容許偏差為0.001，壓力水頭容許偏差為1 cm。模型基本參數(shù)不變，以每次灌水量和灌水時間的試驗數(shù)據(jù)進行模型的驗證。模擬值和實測值的吻合程度采用決定系數(shù)R2和均方根誤差RMSE指標進行評價，計算公式如下：

式中：Si、Mi分別為模擬值和實測值；N 為樣本數(shù)，無量綱。

2 結(jié)果與分析

2.1 模型驗證

通過Hydrus-1D軟件構(gòu)建其模型，輸入各項參數(shù)，設置時間、空間離散化處理參數(shù)等進行運算，并進行模型可靠性驗證分析。由圖2 看出，第二次灌水后灌排聯(lián)合下W35P7處理滴頭處的土壤水分、鹽分模擬值與實測值的R2和RMSE 分別為0.975 和0.286、0.997 和0.217。同時第二次灌水后距離滴頭±15 cm 處土壤水分、鹽分模擬值與實測值跟滴頭處趨勢相同，規(guī)律一致。

通過含水率變化圖可以看出，其土柱土壤水分運動規(guī)律明顯，豎直方向增加趨勢明顯，變化率是由高到低再逐漸升高，從0～20 cm 土層內(nèi)的R2和RMSE 看出，模擬值小于實測值，是由于土柱表層覆膜和取樣數(shù)量較少導致，而20～80 cm土層內(nèi)的R2和RMSE看出，模擬值和實測值趨勢一致且相當吻合，這是土壤在裝填時，經(jīng)過篩分和分層填裝的結(jié)果，而這也側(cè)面印證了符合達西定律、土壤水鹽運動方程的均勻連續(xù)性假設。通過含鹽量變化圖可以看出，0～40 cm 土層內(nèi)含鹽量下降明顯，實測值和模擬值吻合較好，40～80 cm 土層內(nèi)含鹽量變化趨勢減緩，且實測值都小于模擬值，是重力作用的水流運動導致的，同時也是土壤鹽分運動遵循對流彌散理論的結(jié)果。第二次灌水后各處理的土壤水分、鹽分模擬值與實測值的R2和RMSE 以表格的形式展示，如表4 所示，第一次和第三次灌水后的各處理的土壤水分、鹽分模擬值與實測值的R2和RMSE 也具有相似的規(guī)律。由此看出，隨著灌水量的增大，不同處理滴頭處土壤鹽分和水分實測值與模擬值的RMSE均變小，而R2均變大，同樣隨著排水效率的增大，兩者之間也表現(xiàn)出較好的一致性。

表4 灌排聯(lián)合下滴頭處水鹽運移模擬實測值和模擬值驗證Tab.4 Simulated measured value and simulated value verification of water and salt transport at the drip joint under irrigation and drainage

2.2 灌排聯(lián)合模式下土壤水分動態(tài)分布規(guī)律

圖3表示第三次灌水后不同排水措施下滴頭處土壤含水率分布特征變化。由圖3（a）可知，當無排水措施時，不同灌水量下的土壤含水率變化，豎直方向上總體表現(xiàn)為：先高變低再變高后穩(wěn)定的變化趨勢，隨灌水量的增加各土層內(nèi)含水量均增大，0～80 cm土層內(nèi)W25、W35含水率趨勢相同，但是0～40 cm 土層內(nèi)W35含水率比W25平均增大4%，40～80 cm 土層內(nèi)W35含水率比W25平均增大1.5%，而W15含水率曲線因為低灌水量，趨勢與其他兩種變化規(guī)律不同，20～50 cm 土層內(nèi)含水率最大。由圖3（b）可知，當排水管d=5 cm 時，W15、W25含水率曲線趨勢一致，40～80 cm 土層內(nèi)W25含水率比W15平均增大2.5%，W35含水率比W15平均增大3.05%，而W35含水率曲線因為高灌水量，大于20 cm 土層內(nèi)含水率變化趨勢穩(wěn)定。由圖3（c）可知，當排水管d=7 cm時，0～40 cm土層內(nèi)W15、W35含水率曲線變化趨勢相同，W35含水率比W15平均增大4.26%，40～80 cm 土層內(nèi)W25、W35含水率曲線趨勢一致，W35含水率比W25平均增大1.36%。由于排水管的作用，使得40～80 cm 土層內(nèi)W35P7比W35P5平均降低0.89%，W35P7比W35P0平均降低1.56%，其中越接近排水管位置的含水率變化越明顯。3 種灌水量下土壤含水量隨土壤深度的增加呈現(xiàn)不一樣的變化趨勢，排水管d=7 cm 的處理排水效果優(yōu)于排水管d=5 cm，且含水率的變化率都顯著高于無排水管的處理（P＜0.05）。綜上所述，通過第一次灌水和第二次灌水取樣的數(shù)據(jù)也可以得到相同的變化趨勢，同時，排水管上方與排水管兩側(cè)旁各土層內(nèi)含水量具有一致的變化規(guī)律，三種灌水量下排水管周圍各土層內(nèi)含水量呈穩(wěn)定的變化趨勢，灌溉水在重力作用下不斷向下遷移逐漸滲入暗管，這是由于土壤和濾層結(jié)構(gòu)形成兩個不同的界面，砂石濾層導水率較大，濾層與土壤界面處匯聚的水分滲入排水管，過多水分在重力作用下逐漸向下遷移。

2.3 灌排聯(lián)合模式下土壤鹽分動態(tài)分布規(guī)律

圖4表示第三次灌后不同處理下土壤鹽分分布特征及運移規(guī)律，其中W15P0表示在無排水措施下的土壤鹽分運移規(guī)律體現(xiàn)，W35P7表示在最大排水量的土壤鹽分分布規(guī)律體現(xiàn)，不同的灌水量對土柱土壤的濕潤范圍和含水率影響有顯著的差異，而根據(jù)“鹽隨水來，鹽隨水去”的原理，土壤鹽分的運動受飽和水運動影響，且土壤鹽分運動遵循對流彌散理論。由圖圖4（a）可知，在無排水措施下，0～80 cm 土層內(nèi)鹽分含量變化速率緩慢，其中0～40 cm 土層內(nèi)含鹽量降低速率要高于40～80 cm土層內(nèi)含鹽量。灌水量越小40～80 cm土層內(nèi)鹽分含量越大，說明當灌水量較小時，只有少部分鹽分被水流帶入下層土壤，且溶解在地下水中的鹽分較小。在0～20 cm土層內(nèi)含鹽量呈“V”型分布。由圖4（b）可知，在排水管d=7 cm 時，0～80 cm土層內(nèi)鹽分含量變化速率較快，其中70～80 cm土層內(nèi)含鹽量平均要比60～70 cm 土層內(nèi)含鹽量高1.69 g/kg，與無排水措施處理相比，0～40 cm、40～80 cm 土層內(nèi)平均含鹽量分別降低0.527、0.824 g/kg。綜上所述，灌水量相同時，排水措施成為影響水鹽運移的關(guān)鍵，導致土壤含鹽量變化效果不同，總體表現(xiàn)為：排水管d=7 cm＞排水管d=5 cm＞無排水管。排水管以上土層內(nèi)含鹽量變化趨勢一致，垂直方向上沿濕潤峰方向逐漸增加，水平方向上差異不顯著。因此灌水量不同對土壤鹽分的淋洗作用不同。

2.4 灌排聯(lián)合模式下土壤脫鹽效果影響

表5 表示灌排聯(lián)合下不同處理各土層內(nèi)土壤脫鹽率變化，通過方差分析探討不同處理下各土層內(nèi)土壤平均脫鹽率的差異性。當排水管d=7 cm 時，W15P7處理與W25P7處理、W35P7處理相比0～20 cm 和0～80 cm 土層內(nèi)脫鹽率存在顯著性差異，其他各土層內(nèi)脫鹽率在三個處理下均存在顯著性差異，各土層內(nèi)脫鹽率均在W35P7處理下達到最大值。當無排水措施時，即d=0時，灌水量對各土層內(nèi)脫鹽率影響較大，增大灌水量可以提高土壤脫鹽率。說明增加灌水量對鹽堿土鹽分具有較好的淋洗作用，整體來看各土層內(nèi)脫鹽率隨著灌水量的增加而增加，隨著排水管直徑的增大而增加，同樣隨著土壤深度越深土壤脫鹽率變小，表層0～20 cm土層內(nèi)脫鹽率最大。

表5 灌排聯(lián)合下不同處理各土層內(nèi)土壤脫鹽率Tab.5 Desalination rate of soil in different soil layers under different irrigation and drainage combinations

3 結(jié) 論

（1）通過灌排聯(lián)合下的水鹽運移模擬試驗，3種灌水量下土壤含水率隨土壤深度的增加呈現(xiàn)異樣趨勢，排水管d=7 cm的處理排水效果優(yōu)于排水管d=5 cm，且含水率的變化率都顯著高于無排水管的處理（P＜0.05）；灌水量相同時，土壤含鹽量變化效果總體表現(xiàn)為：排水管d=7 cm＞排水管d=5 cm＞無排水管；經(jīng)過3 次灌水淋洗后，各處理0～20 cm 土層內(nèi)土壤鹽分淋洗效果最顯著，灌水定額35 L 及排水管d=7cm 處理下脫鹽率達到最大值83.69%，灌排聯(lián)合模式可有效地降低土壤鹽分含量，可高效地改良鹽堿地。

（2）利用HYDRUS-1D模型對灌排聯(lián)合模式下的水鹽運移規(guī)律進行了模擬，利用實測數(shù)據(jù)和模擬數(shù)據(jù)對比分析，通過R2和RMSE 得出，隨著灌水量、灌水次數(shù)、土壤深度和排水措施變化，土壤水鹽含量的模擬值和實測值有較好的擬合性，模型具有一定的可靠性，這將為干旱區(qū)灌排聯(lián)合模式下的水鹽運移理論研究提供幫助。