王 鍇, 楊澤元, 袁 悅, 陳志軍

(1.長(zhǎng)安大學(xué) 環(huán)境科學(xué)與工程學(xué)院, 陜西 西安 710054; 2.旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室,陜西 西安 710054; 3.陜西省地下水與生態(tài)環(huán)境工程研究中心, 陜西 西安 710054)

毛烏素沙地位于半干旱氣候區(qū)，境內(nèi)包含陜北和內(nèi)蒙能源化工基地，隨著該地區(qū)工農(nóng)業(yè)發(fā)展與人口增長(zhǎng)，用水量短缺形勢(shì)越發(fā)嚴(yán)峻。該地區(qū)主要用水來(lái)源為地下水，降水入滲是地下水補(bǔ)給的主要來(lái)源，研究該地區(qū)降水入滲補(bǔ)給是其地下水資源評(píng)價(jià)的重要基礎(chǔ)工作[1]。前人[2-3]常采用年降水量與區(qū)域入滲補(bǔ)給系數(shù)計(jì)算區(qū)域性降水入滲補(bǔ)給量，很少考慮入滲過(guò)程的滯后性。地下水補(bǔ)給總是滯后于降水事件，某一時(shí)段內(nèi)補(bǔ)給量實(shí)際由該時(shí)段之前的時(shí)段內(nèi)降水所決定。同時(shí)，由于補(bǔ)給滯后的存在，含水層對(duì)于人類活動(dòng)與氣候變化的響應(yīng)往往不會(huì)立即發(fā)生[4]，這給地下水可持續(xù)開(kāi)發(fā)帶來(lái)挑戰(zhàn)。因此研究地下水的補(bǔ)給滯后現(xiàn)象對(duì)于提升地下水資源評(píng)價(jià)精度與地下水資源合理利用具有重要意義。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于陜西省榆林市補(bǔ)浪河鄉(xiāng)，地處毛烏素沙地南緣。年均降水量365 mm，年均蒸發(fā)量1 246 mm。主要潛水含水層為第四系風(fēng)積孔隙含水層，巖性為中砂，不含黏土(表1)。地下水常年埋深1.3～1.5 m，主要補(bǔ)給來(lái)源為大氣降水，占天然補(bǔ)給總量85%以上[16]，主要排泄去路為蒸發(fā)排泄。

表1 研究區(qū)不同層土壤質(zhì)地

1.2 研究方法

在原國(guó)土資源部地下水與生態(tài)陜西榆林野外觀測(cè)基地建立原位試驗(yàn)場(chǎng)。監(jiān)測(cè)地下5，10，30，60，90 cm深處含水率，同時(shí)獲取場(chǎng)地氣象數(shù)據(jù)。包氣帶含水率和溫度由5 TM傳感器監(jiān)測(cè)(美國(guó)Decagon公司)，采集頻率為10 min。氣象數(shù)據(jù)由試驗(yàn)場(chǎng)內(nèi)的波紋比系統(tǒng)(05130-5RMYong，美國(guó)Campbell公司)自動(dòng)監(jiān)測(cè)。地下水位由MiniDiver自動(dòng)監(jiān)測(cè)并記錄數(shù)據(jù)，氣象數(shù)據(jù)和地下水位的采集頻率均為1 h。

以野外原位監(jiān)測(cè)數(shù)據(jù)為基礎(chǔ)研究場(chǎng)地內(nèi)包氣帶含水率與地下水位對(duì)降水入滲的滯后響應(yīng)。根據(jù)場(chǎng)地條件建立水文地質(zhì)概念模型并利用Hydrus-1D軟件建立降水入滲模型。采用野外監(jiān)測(cè)數(shù)據(jù)進(jìn)行模型參數(shù)的識(shí)別與優(yōu)化。

采用相對(duì)平均誤差(AVRE)、相關(guān)系數(shù)(R)作為指標(biāo)進(jìn)行識(shí)別期模型誤差分析。

AVRE=∑(VC1/VO1)/n

(1)

(2)

式中：VCi——模擬值； VOi——原位觀測(cè)值；n——數(shù)據(jù)總數(shù)量。

通過(guò)上述過(guò)程識(shí)別獲取基本參數(shù)，基于原位監(jiān)測(cè)含水率，負(fù)壓，溫度數(shù)據(jù)建立反解模型(共14 d，4 695個(gè)數(shù)據(jù)點(diǎn))。反解計(jì)算采用Hydrus-1D的Inverse solution模塊進(jìn)行，對(duì)van-Genuchten方程中的各參數(shù)θr(殘余含水率)，α,n(水土特征曲線方程的參數(shù))和Ks(飽和滲透系數(shù))進(jìn)行反解，反解方法為最小目標(biāo)函數(shù)法[17]。選取SSR，r2最優(yōu)對(duì)應(yīng)的參數(shù)作為參數(shù)優(yōu)化結(jié)果。其中SSR為反解過(guò)程中目標(biāo)函數(shù)值。r2為擬合值與觀測(cè)值的相關(guān)系數(shù)，r2為1表示擬合最佳。

(3)

注意到,熱帶西北太平洋位于東亞季風(fēng)區(qū),而熱帶東南印度洋位于南亞季風(fēng)區(qū)以南(Wang and Linho,2002),環(huán)流結(jié)構(gòu)有所不同。同時(shí),在氣候?qū)W上夏季位于對(duì)流層低層起自索馬里至西太平洋熱帶地區(qū)的西風(fēng)氣流和對(duì)流層上層的東風(fēng)氣流將西太平洋熱帶地區(qū)與熱帶印度洋地區(qū)聯(lián)系起來(lái)。對(duì)印度洋與太平洋上環(huán)流變化聯(lián)系的研究?jī)H限于赤道或近赤道的緯向區(qū)域,其聯(lián)系的重要紐帶是Walker環(huán)流,而對(duì)于赤道外地區(qū)兩大洋上空環(huán)流變化的聯(lián)系卻未能有較好的研究。那么,位于海洋性大陸東南邊緣和西北邊緣的兩個(gè)地區(qū)上的環(huán)流變化是否存在聯(lián)系呢?這種聯(lián)系是否存在年代際變化呢?其原因是什么呢?本文將回答這些問(wèn)題。

參數(shù)獲取完畢后，以累計(jì)底部通量為指標(biāo)，采用單參數(shù)擾動(dòng)方法進(jìn)行敏感性分析，重點(diǎn)關(guān)注不同參數(shù)改變對(duì)于地下水補(bǔ)給滯后的影響。采用輸入輸出變化率(ROV)來(lái)分析模型計(jì)算對(duì)于參數(shù)的敏感性[15]。

(4)

式中：C(t)——模型實(shí)際計(jì)算結(jié)果(改變某一參數(shù)值的結(jié)果)；Cref(t)——模型參考計(jì)算結(jié)果(所有參數(shù)保持不變的結(jié)果)，此處計(jì)算結(jié)果均為累計(jì)底部通量；P——實(shí)際輸入?yún)?shù)；Pref為參考輸入?yún)?shù)；t——時(shí)間。ROV(t)的絕對(duì)值越大，則模型對(duì)該參數(shù)越敏感。由于介質(zhì)被劃為多層，因此輸入輸出變化率均采用各層平均值。

2 結(jié)果與分析

2.1 包氣帶含水率響應(yīng)滯后

研究區(qū)降水年際變化大，降水主要集中在每年的6—9月份，為豐雨期。由于每年的11月至次年3月是該區(qū)的冰凍期，只發(fā)生少量降雪，研究時(shí)段確定為2014年3—10月并將時(shí)段內(nèi)的降水依據(jù)國(guó)標(biāo)劃分為小雨、中雨、大雨與暴雨[18]。包氣帶深度越淺，降水后含水率響應(yīng)越劇烈，部分降水事件中地表最大含水率甚至接近飽和；同一地下水埋深下，降水量越大，包氣帶響應(yīng)深度越大。以8月26日至9月5日為例(圖1)，8月27日降水0.82 cm(小雨)，僅5 cm深度有明顯響應(yīng)，其他深度含水率均無(wú)明顯變化；8月30日降水7.1 cm(暴雨)，包氣帶所有深度含水率都明顯上升。選取研究時(shí)段內(nèi)主要降水事件，通過(guò)相關(guān)性分析研究降水量與入滲響應(yīng)深度的關(guān)系。如圖1所示，入滲響應(yīng)深度與降水量呈顯著線性相關(guān)，相關(guān)系數(shù)R為0.94>R(16，0.01)=0.59。隨著包氣帶深度增加，含水率響應(yīng)存在明顯滯后；以7月8日降水為例，10 cm處滯后時(shí)間為13 h，30 cm為19 h，90 cm為29 h。

圖1 入滲響應(yīng)深度與降水量相關(guān)關(guān)系

注意到7月3日降水(降水量2.12 cm)歷時(shí)4 h，最大雨強(qiáng)14 mm/h，5，10，30 cm深度含水率均明顯響應(yīng)，但90 cm處含水率只微弱變化；7月8日(降水量4.1 cm)歷時(shí)27 h，最大雨強(qiáng)僅3 mm/h，90 cm處含水率明顯上升。前者降水強(qiáng)度更大，后者降水量更大，這說(shuō)明降水量對(duì)于響應(yīng)深度的影響大于降水強(qiáng)度(圖2)。

圖2 降水量與含水率變化時(shí)間序列

2.2 地下水滯后補(bǔ)給

無(wú)降水的時(shí)段，地下水位保持相對(duì)穩(wěn)定(圖3)，小雨和中雨后地下水位均未響應(yīng)，如6月25日降水5.2 mm，水位并無(wú)變化。大雨暴雨均會(huì)對(duì)地下水形成有效補(bǔ)給(如6月29日、7月3日、7月9日的降水)，使得地下水位上升。試驗(yàn)時(shí)段內(nèi)對(duì)地下水形成有效補(bǔ)給的降水統(tǒng)計(jì)結(jié)果詳見(jiàn)表2。從表2可知，降水量越大，地下水位上升幅度越大。相近降水事件如7月3日降水(降水量2.12 cm)與6月29日降水(共2.24 cm)下，水位埋深與地下水位升幅反相關(guān)。這是由于在蒸發(fā)條件相似時(shí)，水位埋深越大，包氣帶水分虧缺越大，降水后補(bǔ)給量也越小。降水之后，地下水位并不立即上升，一般把從降水開(kāi)始到地下水位上升的時(shí)間間隔稱為補(bǔ)給滯后時(shí)間；由表2可知，6月29日降水(共2.24 cm)補(bǔ)給滯后時(shí)間為4 h，7月3日降水(共2.12 cm)補(bǔ)給滯后時(shí)間為5 h；注意到這一滯后時(shí)間甚至小于30 cm處含水率響應(yīng)時(shí)間，表明降水后存在優(yōu)勢(shì)流通道，降水入滲通過(guò)優(yōu)勢(shì)流通道提前進(jìn)入地下水，引起地下水位響應(yīng)。歷次降水后地下水補(bǔ)給滯后時(shí)間均隨著降水強(qiáng)度的增大而減小，但減小幅度逐漸放緩。這是由于隨著降水強(qiáng)度的增大，入滲率也隨之增大；但當(dāng)降水強(qiáng)度逐漸接近或超過(guò)土壤入滲能力后，入滲率不再變化，相應(yīng)的入滲補(bǔ)給滯后時(shí)間也趨于一定值。需要注意的是，9月10日降水3.06 cm，較9月22日降水5.46 cm差別并不大，但后者地下水位上升劇烈，這是由于后者的前期降水量更大，前期降水有效補(bǔ)給了包氣帶水分，較高的包氣帶含水率促進(jìn)了后期降水入滲。

圖3 地下水位與降水量變化時(shí)間序列

雨 型大 雨暴 雨降水量/mm21.241.430.671.654.622.4降水強(qiáng)度/(mm·h-1)4.241.531.717.92.14.48地下水位升幅/cm1.31.62.99.19.91.8補(bǔ)給滯后時(shí)間/h5118364 地下水埋深/cm雨前129.4129.5130.6145.3133.1121.6雨后128.1127.9127.7136.2123.2119.8

2.3 降水入滲數(shù)值模型建立與參數(shù)敏感性分析

2.3.1 模型建立 以地表為坐標(biāo)原點(diǎn)，向上為正。建立水文地質(zhì)概念模型，地表作為模型的上邊界。地下埋深150 cm處作為模型的下邊界，忽略水平徑流，由此可以概化為變飽和帶一維非穩(wěn)定流模型。含水介質(zhì)依據(jù)顆粒分析和容重概化為3層，其中0—20 cm為第1層，20—80 cm為第2層，80—150 cm為第3層。在Hydrus中選用大氣邊界作為水分運(yùn)移的上邊界條件，可變壓力水頭作為水分運(yùn)移的下邊界條件。熱運(yùn)移的上下邊界條件分別為上下邊界溫度。選用Van-Genuchten方程作為水力特征曲線，考慮持水曲線的滯后效應(yīng)。含水率初始條件采用與初始時(shí)間地下水位相平衡的含水率分布。溫度初始條件設(shè)置為初始時(shí)刻監(jiān)測(cè)值。模型的初始土壤水力學(xué)參數(shù)詳見(jiàn)表3。采用Hydrus中基于顆粒分析和容重?cái)?shù)據(jù)的神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)結(jié)果。

表3 經(jīng)神經(jīng)元預(yù)測(cè)分析得到的土壤水力學(xué)參數(shù)

注：θs，α，n,，Ks，θr均為van-Genuchten方程中參數(shù)。θr為殘余含水率；α和n為水土特征曲線方程的參數(shù)；Ks為飽和滲透系數(shù)；θs為飽和含水率。

2.3.2 模型參數(shù)獲取 選取2014年7月2—15日作為模型參數(shù)的識(shí)別期。識(shí)別期包氣帶深度5，30，60 cm處含水率、溫度模擬值與觀測(cè)值對(duì)比如圖4所示。

圖4 模型識(shí)別期包氣帶不同深處含水率與溫度

模型在識(shí)別期的計(jì)算誤差結(jié)果詳見(jiàn)表4—5，除埋深60 cm處誤差較大外，總體計(jì)算結(jié)果較好。經(jīng)過(guò)參數(shù)反演優(yōu)化，選取SSR=1.854，r2=0.926時(shí)對(duì)應(yīng)的參數(shù)作為優(yōu)化結(jié)果，最終確定的水、熱力學(xué)參數(shù)詳見(jiàn)表6—7。

表4 識(shí)別期模型計(jì)算含水率誤差結(jié)果

注： AVRE為模型計(jì)算值與監(jiān)測(cè)值的相對(duì)平均誤差，R為相關(guān)系數(shù)。

表5 識(shí)別期模型計(jì)算溫度誤差結(jié)果

2.3.3 van-Genuchten方程參數(shù)敏感性分析 由于調(diào)參過(guò)程中涉及的參數(shù)較多，對(duì)van-Genuchten主要的土壤水力學(xué)參數(shù)θr，α，n，Ks進(jìn)行了敏感性分析。通過(guò)對(duì)某一參數(shù)進(jìn)行增減，維持其余參數(shù)不變，確定參數(shù)敏感性[15]。敏感性分析中參數(shù)取用情況詳見(jiàn)表8。

表6 模型識(shí)別及優(yōu)化后所得土壤水力學(xué)參數(shù)

表7 模型識(shí)別及優(yōu)化后所得土壤熱力學(xué)參數(shù)

注：b1,b2,b3為熱傳導(dǎo)系數(shù)λ(θ)的參數(shù)；Cn為固態(tài)物質(zhì)的體積熱容；Cw為水的體積熱容；Di為熱擴(kuò)散系數(shù)；So為固態(tài)物質(zhì)體積百分比。

表8 參數(shù)敏感性分析中所用土壤水力學(xué)參數(shù)

注：考慮到n最小取1，故n的變化為在原始值與1作差的基礎(chǔ)上再增減25%。

采用累積底部通量誤差進(jìn)行參數(shù)敏感性分析。時(shí)間起點(diǎn)為7月2日0：00，圖中下標(biāo)為0的表示該參數(shù)的初始值。-25表示減小25%的參數(shù)值，+25表示增加25%的參數(shù)值(圖5)。各參數(shù)的輸入輸出變化率詳見(jiàn)表9。從表9可知，底部通量對(duì)θs最為敏感，對(duì)Ks最不敏感。α與n的敏感性介于θs與Ks之間，其中n比α更敏感。若以底部通量代替地下水補(bǔ)給量，則地下水補(bǔ)給量的大小主要由θs控制，θs越大，入滲過(guò)程中包氣帶持水能力越強(qiáng)，地下水補(bǔ)給量越小。從圖5中可以看出，θs與Ks對(duì)于降水補(bǔ)給的滯后影響較大，當(dāng)Ks減小或θs增大時(shí)，累計(jì)底部通量拐點(diǎn)時(shí)間(指示著補(bǔ)給開(kāi)始時(shí)間)明顯滯后；其中θs增大25%，滯后時(shí)間增長(zhǎng)10～16 h；Ks減小25%，滯后時(shí)間增長(zhǎng)6～8 h。值得注意的是，若包氣帶剖面處于蒸發(fā)狀態(tài)，Ks的改變幾乎不引起底部通量的變化，此時(shí)Ks的敏感性很弱。

圖5 不同參數(shù)變化對(duì)模擬底部通量的影響

項(xiàng)目α-25n?-25Ks-25θs-25α+25n?+25Ks+25θs+25ROV(tend)1.222.930.442.21.783.760.254.44ROV(tmax)9341550241517136018401322630ROV(tmin)1.212.880.310.00751.782.110.00134.42

注：tend表示結(jié)束時(shí)刻，tmax為ROV(t)最大時(shí)刻，tmin為ROV(t)最小時(shí)刻； 下標(biāo)為-25表示在原參數(shù)基礎(chǔ)上減小25%的參數(shù)值，+25表示增加25%的參數(shù)值。

3 討 論

3.1 降水入滲滯后補(bǔ)給的影響因素探討

由于土地利用變化、氣候變化和工農(nóng)業(yè)發(fā)展，世界不同地區(qū)的土壤水分和地下水正越來(lái)越多地受到影響但呈現(xiàn)出不同規(guī)律。本文確定的地下水滯后補(bǔ)給時(shí)間為4～11 h，小于相近水位埋深下三江平原地區(qū)滯后時(shí)間(28～50 d)[19]，這是由于后者包氣帶巖性更細(xì)，入滲速率較慢。在大埋深地區(qū)，由于入滲補(bǔ)給路徑很長(zhǎng)，地下水補(bǔ)給滯后時(shí)間甚至能達(dá)數(shù)十至數(shù)百年[4，20]。土壤水分分布同樣是影響入滲補(bǔ)給的重要因素，按照活塞流入滲理論，包氣帶含水率響應(yīng)滯后時(shí)間與土壤水分分布的關(guān)系可以表示為[21]：

I(t)=(θ0-θi)Zf

(5)

式中：I(t)——累計(jì)入滲量；Zf——入滲深度；θ0——入滲面上部含水率；θi——剖面初始含水率。

由公式(5)得知，包氣帶含水率響應(yīng)滯后時(shí)間隨包氣帶初始含水率增大而減小，而初始含水率反映了前期降水的影響,因此入滲補(bǔ)給除受到本次降水的影響，還受到前期降水的控制。

3.2 多參數(shù)敏感性分析

本研究采用的是單參數(shù)敏感性計(jì)算，能夠較為精確得出各個(gè)參數(shù)的敏感性。但實(shí)際的土壤參數(shù)之間存在相關(guān)性且總是處于一個(gè)合理范圍。比如θs增大導(dǎo)致降水入滲滯后時(shí)間增大，但是θs增大同樣會(huì)導(dǎo)致Ks增大，后者能促進(jìn)降水入滲；因此，對(duì)于Ks與θs等相關(guān)性較強(qiáng)的參數(shù)采用多因素分析能更為準(zhǔn)確地探討其敏感性。

3.3 未來(lái)研究

由于研究區(qū)的地下水埋深淺，降水入滲補(bǔ)給速率相對(duì)較快，包氣帶淺部補(bǔ)給滯后不明顯，通過(guò)監(jiān)測(cè)數(shù)據(jù)確定補(bǔ)給滯后時(shí)間的精度并不高。此外，研究時(shí)段內(nèi)降水頻次高。因此后續(xù)研究應(yīng)當(dāng)考慮累計(jì)降水對(duì)于降水入滲的影響并提高淺部包氣帶監(jiān)測(cè)精度。研究區(qū)處在固定沙丘上，地表平緩，側(cè)向徑流微弱，故采用一維滲流模型；但場(chǎng)地內(nèi)稀疏分布著沙柳，干旱時(shí)期沙柳吸收水分較強(qiáng)，對(duì)土壤水分分布存在影響[22]。同時(shí)地下水補(bǔ)給滯后時(shí)間小于包氣帶含水率響應(yīng)滯后時(shí)間，這表明場(chǎng)地中存在不可忽視的優(yōu)勢(shì)流通道，可能的原因是場(chǎng)地內(nèi)包氣帶巖性不均一或植被根系吸水[23]。因此后續(xù)研究也應(yīng)當(dāng)針對(duì)植被吸水做進(jìn)一步分析并考慮使用雙滲透模型解決優(yōu)勢(shì)流問(wèn)題。

4 結(jié) 論

結(jié)合原位監(jiān)測(cè)數(shù)據(jù)，分析了毛烏素沙地地下水淺埋區(qū)降水入滲補(bǔ)給規(guī)律，利用Hydrus-1D軟件建立了降水入滲數(shù)值模型，獲取了模型水力學(xué)參數(shù)和熱力學(xué)參數(shù)；分析了van-Genuchten方程主要水力學(xué)參數(shù)敏感性；為研究包氣帶和地下水位對(duì)降水入滲的響應(yīng)機(jī)制、分析降水入滲補(bǔ)給機(jī)理奠定了基礎(chǔ)。

(1) 包氣帶含水率響應(yīng)深度隨著降水量增大而增大，其中小雨型響應(yīng)深度為3～10 cm，中雨型為30～60 cm，大雨型為60～90 cm，暴雨型均大于90 cm。入滲響應(yīng)深度與降水量的關(guān)系可表示為：D=35.5P-17.1，其中D為入滲響應(yīng)深度(cm)，P為降水量大小。隨著前期累計(jì)降水量增大，降水入滲深度也增大。

(2) 研究區(qū)內(nèi)地下水補(bǔ)給滯后時(shí)間約為4～11 h。包氣帶水分和地下水滯后補(bǔ)給受地下水埋深、包氣帶巖性和含水率分布共同影響，其中包氣帶水分滯后補(bǔ)給時(shí)間隨著初始含水率增大而減小。

(3) 建立了該地變飽和帶水汽熱耦合運(yùn)移模型。采用模型識(shí)別與驗(yàn)證獲取了包氣帶水力學(xué)參數(shù)與熱力學(xué)參數(shù)，以此作為參數(shù)反演優(yōu)化的初始輸入值；參數(shù)優(yōu)化最佳結(jié)果θr為0.03，θs為0.34，Ks為23 cm/h,各層n為2.0～2.2，α為0.060～0.074。

(4) 與數(shù)值模型有關(guān)的各項(xiàng)水力學(xué)參數(shù)中，敏感性自強(qiáng)到弱的排序是θs，n，α，KS，底部通量對(duì)于θs最為敏感；降水補(bǔ)給滯后時(shí)長(zhǎng)與θs呈正相關(guān)關(guān)系。