基于HYDRUS模型的華北平原小麥種植區(qū)水鹽運(yùn)移模擬*

李 琦,李發(fā)東,**,張秋英,喬云峰,杜 錕,朱 農(nóng),楊 廣,李俊峰,何新林

(1.石河子大學(xué)水利建筑工程學(xué)院 石河子 832000;2.中國(guó)環(huán)境科學(xué)研究院 北京 100012;3.中國(guó)科學(xué)院地理科學(xué)與資源研究所 北京 100101;4.中國(guó)科學(xué)院大學(xué) 北京 100049)

農(nóng)田土壤水是限制農(nóng)作物生長(zhǎng)必需的因素之一,缺少有效的灌溉制度導(dǎo)致的土壤干旱,易導(dǎo)致根區(qū)環(huán)境惡化、作物吸水困難甚至死亡,而適當(dāng)灌溉可在一定程度上減輕缺水對(duì)作物生長(zhǎng)造成的負(fù)面影響[1-5]。目前灌溉農(nóng)業(yè)消耗著世界上70%的淡水用水量,我國(guó)農(nóng)業(yè)灌溉水平相比其他發(fā)達(dá)國(guó)家較為落后,其中華北地區(qū)灌溉水利用率僅為0.587[6-10]。華北平原農(nóng)田灌溉用水多來(lái)自于地下水,對(duì)地下水使用與管理造成較大壓力且地下水大量開(kāi)采造成的環(huán)境問(wèn)題十分突出。

目前水分運(yùn)移預(yù)測(cè)常用的Gardner模型和van Genuchten模型具有較好的模擬效果,后者因更高的精確度被廣泛使用[11-14]。HYDRUS-1D模型由美國(guó)鹽土實(shí)驗(yàn)室研發(fā),其中內(nèi)嵌van Genuchten模型,可用于分析各種飽和多孔介質(zhì)中一維垂向的水流和溶質(zhì)運(yùn)移,亦可用于模擬水、熱和多種溶質(zhì)在不同飽和介質(zhì)中的運(yùn)動(dòng)[15-17]。該模型提供了豐富的參數(shù)設(shè)置及邊界條件,可根據(jù)實(shí)際情況進(jìn)行選擇。已有學(xué)者使用HYDRUS-1D模型模擬土壤水鹽運(yùn)移狀態(tài)后發(fā)現(xiàn),鹽分隨灌水而向下層運(yùn)動(dòng),且鹽分含量隨灌水定額的增加而減小[18-20];而潘延鑫等[21]使用HYDRUS模型在陜西鹵泊灘進(jìn)行試驗(yàn)并發(fā)現(xiàn),500 m3·hm?2的灌溉定額可以有效控制當(dāng)?shù)赝寥婪e鹽。王志坤[22]使用HYDRUS-1D模型對(duì)相同灌水定額、不同灌水模式條件下1 m土體內(nèi)鹽分含量做出對(duì)比,發(fā)現(xiàn)溝灌模式下土壤平均含鹽量低于畦灌模式;采用不同礦化度的灌水水質(zhì)可能會(huì)使作物嚴(yán)重減產(chǎn),但土壤積鹽率會(huì)隨灌水礦化度的增加而增加[23-24]。González等[25]模擬出充分及不充分灌溉下土壤水分運(yùn)動(dòng)情況,其均方根誤差為0.014～0.025 cm3·cm?3。探索農(nóng)田水鹽運(yùn)移規(guī)律,對(duì)幫助作物生長(zhǎng)及緩解水資源危機(jī)均具有重要意義。因此,使用HYDRUS-1D模型探究和模擬灌溉條件下華北平原地區(qū)農(nóng)田土壤水鹽運(yùn)移規(guī)律,對(duì)于提高灌溉水利用效率,緩解該區(qū)域農(nóng)田灌溉用水壓力具有重要意義。

我國(guó)HYDRUS-1D模型的應(yīng)用研究多集中于西北干旱地區(qū),且多數(shù)模擬試驗(yàn)于室外農(nóng)田開(kāi)展,灌水定額難以精確控制,而針對(duì)華北平原小麥(Triticum aestivum)根區(qū)范圍內(nèi)土壤水鹽運(yùn)移模擬研究較為空缺。本研究以華北地區(qū)中國(guó)科學(xué)院禹城綜合試驗(yàn)站為試驗(yàn)點(diǎn),選擇當(dāng)?shù)氐湫投←溵r(nóng)田土壤為研究對(duì)象,采用野外觀測(cè)和室內(nèi)土柱試驗(yàn)相結(jié)合的方法,借助HYDRUS-1D軟件建立水鹽運(yùn)移模型,明確模擬軟件對(duì)研究區(qū)水鹽運(yùn)移研究的適配性,探究華北平原農(nóng)田水鹽運(yùn)移的規(guī)律,為此模型驗(yàn)證提供數(shù)據(jù)支撐,并在優(yōu)化灌溉管理制度方面建立一定理論基礎(chǔ)。

1 研究區(qū)概況及研究方法

1.1 研究區(qū)概況

試驗(yàn)在山東省禹城市中國(guó)科學(xué)院禹城綜合試驗(yàn)站開(kāi)展,禹城市地處華北平原中部(36°41′36″～37°12′13″N,116°22′11″～116°45′00″E),屬暖溫帶大陸季風(fēng)氣候,年均氣溫13.3 ℃,年均降水量555.5 mm,其中6—9月總降水量約占全年降水量的76%,月際分布不均勻[26]。境內(nèi)主要地表徑流為徒駭河、趙牛河及潘莊引黃干渠,主要耕作方式為冬小麥-夏玉米(Zea mays)輪作,持續(xù)的土壤改良使得全區(qū)糧食產(chǎn)量從1980年的1.16×105t增加到2015年的8.28×105t,成為華北平原地區(qū)典型的產(chǎn)糧地[27]。

1.2 水分運(yùn)移模擬研究

研究利用HYDRUS-1D模型模擬不同深度土壤體積含水率實(shí)測(cè)值與模擬值,闡明土壤含水率變化規(guī)律,探究土壤水分變動(dòng)的驅(qū)動(dòng)因素。選用均方根誤差(RMSE)、平均誤差(AE)和納什效率系數(shù)(NSE)評(píng)價(jià)HYDRUS-1D模型對(duì)研究區(qū)水分運(yùn)移模擬的適配性并分析模型出現(xiàn)偏差的原因,其中NSE高于0.75時(shí)一般認(rèn)為模擬結(jié)果良好[28]。

1.2.1 土壤含水率實(shí)測(cè)數(shù)據(jù)獲取

選取試驗(yàn)站內(nèi)長(zhǎng)期種植的冬小麥試驗(yàn)田作為試驗(yàn)地,采用土壤水分廓線系統(tǒng)(Sentek Pty Ltd.,澳大利亞)測(cè)定土壤體積含水率,使用CR200X數(shù)據(jù)采集器(Campbell Scientific,美國(guó))采集并記錄數(shù)據(jù)。土壤水分監(jiān)測(cè)系統(tǒng)埋設(shè)深度分別為10 cm、20 cm、30 cm、40 cm和60 cm,頻率為每30 min采集數(shù)據(jù)一次。

1.2.2 土壤水分運(yùn)移的基本方程

土壤體積含水率的計(jì)算采用一維Richards方程:

式中:θ為土壤體積含水率(cm3·cm?3);t為時(shí)間(d);K為非飽和土壤導(dǎo)水率(cm·d?1);h為壓力水頭(cm);α為水流方向和垂直方向的夾角(°),在垂向運(yùn)動(dòng)中α=0°;S表示作物對(duì)土壤水的吸收作用(d?1);z表示垂直高程(cm),以地表為0且向上為正方向。本次研究不考慮作物吸收,故S=0。

選用van Genuchten-Mualem模型擬合公式(1)中涉及參數(shù):

式中:m、n均為土壤水分保持函數(shù)中的參數(shù),m=1?1/n,n＞1;sθ為土壤飽和含水率(cm3·cm?3);rθ為土壤殘余含水率(cm3·cm?3);Ks為飽和導(dǎo)水率(cm·d?1);α代表進(jìn)氣吸力的倒數(shù)(cm?1);l代表空隙連通性,此值一般取0.5[29];Se代表相對(duì)飽和度。

1.2.3 模型參數(shù)設(shè)置

鑒于冬小麥的主根區(qū)為地表以下40 cm[30],故本次水分運(yùn)移模擬范圍選取地表以下80 cm,空間步長(zhǎng)取1 cm,于10 cm、20 cm、30 cm、40 cm和60 cm處添加觀測(cè)點(diǎn)。模擬時(shí)長(zhǎng)50 d,初始時(shí)間步長(zhǎng)取0.05 d,最小時(shí)間步長(zhǎng)和最大時(shí)間步長(zhǎng)分別取10?5d和1 d[31],根據(jù)實(shí)測(cè)數(shù)據(jù),上、下邊界初始含水率分別設(shè)置為0.2 cm3·cm?3、0.48 cm3·cm?3。由于研究區(qū)地下水埋深較大,多年浮動(dòng)于4～5 m間[32],對(duì)模擬土層影響可忽略不計(jì),故分別選擇可積水的大氣邊界、自由排水邊界為上、下邊界建立模型。運(yùn)動(dòng)方程中其余參數(shù)根據(jù)模型土壤數(shù)據(jù)庫(kù)中自帶12種土壤特性參數(shù),及試驗(yàn)地土壤容重和土壤顆粒組成由模型預(yù)測(cè)得到,經(jīng)與實(shí)測(cè)值校準(zhǔn)優(yōu)化得到參數(shù),如表1所示。土壤容重及飽和體積含水率采用環(huán)刀法測(cè)定,殘余體積含水率采用烘干法稱重測(cè)定。

1.3 鹽分運(yùn)移模擬研究

為更精確地控制灌溉水量,此次鹽分運(yùn)移模擬選擇于室內(nèi)開(kāi)展土柱試驗(yàn),設(shè)置定額灌水量,對(duì)比不同時(shí)刻土壤剖面含鹽量實(shí)測(cè)值與模擬值,觀測(cè)土壤剖面含鹽量隨時(shí)間的變化規(guī)律及分布特征;為探索不同灌水量對(duì)土壤鹽分運(yùn)移規(guī)律的影響,進(jìn)行3次試驗(yàn),設(shè)置多種灌水量:8 L、16 L、24 L(對(duì)應(yīng)農(nóng)田灌溉量6.37 mm、12.73 mm、19.10 mm),驗(yàn)證HYDRUS-1D模型對(duì)研究區(qū)鹽分運(yùn)移模擬的可行性并討論分析模型出現(xiàn)偏差的原因。

表1 試驗(yàn)地土壤基本物理性質(zhì)和水力學(xué)參數(shù)Table 1 Basic physical properties and hydraulic parameters of the experimental soil

1.3.1 土壤中鹽分運(yùn)移的基本方程

本次研究選擇Equilibrium模型。溶質(zhì)運(yùn)動(dòng)方程如下:

式中:c為土壤溶液中溶質(zhì)濃度(mg·cm?3);q為水流通量(cm·d?1);D為水動(dòng)力彌散系數(shù)(cm2·d?1);θ為土壤體積含水率(cm3·cm?3);t為時(shí)間(d);z代表垂直高程(cm),以地表為0且向上為正方向。

1.3.2 模型參數(shù)設(shè)置

選擇上邊界為濃度通量邊界,下邊界為零濃度梯度邊界,鹽分含量以液相濃度計(jì)算。土壤中溶質(zhì)的運(yùn)移主要取決于垂向彌散度的大小,通常室外田間試驗(yàn)的垂向彌散度可取5～20 cm,本次模擬選擇取10 cm[33],總模擬時(shí)長(zhǎng)取48 h,根據(jù)實(shí)測(cè)數(shù)據(jù),設(shè)置土壤剖面初始含鹽量為0.17 g·kg?1。

1.3.3 土柱試驗(yàn)

本次土壤鹽分運(yùn)移研究制作土柱高Z=100 cm、直徑D=40 cm,在底層鋪放濾砂,防止堵塞排水口,按照1.33 g·cm?3容重分層回填土壤,并在土壤表層鋪放碎石,防止灌水時(shí)水流過(guò)激,灌溉用水取自研究區(qū)地下水,礦化度為0.9 g·L?1。分別在距離土體表面10 cm、20 cm、40 cm、60 cm處開(kāi)直徑為1 cm的孔用以埋設(shè)傳感器并觀測(cè)數(shù)據(jù)。采用CS655水溫鹽傳感器(Campbell Scientific,美國(guó))、CR1000數(shù)據(jù)采集器(Campbell Scientific,美國(guó))監(jiān)測(cè)土壤溫度(測(cè)量范圍:?10～70 ℃;精度:±0.5 ℃)、電導(dǎo)率(測(cè)量范圍:0～8 dS·m?1;精度:±5%)及含水率(測(cè)量范圍:5%～50%;精度:±3%),裝置示意圖如圖1所示。

此次試驗(yàn)將待測(cè)土樣按照土水比1∶5浸提并過(guò)濾,視濾液電導(dǎo)率為土壤電導(dǎo)率(EC),含鹽量(Q)使用蒸干法測(cè)定,并如圖2所示擬合土壤電導(dǎo)率-含鹽量標(biāo)準(zhǔn)曲線,得出Q=0.7775EC+0.0066,R2=0.8977,相關(guān)性較高。后可根據(jù)擬合公式由土壤電導(dǎo)率計(jì)算出土壤含鹽量[34]。

2 結(jié)果與分析

2.1 土壤水分運(yùn)移模擬

根據(jù)所建模型,得到土壤不同深度處含水率隨時(shí)間變動(dòng)曲線如圖3所示?？梢钥闯鰷\層(10 cm)土壤體積含水率變化最為劇烈,體積含水率在0.5 d內(nèi)從0.23 cm3·cm?3增加到0.42 cm3·cm?3。而深層(60 cm)土壤體積含水率在2.3 d內(nèi)僅從0.41 cm3·cm?3增加到0.44 cm3·cm?3。

接受補(bǔ)給后,土壤各含水率隨深度依次升高,且補(bǔ)給對(duì)深層土壤的影響遠(yuǎn)小于對(duì)淺層土壤。這說(shuō)明受到蒸散發(fā)作用、淺層土壤及作物吸收、土壤導(dǎo)水率和吸水性能的限制,外界補(bǔ)給對(duì)深層土壤含水率影響很小。土壤體積含水率越大,吸力值越小,當(dāng)土壤含水率趨于飽和時(shí),吸力值趨于零。

如表2所示,水分運(yùn)移模擬結(jié)果的均方根誤差平均值為0.0176 cm3·cm?3,納什效率系數(shù)平均值為0.826,變異系數(shù)為0.0560,表示此模型能較好且較穩(wěn)定地模擬研究區(qū)土壤內(nèi)部水分的運(yùn)移過(guò)程,可作為研究區(qū)水分運(yùn)移理論研究的有效工具。

表2 不同土層深度土壤水分含量模擬結(jié)果驗(yàn)證Table 2 Verification of simulation results of soil water contents at different soil depths

2.2 土壤鹽分運(yùn)移模擬

根據(jù)所建模型得到土壤鹽分運(yùn)移結(jié)果如圖4所示,12 h時(shí)含鹽量峰值處于14 cm處,為0.715 g·kg?1;24 h時(shí)峰值位置下移至18 cm處,為0.649 g·kg?1;48 h時(shí)峰值位置下移至23 cm處,為0.584 g·kg?1,較36 h前降低18.3%,各時(shí)刻鹽分含量在垂向上整體呈現(xiàn)先增大后減少的分布規(guī)律,峰值所處位置的下移速度及峰值的減小速度均逐漸降低。

各時(shí)刻土壤鹽分運(yùn)移模擬結(jié)果驗(yàn)證如表3所示,鹽分運(yùn)移模擬結(jié)果的均方根誤差平均值為0.099 g·kg?1,平均誤差平均值為0.174 g·kg?1,但納什效率系數(shù)的變異系數(shù)達(dá)到9.71,表明此模型在部分時(shí)刻能較好地模擬研究區(qū)土壤鹽分的運(yùn)移過(guò)程。

3種灌水量入滲48 h后土壤剖面鹽分分布結(jié)果如圖5所示,并如圖6對(duì)實(shí)測(cè)值和模擬值進(jìn)行相關(guān)性分析,得到?jīng)Q定系數(shù)R2為0.8216,相關(guān)性較高。當(dāng)灌水量增加至16 L、24 L,入滲48 h后,其含鹽量峰值分別升高為0.720 g·kg?1、0.786 g·kg?1,相較于灌水量8 L時(shí)分別增加23.3%和34.6%;峰值所處位置分別下移至26 cm、29 cm深度處,說(shuō)明灌溉水量的增加會(huì)顯著加劇土壤中鹽分的累積效應(yīng)。

表3 各入滲時(shí)間后不同深度土壤鹽分含量模擬結(jié)果驗(yàn)證Table 3 Verification of simulation results of salt contents at different soil depths after each infiltration time

3 討論

此次研究結(jié)果證實(shí),HYDRUS-1D模型可較為精確地模擬研究區(qū)農(nóng)田土壤剖面水分及鹽分運(yùn)移,5種深度的土壤含水率在接受補(bǔ)給后的整體變動(dòng)趨勢(shì)與前人研究所得結(jié)論基本一致[19,25],其中淺層土壤水鹽模擬結(jié)果由于受到強(qiáng)烈的氣象因素及補(bǔ)給影響,與實(shí)測(cè)值偏差較為明顯,隨著深度增加,外界因素的影響逐漸弱化,模擬結(jié)果趨于穩(wěn)定[35]。但土柱試驗(yàn)所得鹽分分布規(guī)律與羅朋等[36]的研究結(jié)論有所不同,據(jù)其研究所示,入滲后土壤剖面鹽分的分布呈現(xiàn)持續(xù)增長(zhǎng)的趨勢(shì),即下層土壤的含鹽量大于上層土壤。但此次研究表明,入滲后土壤剖面含鹽量應(yīng)隨深度呈現(xiàn)先增長(zhǎng)后減少的規(guī)律,即在剖面某深層處出現(xiàn)鹽分含量的峰值,且灌水量越大,同一深度處鹽含量越大,含鹽量峰值越大,峰值所處位置越深。鑒于土壤水分運(yùn)移會(huì)攜帶鹽分,故本研究將此次模型計(jì)算所得水通量與含鹽量變化整理如圖7所示,結(jié)合鹽分運(yùn)移結(jié)果可以看出,土壤剖面內(nèi)水流通量隨時(shí)間逐漸減少,鹽分位移也隨之減緩,含鹽量累積逐漸達(dá)到最大值。這是因?yàn)楫?dāng)補(bǔ)給遠(yuǎn)大于蒸發(fā)時(shí),重力成為水鹽運(yùn)移的主導(dǎo)因素,鹽分跟隨水分下移,且隨著水分運(yùn)移的穩(wěn)定而逐漸停止,即會(huì)在土壤濕潤(rùn)鋒附近出現(xiàn)鹽分累積的現(xiàn)象,鄭冬梅[37]的研究結(jié)果與這一推測(cè)相符。

本研究中,補(bǔ)給初期以及土柱淺層的土壤鹽分運(yùn)移模擬結(jié)果與實(shí)測(cè)值偏差較大。這與淺層土壤接受補(bǔ)給后內(nèi)部水鹽運(yùn)移過(guò)于劇烈有關(guān),導(dǎo)致模型難以精確計(jì)算;同時(shí)由于表層土壤受耕作、根系等人類活動(dòng)影響,理化性質(zhì)變異性較大,導(dǎo)致模擬值與實(shí)測(cè)值偏差增大[26]。另一導(dǎo)致模型出現(xiàn)誤差的原因筆者分析為土壤環(huán)境溫度變動(dòng),故繪制土柱內(nèi)各深度處實(shí)測(cè)土壤含鹽量與實(shí)測(cè)土壤溫度隨時(shí)間的變化趨勢(shì)圖(圖8),對(duì)比發(fā)現(xiàn)土壤各深度含鹽量會(huì)受溫度影響變化且整體與溫度呈正相關(guān)。此現(xiàn)象與樊會(huì)敏等[38]的研究結(jié)論相吻合。出現(xiàn)這種現(xiàn)象是因?yàn)?土壤接受補(bǔ)給后,濕潤(rùn)鋒下移并連通上下層土壤毛管水,每日9:00后氣溫開(kāi)始爬升,下層土壤溫度升高,蒸發(fā)導(dǎo)致下層土壤水分?jǐn)y帶著鹽分依靠毛管力上升并滯留于上層[39],故相對(duì)于蒸發(fā)強(qiáng)度較弱時(shí),在土壤溫度變化劇烈的時(shí)段內(nèi),模型模擬結(jié)果會(huì)與實(shí)測(cè)值產(chǎn)生更大的偏差。

小麥需水量最大的時(shí)間段在拔節(jié)期到開(kāi)花期(從3月底至4月底),日均需水量約5.71 mm;越冬期到返青期(11月底至翌年3月初)需水量較小,日均需水量不到0.25 mm,而小麥在萌芽期對(duì)土壤含鹽量較為敏感,其他時(shí)期耐鹽性較好,但土壤鹽分含量仍不得高于6 dS·m?1,否則會(huì)導(dǎo)致作物蒸騰量減少,吸水困難[40-41]。針對(duì)此類問(wèn)題,可借助HYDRUS模型,通過(guò)改變初始條件、豐富模塊參數(shù),準(zhǔn)確且快速地預(yù)測(cè)土壤中水分及鹽分的動(dòng)態(tài)變化,有效避免不合理灌溉所帶來(lái)的水資源浪費(fèi)及作物根區(qū)環(huán)境惡化等問(wèn)題。

4 結(jié)論

以中國(guó)科學(xué)院禹城綜合試驗(yàn)站為研究區(qū),結(jié)合野外長(zhǎng)期觀測(cè)所得土壤水分運(yùn)移數(shù)據(jù)及室內(nèi)土柱試驗(yàn)所得鹽分運(yùn)移數(shù)據(jù),利用HYDRUS-1D模型研究該地區(qū)農(nóng)田水鹽運(yùn)移規(guī)律。得出以下結(jié)論:

1)實(shí)地監(jiān)測(cè)的土壤初始體積含水率隨土壤深度的增加而增加,發(fā)生補(bǔ)給時(shí),土壤體積含水率隨深度陸續(xù)增加,其中淺層土壤的波動(dòng)最為劇烈,但受到淺層土壤及作物的吸收、蒸散發(fā)作用和吸水性能的限制,補(bǔ)給對(duì)研究區(qū)深層(60 cm)土壤含水率影響較小。

2)借助模型對(duì)研究區(qū)土壤水鹽運(yùn)移時(shí)空分布規(guī)律進(jìn)行模擬,所得結(jié)果與實(shí)測(cè)值吻合度較高,表明模型具有較好的仿真性,可用作研究區(qū)水鹽運(yùn)移理論研究的一種有效手段。

3)增加灌水量導(dǎo)致土柱土壤內(nèi)鹽分向下位移,表現(xiàn)出顯著的農(nóng)田洗鹽效果。發(fā)生補(bǔ)給時(shí),土壤鹽分隨水分下移,且隨時(shí)間的推移,水分通量的減少對(duì)鹽分持續(xù)運(yùn)移的影響顯著,導(dǎo)致土壤剖面在不同時(shí)刻存在不同深度的積鹽區(qū)。