基于改進(jìn)的HYDRUS-2D模型暗管排水過程模擬

劉 義,曾文治，敖 暢,雷國(guó)慶,伍靖?jìng)?黃介生

(武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

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

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

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

1 材料與方法

1.1 試驗(yàn)設(shè)計(jì)

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

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

1.2 模型的基本原理

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

在HYDRUS-2D模型中，假定土壤為具有各向異性結(jié)構(gòu)的多孔介質(zhì)，并且忽略空氣和溫度對(duì)土壤水和鹽運(yùn)移的影響。用二維飽和-非飽和Richards方程來描述HYDRUS-2D中的土壤水分運(yùn)動(dòng):

(1)

表1 土壤的物理性質(zhì)

圖1 室內(nèi)試驗(yàn)裝置示意圖Fig.1 Schematic diagram of laboratory experiment device

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

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

(2)

(3)

模型中使用二維對(duì)流彌散方程來描述土壤鹽分的運(yùn)移:

(4)

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

1.3 邊界設(shè)置和初始條件

為了改進(jìn)HYDRUS-2D在暗管排水中的應(yīng)用，考慮了3種暗管邊界的概化方式(圖2)。第一種方案，設(shè)置暗管邊界為原始的滲透邊界(PSB)作為其它兩種方案的對(duì)照組。第二種方案，在暗管周圍1cm范圍內(nèi)添加虛擬土層(VSL)，其作用是模擬無紡布和波紋管的結(jié)合界面。第三種方案，考慮到暗道的實(shí)際開孔率非常小，將暗管邊界設(shè)置為無通量和自由滲透相交替的邊界，用來模擬暗管的實(shí)際開孔面積(AHA)。因?yàn)樵囼?yàn)中使用的淋洗槽的周邊和底部是不透水的，所以將模擬區(qū)域的左右邊界和下邊界設(shè)置為無通量邊界。在本試驗(yàn)中，由于上邊界的水頭隨時(shí)間變化是逐漸減小的，所以根據(jù)實(shí)測(cè)水頭數(shù)據(jù)，上邊界可簡(jiǎn)化為變水頭邊界，當(dāng)上邊界無積水存在時(shí)則變?yōu)榇髿膺吔纭?/p>

模擬區(qū)域內(nèi)各層土壤的初始含水率和含鹽量在填槽前進(jìn)行測(cè)量，假設(shè)土壤含水率和含鹽量在土壤剖面中呈線性分布，具體見表2。

1.4 模型的率定與驗(yàn)證

本次數(shù)值模擬以“小時(shí)(h)”為時(shí)間單位，模擬周期分別為90 h和60 h。其中第一次試驗(yàn)用于模型的率定，第二次試驗(yàn)用于模型的驗(yàn)證。輸入的土壤參數(shù)由HYDRUS-2D內(nèi)嵌的人工神經(jīng)網(wǎng)絡(luò)計(jì)算工具(Rosetta)根據(jù)測(cè)量的每層土壤的容重和粒徑分布來計(jì)算。在參數(shù)的率定方面，首先對(duì)暗管的排水和排鹽過程進(jìn)行率定，然后對(duì)土壤剖面的最終含鹽量和含水率進(jìn)行率定。在模擬過程中，首先提取暗管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的深度處設(shè)置觀察點(diǎn)，并提取觀測(cè)點(diǎn)土壤剖面的最終含鹽量和含水率。最后，通過PEST的多目標(biāo)校正法將實(shí)測(cè)值與模擬值進(jìn)行對(duì)比。當(dāng)模擬值與實(shí)測(cè)值足夠接近時(shí)，則選取率定的參數(shù)值作為模型最終輸入?yún)?shù)值。上述過程由FORTRAN編程完成。率定后的土壤水力參數(shù)如表3所示。在參數(shù)率定之后，將第二個(gè)試驗(yàn)用于模型的驗(yàn)證。

1.5 評(píng)價(jià)指標(biāo)

均方根誤差(RMSE)和決定系數(shù)(R2)這兩個(gè)指標(biāo)用于評(píng)估模型的模擬精度：

(5)

(6)

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

表3 三種邊界率定的土壤水力參數(shù)

2 結(jié)果與分析

2.1 模型模擬結(jié)果

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

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

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

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

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

2.2 虛擬土層相關(guān)參數(shù)的確定

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

Ks=0.019RW+0.0058R2=0.99

(7)

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

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

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

3 討 論

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

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

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

4 結(jié) 論

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

2)在3種模型的各項(xiàng)模擬結(jié)果中，VSL在對(duì)暗管排水和排鹽過程的模擬最為準(zhǔn)確。VSL能夠很好地反映暗管的排水規(guī)律和水鹽運(yùn)移情況。因此，VSL可以對(duì)干旱地區(qū)應(yīng)用暗管排鹽技術(shù)的理論研究提供幫助。

3)基于VSL的模擬試驗(yàn)，發(fā)現(xiàn)虛擬土層的環(huán)寬(RW)與飽和導(dǎo)水率(Ks)之間也存在定量的線性關(guān)系。在不同的模擬情景中，應(yīng)用此公式可以確定虛擬土層中的相關(guān)模型參數(shù)。

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