基于數(shù)值模擬的地下水位變化下灌水頻次方案研究

林鵬飛，王水獻，蔣夢瑤，楊 麗

(蘭州大學(xué)資源環(huán)境學(xué)院，蘭州 730000)

土壤水，作為水循環(huán)過程中的一個重要環(huán)節(jié)，是地表水、地下水和大氣水的重要聯(lián)系紐帶。在植物及農(nóng)作物的生長過程中，土壤水同樣扮演著重要的角色，因此對土壤水的研究也就顯得尤為重要[1,2]。

由于土壤水的不可見性使其觀測較困難，研究中通常需借助模型對土壤水的入滲過程進行模擬。HYDRUS-1D作為用于模擬變飽和多孔介質(zhì)中水分、能量、溶質(zhì)運移的數(shù)值模型[3]，便可以很好地模擬土壤水分在土壤中的分布，時空變化，運移規(guī)律，分析人們普遍關(guān)注的農(nóng)田灌溉[4-6]、田間施肥[7-9]、環(huán)境污染[10]等實際問題。如卞建民[11]等在野外調(diào)查和資料收集的基礎(chǔ)上,利用HYDRUS-1D模型模擬分析了旱田(玉米地)灌溉條件下地下水入滲補給過程，并發(fā)現(xiàn)旱田灌溉條件下，土壤水與地下水會發(fā)生雙向的、動態(tài)的水量頻繁交換。Xuezhi Tan[12]等基于交替濕潤和干燥(AWD)灌溉和持續(xù)灌溉(CFI)設(shè)計了8種不同試驗情形，運用HYDRUS-1D進行模擬，發(fā)現(xiàn)運用HYDRUS-1D正向?qū)С龅膲毫λ^模擬結(jié)果并不理想，而反向模擬出的結(jié)果較好，同時得出了交替濕潤和干燥(AWD)灌溉對比持續(xù)灌溉(CFI)而言滲漏量更低的結(jié)論。

本試驗將研究區(qū)選在了焉耆盆地下五號渠村，已有研究[13]分析了焉耆盆地灌水量與灌水頻次對籽瓜性狀影響，但尚未有研究以土壤水的波動性和土壤水的下滲量作為聯(lián)合評價指標，分析在地下水位和灌水頻次正交改變條件下干旱綠洲灌區(qū)的最優(yōu)灌溉方案。因此，本文便以新疆焉耆盆地下五號渠村作為典型干旱綠洲灌區(qū)，結(jié)合HYDRUS-1D模型模擬，綜合考慮灌水頻次和地下水位變化對土壤水分規(guī)律的影響，分析并得到了滴灌條件下該區(qū)最優(yōu)灌水方案。

1 材料與方法

1.1 研究區(qū)概況

焉耆盆地位于新疆塔里木盆地東北側(cè)，屬于南天山山脈中的中生代斷陷盆地，是一個較為典型的綠洲-荒漠交錯地區(qū)。該地位于東經(jīng) 86°39′～88°20′，北緯 41°23′～43°31′，總面積約為1.3 萬km2(其中灌區(qū)面積約為0.56 萬km2)，其地勢西高東低，北高南低，東西方向長約170 km，南北方向?qū)捈s60～90 km，總體表現(xiàn)為四周向盆地傾斜的地貌形態(tài)。地面海拔一般在1 050～1 200 m，最低處為我國最大的內(nèi)陸淡水湖——博斯騰湖。焉耆盆地在氣候上屬于暖溫帶大陸性干旱氣候，熱量和光照豐富，多年平均氣溫為8.6 ℃，年降水量為50～80 mm，年蒸發(fā)能力為2 000～2 450 mm[14-16]。

1.2 數(shù)據(jù)來源

本次研究主要獲取了焉耆盆地下五號渠村土壤水分監(jiān)測數(shù)據(jù)，此外還使用了焉耆縣氣象數(shù)據(jù)，主要包括焉耆氣象站實測蒸散發(fā)數(shù)據(jù)、最低和最高氣溫、平均氣溫、相對濕度、風(fēng)速、日照時數(shù)和降雨。根據(jù)焉耆盆地番茄灌溉制度得到了模擬期作物灌溉數(shù)據(jù)。為了掌握試驗區(qū)土壤的分層情況，獲取了焉耆盆地下五號渠村土壤巖性數(shù)據(jù)。為了確定建立模型的邊界條件，獲取了試驗區(qū)地下水位變化數(shù)據(jù)。

1.3 蒸散發(fā)的計算

根據(jù)1998年聯(lián)合國糧農(nóng)組織(FAO)推薦的Penman-Monteith(PM)公式[17]，計算得出了研究區(qū)參考作物蒸發(fā)蒸騰量，計算公式如下：

(1)

式中：ET0為參考作物蒸發(fā)蒸騰量，mm/d；Δ為飽和水汽壓溫度曲線上的斜率，kPa/K；Rn為凈太陽輻射，J/(m2·d)；G為土壤熱通量，J/(m2·d)；ea為飽和水汽壓，kPa；ed為實際水汽壓，kPa；γ為濕度計常數(shù)，kPa/K；U2為距地表2 m處的風(fēng)速，m/s。

為了驗證本次研究在HYDRUS-1D模型模擬中選取Penman-Monteith(PM)計算蒸散發(fā)的準確性，將1998年聯(lián)合國糧農(nóng)組織(FAO)推薦的Penman-Monteith(PM)公式計算的參考作物蒸騰蒸發(fā)量與氣象站實際觀測的大型蒸發(fā)皿蒸發(fā)量進行對比，對比結(jié)果圖如圖1所示。

圖1 參考作物蒸騰蒸發(fā)量計算值與實測值

計算相對誤差為-0.053 7，計算獲得的參考作物蒸騰蒸發(fā)量變化趨勢和實際蒸騰蒸發(fā)量變化趨勢也大致相同，可見模型中選取Penman-Monteith(PM)公式計算蒸散發(fā)準確度較高。

1.4 試驗設(shè)計

試驗設(shè)計的原方案為試驗時間從2010年7月3日至8月24日，平均地下水位1 051.07 m，總灌水量為930 m3/hm2，灌水頻次為3次，灌水時間為7月20日、7月26日、8月7日，3次灌水量分別為360、330、240 m3/hm2。

結(jié)合試驗設(shè)計的原方案，運用HYDRUS-1D進行數(shù)值模擬，以驗證模型模擬的精度，并在此基礎(chǔ)上設(shè)計了改變地下水位、改變灌水頻次(方案一、方案二、方案三)、綜合考慮改變地下水位和灌水頻次的試驗方案。

2 HYDRUS模型

2.1 模型的建立

模型基本方程：采用Richards方程(1931年)作為水流運動的控制性方程，表征土壤水力特性采用van Genuchten 方程(1980年)。

試驗時間從2010年7月3日至8月24日，模擬時間共計53 d，模擬時段為番茄成熟期，采用變時間步長剖分方式，根據(jù)收斂迭代次數(shù)調(diào)整時間步長，允許最小時間步長為0.001 d，最大時間步長為1 d。

模型模擬土層深度為250 cm，利用HYDRUS-1D模型將土壤剖面離散化分為251個節(jié)點，每個節(jié)點間距1 cm。再根據(jù)試驗地實測的土壤質(zhì)地巖性數(shù)據(jù)，將土壤劃分為沙壤土(0～65 cm)，壤土(65～85 cm)，粉壤土(85～130 cm)，黏土(130～250 cm)4層。通過4層不同土壤巖性賦予不同的土壤水力參數(shù)值。

2.2 臨界條件

初始條件：根據(jù)試驗區(qū)數(shù)據(jù)，初始邊界條件確定為使用觀測的土壤水分含量數(shù)據(jù)。

上邊界條件：由于建立的模型上邊界直接與大氣接觸，且試驗區(qū)基本不產(chǎn)生徑流，故上邊界確定為大氣邊界。

下邊界條件：試驗區(qū)為淺埋深區(qū)，平均埋深1.93 m，模擬時段地下水埋深均小于2.5 m，而模擬土層深度為2.5 m，模擬時間內(nèi)地下水水位隨時間變化，故下邊界條件確定為變水頭邊界。

2.3 參數(shù)敏感性分析

HYDRUS-1D中采用的參數(shù)調(diào)整主要是針對土壤水力特性參數(shù)的調(diào)整，包括對土壤殘余含水量Qr，土壤飽和含水量Qs，土壤空氣進氣值倒數(shù)α，土壤孔徑指數(shù)n，土壤飽和導(dǎo)水率Ks的調(diào)整。在參數(shù)調(diào)整過程中，主要針對改變上述參數(shù)中某一參數(shù)，參數(shù)的變化ΔX/X取增加5%、減少5%、增加10%、減少10%，然后計算模擬結(jié)果的變化ΔY/Y，通過二者的比值獲得參數(shù)敏感性系數(shù)S[18,19]。計算公式如下：

(2)

式中：Y為模型結(jié)果值；X為選取的參數(shù)值；ΔY和ΔX均為兩者的改變量。最后在率定過程中主要調(diào)對應(yīng)參數(shù)敏感性系數(shù)S大的參數(shù)，本次研究中取調(diào)整S大于0.2的土壤水力參數(shù)。

最終通過敏感性系數(shù)計算，發(fā)現(xiàn)30 cm土層深度土壤水分含量主要對土壤水力參數(shù)Qs1、α1、n1敏感，50 cm土層深度土壤水分含量主要對土壤水力參數(shù)Qs1、α1、n1敏感，80 cm土層深度土壤水分含量主要對土壤水力參數(shù)Qs2、α2、n2敏感，120 cm土層深度土壤水分含量主要對土壤水力參數(shù)Qs3、n3敏感。故在模型校正中，主要調(diào)整參數(shù)Qs、α、n來提高模型模擬的精度(下標1、2、3、4分別代表砂壤土層、壤土層、粉壤土層、黏土層4個土層)。

2.4 誤差分析

采用均方根誤差(Root Mean Square Error，RMSE)和相對誤差(Relative Error，RE)兩個評價指標來評估模擬精度，計算公式如下:

(3)

(4)

式中：si和oi分別為第i個樣本的模擬值和實測值；N為樣本個數(shù)；均方根誤差(RMSE)的值和相對誤差(RE)的絕對值越小，則反映模擬精度越好。

3 結(jié)果和分析

3.1 模型的驗證和校正

通過調(diào)參，最終得到的土壤水分含量模擬結(jié)果如圖2所示。

圖2 不同深度土壤水分含量的實測值與模擬值對比圖

由圖2不同深度土壤水分含量的實測值與模擬值對比中可以看出，30、80、120 cm的土壤水分含量的對比值均較好的貼合在1:1線兩側(cè)，而50 cm土壤水分含量的對比值中顯示實測值大部分大于模擬值，但總體上模擬和實測的對比值均貼合1∶1線。

再運用誤差分析定量分析模型的模擬精度，得出在30 cm土層深度的土壤水分含量實測值和模擬值的RE為0.037，RMSE為0.020，50 cm土層深度的土壤水分含量實測值和模擬值的RE為-0.074，RMSE為0.028，80 cm土層深度的土壤水分含量實測值和模擬值的RE為-0.042，RMSE為0.021，120 cm土層深度的土壤水分含量實測值和模擬值的RE為0.025，RMSE為0.033。從誤差分析的結(jié)果上看，各土層深度的模擬結(jié)果均較好，HYDRUS-1D的模擬結(jié)果可適用于試驗區(qū)域土壤水分規(guī)律的模擬。

3.2 不同地下水位下土壤水分變化

在淺埋深地區(qū)，地下水位的變化往往對土壤水分的變化有著顯著影響，從而影響到灌溉區(qū)灌溉作物的生長。為此，本文在對原方案進行模型模擬的基礎(chǔ)上研究了地下水位變化對不同深度土壤水分的影響，設(shè)計了如下方案：地下水位上升0.2 m，地下水位上升0.5 m，地下水位下降0.2 m，地下水位下降0.5 m。發(fā)現(xiàn)當(dāng)?shù)叵滤簧仙龝r，每一天的不同深度土壤水分均呈上升趨勢，當(dāng)?shù)叵滤幌陆禃r，每一天的不同深度土壤水分均呈下降趨勢。為了更好地定量了解水位變化對土壤水分含量變化的影響，對不同地下水位變化方案下的土壤水分含量取平均值進行比較，結(jié)果如圖3所示。

圖3 地下水位變化條件下的土壤水分含量

由圖3數(shù)據(jù)可以發(fā)現(xiàn)，當(dāng)?shù)叵滤幌陆禃r，不同深度土壤水分含量變化雖有所下降但下降幅度均不明顯，可以認為地下水位下降對土壤水分含量變化無影響，如地下水位下降0.5 m，30 cm土層深度的土壤水分含量降幅0.001 5，50 cm土層深度的土壤水分含量降幅0.001 9，80 cm土層深度的土壤水分含量降幅0.002，120 cm土層深度的土壤水分含量降幅0.010。當(dāng)?shù)叵滤簧仙龝r，不同深度土壤水分含量均有一定程度的上升，且不同深度的土壤水分含量增幅有所不同，如地下水位上升0.5 m，30 cm土層深度的土壤水分含量增幅0.011 6，50 cm土層深度的土壤水分含量增幅0.016 8，80 cm土層深度的土壤水分含量增幅0.018，120 cm土層深度的土壤水分含量增幅0.064 8，因此，土層深度越深，土壤水分含量增幅越大。從理論上分析原因，可能是在淺埋深地區(qū)，當(dāng)?shù)叵滤幌陆禃r，模擬的土層深度和地下水距離較遠，受地下水影響小，而當(dāng)?shù)叵滤簧仙龝r，模擬的土層深度與地下水距離十分接近，會得到地下水補給且對土壤水分含量影響相對較大。因此，在淺埋深地區(qū)不同土層深度的土壤水分含量受地下水位變化的影響程度與其接觸地下水的距離有關(guān)。

3.3 不同灌水頻次下土壤水分變化

為了探求試驗區(qū)在同等水量下最優(yōu)灌水方案，在原方案的基礎(chǔ)上，首先僅改變灌水頻次，設(shè)計了如下3個灌水頻次方案。

方案一：灌水頻次為5次，灌水時間為7月3日、7月20日、7月26日、8月7日、8月20日，每次灌水186 m3/hm2。

方案二：灌水頻次為7次，灌水時間為7月3日、7月10日、7月20日、7月26日、8月7日、8月14日、8月20日，每次灌水132.85 m3/hm2。

方案三：灌水頻次為9次，灌水時間為7月3日、7月9日、7月15日、7月21日、7月27日、8月2日、8月8日、8月14日、8月20日，每次灌水103.33 m3/hm2。

將3種灌水頻次方案下的模擬結(jié)果與原方案的模擬結(jié)果進行對比，得出各模擬方案下土壤水分含量的標準偏差見表1。

土壤水分含量的標準偏差反映土壤水分的波動情況，標準偏差越小，說明土壤水分維持越穩(wěn)定，土壤水分波動越小，對作物生長也就越有利，由表1中數(shù)據(jù)可以得出，在淺層土壤(30 cm和50 cm土層深度，也是根系層深度活動范圍，是影響作物生長的主要土層深度)，隨著灌水頻次的增加，土壤水分波動性越小。在深層土壤，灌水頻次的增加對土壤水分的波動性影響不大。綜上所述，在僅考慮土壤水分波動性情況下，以方案三最優(yōu)，方案二次之，方案一與原方案不推薦采用。

表1 不同灌水頻次方案下土壤水分含量的標準偏差

3.4 不同灌水頻次下土壤水分下滲量

底部邊界通量是單位時間內(nèi)通過單位面積的水量，包括地下水向上補給量和土壤水的深層滲漏量[20],底部邊界通量的大小可以反映土壤水下滲情況，對于確定合理的灌溉方案有一定的借鑒作用。本文便對不同灌水頻次下底部邊界通量的變化進行研究，得到了不同方案下底部邊界通量的變化如圖4所示。

圖4 不同方案下底部邊界通量日變化和地下水埋深

由圖4可以看出，不同灌水頻次方案下底部邊界通量的變化趨勢是一致的，即改變灌水頻次不影響底部邊界通量的變化趨勢。而影響土壤水深層滲漏和地下水向上補給的主要因素為地下水埋深，當(dāng)?shù)叵滤裆钭儨\時，則地下水開始向上補給；當(dāng)?shù)叵滤裆钭兩顣r，則主要是土壤水向下層滲漏。綜上所述，影響土壤水下滲補給趨勢的主要因素為地下水埋深，與試驗中的灌水頻次無關(guān)。

再分析不同灌水頻次方案下的底部邊界通量的總量，得到的結(jié)果見表2(表2和表4結(jié)果方向均向下，即所得結(jié)果均為向下滲漏)。

由此可見，在方案二下總底部邊界通量最小，對應(yīng)的下滲量最小。故而在僅考慮下滲的情況下方案二水分利用效率更高，為最優(yōu)灌水頻次方案。再綜合考慮土壤水分變化和土壤水分下滲，最優(yōu)灌水頻次方案為方案二。

表2 不同灌水頻次方案下的下滲情況

3.5 不同灌水頻次和地下水位下的土壤水分及下滲量

不同灌水頻次和不同地下水位都會對土壤水分的變化都會產(chǎn)生影響，進而對作物生長產(chǎn)生影響，故而在設(shè)計試驗方案時需綜合考慮二者的影響，本次試驗在確定了最優(yōu)灌水頻次方案的基礎(chǔ)上設(shè)置不同地下水位方案，綜合二者正交方案，以確定最終最優(yōu)方案。由于地下水位上升會導(dǎo)致土壤鹽分的增加，故在方案設(shè)計中僅考慮地下水位下降的情況，在方案二地下水位不變方案的基礎(chǔ)上，設(shè)計了地下水位下降0.2 m，地下水位下降0.5 m兩種方案。同樣地，以土壤水的波動性和土壤水的下滲量作為評價指標，得出了正交方案下的土壤水分含量標準偏差和下滲情況見表3和表4。

表3 地下水位變化下土壤水分含量的標準偏差

表4 地下水位變化下的下滲情況

由標準偏差數(shù)據(jù)可以發(fā)現(xiàn)，隨著地下水位的降低，不同深度的土壤水的波動性開始增加，故而在地下水位不變的時候土壤水分的波動性最小。對比不同地下水位的下滲情況，發(fā)現(xiàn)隨著地下水位的降低，土壤水的下滲量也逐漸增加，地下水位不變的方案下總下滲量最小。綜上所述，在綜合考慮了灌水頻次和地下水位變化的情況下，本次研究推薦的最優(yōu)灌水方案為地下水位不變時的方案二。

4 結(jié) 語

(1)HYDRUS-1D可以很好地模擬干旱綠洲灌區(qū)的土壤水分運移規(guī)律。從模擬結(jié)果圖上看，各土層深度的土壤水分含量模擬和實測的對比值均貼合1∶1線。從誤差分析結(jié)果上看，各土層深度的模擬結(jié)果誤差都較小，模擬結(jié)果良好。

(2)在淺埋深地區(qū)，不同深度的土壤水分含量會隨地下水位的上升而上升，同時，不同深度土壤水分含量受地下水位變化的影響程度與其接觸地下水的距離有關(guān)，與地下水的距離越遠，受地下水位影響越小。

(3)不同灌水頻次下土壤水分的波動情況是不同的，灌水頻次的變化主要影響淺層土壤水分的波動性，且隨著灌水頻次的增加土壤水分波動性減小。

(4)不同灌水頻次底部邊界通量的大小也是不同的，但底部邊界通量的變化趨勢不受灌水頻次的影響，主要受地下水埋深變化的影響，通過比較不同灌水頻次下的下滲量的大小，發(fā)現(xiàn)方案二下的土壤水下滲量最小，為12.357 cm/d。

(5)綜合考慮地下水位和灌水頻次的變化，設(shè)計了正交方案，以土壤水分波動情況和土壤水的下滲情況作為評價指標，本試驗推薦采用的最佳灌水方案為地下水位不變，平均地下水位1 051.07 m，灌水7次，每次灌水132.85 m3/hm2。