基于改進的HYDRUS-2D模型暗管排水過程模擬

劉 義,曾文治，敖 暢,雷國慶,伍靖偉,黃介生

(武漢大學水資源與水電工程科學國家重點實驗室，湖北 武漢 430072)

土壤鹽漬化是一個危害全球農業環境的難題，它不僅能夠抑制作物的生長，降低作物的產量和品質，同時還能降低肥料的利用率，從而增加土壤養分流失的風險[1-2]。暗管排水通常被認為是治理土壤鹽漬化的重要水利工程措施[3-4]。暗管的埋深、間距和外包材料等是影響暗管排水排鹽效果的重要因素，國內外學者對此展開了一系列的研究[5-7]。張展羽等[8]在濱海地區模擬了暗管在不同埋深和不同間距下的排水排鹽效果，確定了適宜在濱海地區埋設暗管的技術參數。Bonaiti等[9]對意大利東北部的地下灌溉和暗管排水的效果進行了研究，表明當管道之間的距離為8 m時，暗管排水效果最佳。王洪義等[10]在地下水埋深較淺的大慶鹽堿地開展了不同埋深下的暗管排水降漬試驗，結果表明暗管埋深0.8 m是當地暗管最佳工程布設參數。王振華等[11]研究了滴灌淋洗條件下的暗管排鹽效果，發現當地下水埋深較大時，暗管間距為15 m排鹽效果最佳。Mastrocicco等[12]的研究中也得到了類似的研究結果。Wang等[13]得出淺層暗管排水能夠有效地降低耕作層土壤的含鹽量，從而提高作物的產量。錢穎志等[14]研究了新疆地區田間的暗管排水試驗，得出暗管上方土壤的脫鹽率隨著埋深的增加而增加。Tao等[15]通過在暗管周圍添加具有較高滲透性能的過濾材料改進了暗管的排水能力。Filipovic等[16]研究證明了在暗管周圍包裹礫石能夠顯著提高暗管的排水效果，主要表現在通過暗管排水能夠顯著減少地表積水和地表徑流。然而，以上的這些試驗通常需要耗費大量的人力物力并且周期一般較長，其結果通常是針對當地的土壤環境，所以缺乏普遍適用性。此外，由于涉及暗管布局的參數較多，很難對這些影響因素進行系統定量試驗。因此，有必要通過數值模型快速、系統地研究暗管排水下的水鹽運移過程，從而優化暗管的布局。

目前，考慮暗管布局的常用模型有HYDRUS-2D、DRAINMOD、SWAT等。相比之下，HYDRUS-2D具有靈活的邊界條件可以較為準確地模擬二維非飽和土中的水鹽運移過程。因此，許多研究將HYDRUS-2D應用于暗管排水的模擬。石培君等[17]運用HYDRUS-2D模擬了膜下滴灌條件下的暗管排水排鹽效果，發現模擬結果與試驗結果吻合較好。Castanheira等[18]通過HYDRUS-2D模擬暗管排水試驗發現，利用暗管排水可以有效降低土壤鹽堿化的風險。Adobea等[19]使用HYDRUS-2D通過對田間土壤含水率變化的模擬，確定了一年中可用于田間作業的工作天數。李顯溦等[20-21]通過田間試驗，基于土壤剖面水鹽的率定結果測試了HYDRUS-2D在模擬暗管排水中的適用性。盡管上述研究中模型的模擬精度是可以接受的，但是他們往往僅從土壤剖面水鹽含量的角度來率定模型，而忽略了排水和排鹽過程，這可能會導致HYDRUS-2D在模擬暗管排水方面具有一定的誤差。此外，以往的研究中暗管邊界主要采用的是滲透邊界，然而滲透邊界的排水能力遠大于實際暗管的排水能力，因此在HYDRUS-2D中僅利用滲透邊界很難模擬暗管的排水、鹽過程。所以有必要為暗管開發新的邊界條件以提高HYDRUS-2D的模擬精度。

本研究旨在改進HYDRUS-2D模型在模擬暗管排水過程中的性能。為此，本文提出了虛擬土層(VSL)和實際開孔面積(AHA)兩種方案來代替HYDRUS-2D暗管的原有滲透邊界(PSB),并且利用試驗數據對以上暗管邊界在HYDRUS-2D的模擬精度進行了評價。

1 材料與方法

1.1 試驗設計

試驗在武漢大學水資源與水電工程科學國家重點試驗室進行。試驗設備是由鋼板焊接而成的淋洗槽，其長度、寬度和高度分別為150、51 cm和95 cm(圖 1)。為了控制水深并確保土壤在大水漫灌中得到持續的淋洗，試驗前在淋洗槽的前面和兩個側面距頂部3.5 cm 處鉆了6個溢流孔(φ=7.5 mm)。暗管埋設在淋洗槽的中心距頂部65 cm的位置，開孔率為3%。坡度為2‰，管道的外徑為54 mm，內徑為50 mm。在管道外包裹一層可滲透的無紡布，用來防止土壤顆粒堵塞暗管的排水孔。

試驗所用土壤樣品采集于中國新疆阿拉爾(40°33′58″N，79°49′32″E)，該地區近年來主要采用暗管排水作為土壤鹽漬化的防治技術。土壤粒徑分布通過激光粒度分析儀(Mastersizer 2000，Malvern儀器有限公司，英國)測得；土壤容重通過環刀法采集原狀土壤測得(表1)。試驗前將土壤自然風干，并用2 mm土壤篩除去雜質進行預處理。在室內試驗中為了模擬田間土壤的層狀分布，將填充至淋洗槽中的土壤分成四層。在填土過程中每隔5 cm將土壤壓實1次，保證每次的土壤壓實度相同。試驗共分兩次進行，在兩次試驗中分別向該淋洗槽中注入200 mm和95 mm的水。通過自記式流量計測得暗管的排水量，并且每30 min收集排水樣品以測得電導率。在每次試驗結束后用土鉆采集10、15、20、25、30、35、40、45、50、55、60、65、70 cm和75 cm深度處的土壤樣品，用來測量土壤剖面中水和鹽的分布。土壤體積含水率用烘箱干燥法測得(105℃， 8 h)，土壤電導率用電導率儀(DDS-307，上海雷磁儀器廠)測得。

1.2 模型的基本原理

采用HYDRUS-2D模型模擬暗管排水條件下土壤水鹽的動態變化。HYDRUS-2D可以基于有限元流動方程的數值解來模擬二維或三維軸對稱流動。在前人的研究中HYDRUS-2D常被用于模擬土壤水分和溶質運移，并且模型的模擬結果在許多情況下均表現良好[22-23]。

在HYDRUS-2D模型中，假定土壤為具有各向異性結構的多孔介質，并且忽略空氣和溫度對土壤水和鹽運移的影響。用二維飽和-非飽和Richards方程來描述HYDRUS-2D中的土壤水分運動:

(1)

表1 土壤的物理性質

圖1 室內試驗裝置示意圖Fig.1 Schematic diagram of laboratory experiment device

式中，θ為土壤體積含水率(cm3·cm-3)；K(θ)為非飽和土壤導水率(cm·h-1)；t為從開始入滲到排水結束的時間(h)；z為土壤深度(cm)。

模型采用土壤水力參數vanGenuchten模型來描述土壤水力特性：

(2)

(3)

模型中使用二維對流彌散方程來描述土壤鹽分的運移:

(4)

式中，C為溶質濃度(g·cm-3)；t為從開始入滲到排水結束的時間(h)；qi為在xi方向的達西流速(cm·h-1)；Dij飽和-非飽和水動力彌散系數(cm2·h-1)，i,j分別表示x,z方向；xi,xj為空間坐標。

1.3 邊界設置和初始條件

為了改進HYDRUS-2D在暗管排水中的應用，考慮了3種暗管邊界的概化方式(圖2)。第一種方案，設置暗管邊界為原始的滲透邊界(PSB)作為其它兩種方案的對照組。第二種方案，在暗管周圍1cm范圍內添加虛擬土層(VSL)，其作用是模擬無紡布和波紋管的結合界面。第三種方案，考慮到暗道的實際開孔率非常小，將暗管邊界設置為無通量和自由滲透相交替的邊界，用來模擬暗管的實際開孔面積(AHA)。因為試驗中使用的淋洗槽的周邊和底部是不透水的，所以將模擬區域的左右邊界和下邊界設置為無通量邊界。在本試驗中，由于上邊界的水頭隨時間變化是逐漸減小的，所以根據實測水頭數據，上邊界可簡化為變水頭邊界，當上邊界無積水存在時則變為大氣邊界。

模擬區域內各層土壤的初始含水率和含鹽量在填槽前進行測量，假設土壤含水率和含鹽量在土壤剖面中呈線性分布，具體見表2。

1.4 模型的率定與驗證

本次數值模擬以“小時(h)”為時間單位，模擬周期分別為90 h和60 h。其中第一次試驗用于模型的率定，第二次試驗用于模型的驗證。輸入的土壤參數由HYDRUS-2D內嵌的人工神經網絡計算工具(Rosetta)根據測量的每層土壤的容重和粒徑分布來計算。在參數的率定方面，首先對暗管的排水和排鹽過程進行率定，然后對土壤剖面的最終含鹽量和含水率進行率定。在模擬過程中，首先提取暗管5、10、15、20、25、30、35、40、45、50、55、60、65、70 h和75 h的累積排水量，以及10、20、30、40、50、60 h和70 h的累積排鹽量來率定暗管的排水和排鹽過程。然后，在10、15、20、25、30、35、40、45、50、55、60、65、70 cm和75 cm的深度處設置觀察點，并提取觀測點土壤剖面的最終含鹽量和含水率。最后，通過PEST的多目標校正法將實測值與模擬值進行對比。當模擬值與實測值足夠接近時，則選取率定的參數值作為模型最終輸入參數值。上述過程由FORTRAN編程完成。率定后的土壤水力參數如表3所示。在參數率定之后，將第二個試驗用于模型的驗證。

1.5 評價指標

均方根誤差(RMSE)和決定系數(R2)這兩個指標用于評估模型的模擬精度：

(5)

(6)

表2 土壤含水率和含鹽量的初始值

表3 三種邊界率定的土壤水力參數

2 結果與分析

2.1 模型模擬結果

3種暗管邊界的模擬結果如圖3和圖4所示。對于VSL和AHA來說，在模擬初期，暗管出流量的模擬值均小于暗管出流量的實測值，而在模擬后期，暗管出流量的模擬值則略大于實測值。但是，總體來看兩種模型的排水歷程與實測排水過程較為一致。然而，在PSB中則出現了相反的模擬結果。盡管PSB的最終排水量與實測值一致，但是排水過程的時長卻遠小于實測排水過程。對比模擬值與實測值，VSL的均方根誤差RMSE為9.48 L，小于AHA和PSB中相對應的值，并且VSL的決定系數(R2)為0.930，在3個模型中最高，因此，VSL模擬排水過程時精度最高。排鹽與排水過程相類似(圖3b和圖4b)，PSB的排鹽過程均比其它2個模型結束得早。此外，在3個模型中，VSL在排鹽過程的率定和驗證階段的RMSE最小而且R2最高，這表明VSL在對暗管排鹽過程的模擬中表現最佳。

圖3c、3d和圖4c、4d顯示了3種邊界對土壤剖面的最終含鹽量和土壤含水率的模擬結果。在模型的率定與驗證的過程中，模型均略微低估了0～40 cm

圖3 PSB、VSL和AHA在模型率定階段的表現Fig.3 The performance of PSB, VSL and AHA during model calibration stage

圖4 VSL和AHA在模型驗證階段的表現Fig.4 The performance of VSL and AHA during model validation stage

深度處的土壤脫鹽率，但40～80 cm深度處土壤剖面含鹽量的模擬值與測量值吻合較好。3種模型的含鹽量模擬精度的RMSE和R2范圍分別為0.757～1.075 g·kg-1和0.964～0.995，這表明模型在模擬土壤的鹽分分布時具有較高的精度。對于土壤含水量，3種模型的RMSE均很小，說明3種模型對土壤剖面水分分布的模擬精度較高。總的來說，與PSB和AHA相比，VSL在模擬排水過程、排鹽過程和土壤剖面水鹽方面表現最好。

2.2 虛擬土層相關參數的確定

由前述所得，VSL是改善HYDRUS-2D用于模擬暗管排水的最佳方案。其中，VSL的飽和滲透系數(Ks)和虛擬土層的圓環寬度(RW)是影響模型模擬精度的兩個關鍵參數。為此，我們分析了暗管排水過程、排鹽過程、土壤剖面含鹽率和含水率對參數Ks和RW的敏感性。當參數Ks和RW變化50%分別導致暗管排水過程和排鹽過程變化32.04%～35.82%，24.63%～27.91%；而土壤剖面鹽分分布和水分分布變化范圍為0.11%～8.93%。這些結果表明暗管排水過程和排鹽過程對參數Ks和RW有較高的敏感性，土壤剖面鹽分和水分對這兩個參數的敏感性較低。可見，Ks和RW的變化對排水、排鹽過程影響最大。為了定量表征它們之間的關系，將虛擬土層環寬(RW)分別設置為0.5、1、1.5、2 cm和2.5 cm，然后用PEST率定相應的Ks。發現在不同的RW情景下，通過調整Ks，該模型能夠很好地模擬土壤剖面的鹽分和水分以及暗管的排水和排鹽過程。將所得的RW與其對應的Ks進行線性擬合(圖5)，擬合效果較好，關系如下:

Ks=0.019RW+0.0058R2=0.99

(7)

式中，Ks為虛擬土層的飽和滲透系數(cm·h-1);RW為虛擬土層的寬度(cm)。

該擬合公式的R2為0.99，接近于1，可見該擬合公式能較好地反映虛擬土層的飽和滲透系數與其圓環寬度之間的定量關系。

圖5 虛擬土層的飽和導水率(Ks)與環寬(RW)之間的關系Fig.5 Relationship between saturated hydraulic conductivity(Ks) and ring width (RW) of virtual soil layer

3 討 論

在數值模擬中，邊界設置是建模的關鍵，它關系到模型的模擬精度[24]。目前，HYDRUS-2D中還沒有專門用于反映暗管排水性能的邊界。在以往的研究中，研究者們主要應用滲透邊界(PSB)來模擬暗管排水[17, 20, 25]。然而，當周圍土壤飽和時，在重力水的作用下，滲透邊界會立即排水。滲透邊界的排水面積大于試驗中所用暗管的開孔面積，所以PSB的排水和排鹽過程可以較快完成，這導致模型模擬的排水過程比實測過程要短。因此，在HYDRUS-2D中僅應用滲透邊界只能模擬暗管排水排鹽過程的總趨勢，不能模擬暗管的整個排水過程的動態變化量，這與李顯溦等[20]的研究結果相近。考慮暗管的實際排水面積、物理結構和開孔率較小(3%)等特性，提出了2種方式(VSL和AHA)來概化暗管的邊界，以改進HYDRUS-2D在模擬暗管排水時的精度。

AHA與VSL均能準確模擬土壤剖面中的水分和鹽分分布(圖3c、3d和 圖4c、4d)。然而在模擬排水量與排鹽量方面，AHA和VSL的最終模擬值均高于實際值。這是由于試驗中所用各層土壤的初始含水率和含鹽量可能存在分布不均的現象，而模擬中采用的土壤初始含水率和含鹽量在各層中是均勻分布的，導致了土壤初始水量與鹽量的模擬值大于實際值[26-27]。從而在相同灌水量的條件下，兩種模型的最終排水量與排鹽量大于實測值。在本研究中，為了表征暗管的實際開孔率和實際排水面積，我們設置了AHA方案。與實測值相比，AHA的排水和排鹽過程具有一定的滯后性。其中的原因是，與PSB相比,盡管AHA的排水面積與暗管的實際排水面積相同，但是AHA所應用的部分滲透邊界的導水率仍與PSB相同且大于實際情況，所以較高的導水率仍會導致暗管的排水速率大于實測值[26]。為了降低AHA的排水速率，上層土壤的Ks被調整為0.12 cm·h-1(表3)，這比下層土壤(1.07～2.07 cm·h-1)小得多。由于飽和導水率較小，上層土壤的透水量減少，下層土壤的飽和速率也隨之降低[28]。暗管周圍土層的飽和速率降低，所以AHA中的滲透邊界的排水時間延后，整體上呈現出一定的滯后性。在模型的驗證階段，AHA的暗管出流時間比模型的率定階段要早，這是由于第二次試驗的土壤各層的初始含水量大于第一次試驗，造成了暗管周圍土壤飽和更快。然而就VSL而言，暗管的排水和排鹽的過程與實測過程幾乎一致。與AHA和PSB的模擬精度相比,VSL的R2最大，RMSE最小，可見添加虛擬土層后模型的精度有了很大的改進。這是因為盡管滲透邊界的排水能力較強，但是虛擬土層的飽和導水率Ks本身較低，流經虛擬土層的水流的速度相對較小[29]，所以最終滲透邊界流出的水量較小，接近實際流出水量。可見，可以通過減少虛擬土層的飽和導水率來削弱滲透邊界排水能力，無需對上層土壤的飽和導水率Ks做出較大的改動。從這點來看，VSL更能反映暗管排水過程的實際情況。

本文的研究結果還表明虛擬土層的圓環寬度(RW)與Ks有很好的線性關系。這是由于虛擬土層本身的Ks較低，隨著土層的擴展，滲透邊界附近的水流流入量會越來越小[27]，應該增大暗管周圍虛擬土層的飽和滲透系數(Ks)，使得模擬值與實測值基本相符。虛擬土層區別于實際土壤，虛擬土層是用來表征暗管性能的特殊“土層”，本研究中的虛擬土層所選用的參數能夠較為準確地模擬暗管的真實出流情況。然而在不同的試驗或二維模擬區域中，虛擬土層的參數和環寬往往是不同的。為了簡化虛擬土層參數的選取，公式(7)僅考慮了環寬改變時對飽和導水率的影響并建立了線性關系，利用此公式能夠快速確定不同模擬情景中的虛擬土層參數的取值，不同情景中的模擬結果也驗證該公式的準確性。此外，在建模階段，虛擬土層的環寬(RW)不宜超過埋設暗管的土層厚度，否則可能會對土壤的物理性質和虛擬土層飽和導水率的選取產生較大的影響，并最終導致模型的模擬精度較低。

4 結 論

1)根據暗管的物理結構、實際排水面積和排水速率，提出了兩種策略(VSL和AHA)來概化HYDRUS-2D中的暗管邊界，并通過2次室內暗管排水試驗的實測值對模擬的排水過程、排鹽過程、土壤剖面水分和鹽分分布的精度進行了評估。模擬值與實測值具有很好的一致性。

2)在3種模型的各項模擬結果中，VSL在對暗管排水和排鹽過程的模擬最為準確。VSL能夠很好地反映暗管的排水規律和水鹽運移情況。因此，VSL可以對干旱地區應用暗管排鹽技術的理論研究提供幫助。

3)基于VSL的模擬試驗，發現虛擬土層的環寬(RW)與飽和導水率(Ks)之間也存在定量的線性關系。在不同的模擬情景中，應用此公式可以確定虛擬土層中的相關模型參數。

在本文的研究中，僅考慮了土壤剖面鹽分和水分分布的靜態值來率定模型而忽略了過程中的動態變化，并且邊界研究僅限于暗管邊界。下一步的研究中需要考慮土壤水鹽的動態變化和上邊界的概化形式以進一步提高模型的模擬精度。