王強(qiáng)民,董書寧,王文科,王 皓

(1.煤炭科學(xué)研究總院，北京 100013; 2.中煤科工集團(tuán)西安研究院有限公司，陜西 西安 710054; 3.陜西省煤礦水害防治技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054; 4.長安大學(xué) 環(huán)境科學(xué)與工程學(xué)院，陜西 西安 710054; 5.旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054)

榆神礦區(qū)是我國重要的煤炭生產(chǎn)基地和高強(qiáng)度開采區(qū)，同時也是最為典型的生態(tài)環(huán)境脆弱區(qū)[1]，采礦活動、植被恢復(fù)及氣候變化控制著礦區(qū)的整體生態(tài)格局[2-4]。受采煤地下水位下降、采煤塌陷及礦區(qū)土地類型轉(zhuǎn)換等因素影響，煤炭開采過程中出現(xiàn)了植被退化的現(xiàn)象，為避免植被的退化及死亡，礦區(qū)開展了高強(qiáng)度的植被恢復(fù)工程。研究表明，截至2006年底，神東礦區(qū)累計(jì)投入生態(tài)環(huán)境建設(shè)經(jīng)費(fèi)2.69億元，使得由開發(fā)初期(1985年)植被覆蓋率3%～11%增長至59.4%以上(2005年)[5]，礦區(qū)和井田尺度的植被覆蓋度大幅提升[6-8]，植被對降水的截留量和蒸騰量也隨之增大[9-10]，造成降水對地下水的補(bǔ)給量減小。而研究區(qū)的主要地貌類型為風(fēng)積沙，透水性好，降雨入滲補(bǔ)給量占地下水總補(bǔ)給量的78.5%[11]，且區(qū)內(nèi)河流40%～70%的流量來自地下水的補(bǔ)給[12]，水文生態(tài)環(huán)境與地下水的關(guān)系密切。因此，分析高強(qiáng)度植被恢復(fù)對地下水補(bǔ)給的影響已成為生態(tài)脆弱礦區(qū)亟待解決的關(guān)鍵科學(xué)問題之一。

目前，在生態(tài)脆弱礦區(qū)地下水和植被關(guān)系的研究方面，主要集中在地下水狀態(tài)變化對植被的影響[13-14]。例如通過開展野外原位監(jiān)測試驗(yàn)，研究了采煤塌陷區(qū)土壤水的變異規(guī)律及生態(tài)效應(yīng)[13];結(jié)合原位試驗(yàn)和數(shù)值模擬，分析了植被對采煤地下水位下降的響應(yīng)機(jī)制[14]。近年來，針對植被變化對地下水的影響也開展了部分研究[9,15-16]，但研究整體偏少。例如，張俊等通過建立區(qū)域尺度剖面二維模型，分析了鄂爾多斯盆地植被變化對地下水流系統(tǒng)的影響[15];張建軍研究了黃土高原水土保持對土壤水分的影響，認(rèn)為人工刺槐的耗水量大，甚至大于當(dāng)?shù)亟邓?，?dǎo)致深層土壤干化[16]，地下水很難得到補(bǔ)給[9]。但是，以上研究主要集中在黃土高原區(qū)，而榆神礦區(qū)植被和地下水的相互反饋?zhàn)饔幂^黃土高原區(qū)強(qiáng)烈，研究結(jié)果不能很好的指導(dǎo)生態(tài)脆弱礦區(qū)植被恢復(fù)實(shí)踐。

筆者結(jié)合野外裸土區(qū)和植被覆蓋區(qū)2個原位試驗(yàn)點(diǎn)的觀測數(shù)據(jù)，以2種模式下的降水—土壤水—地下水轉(zhuǎn)化過程為研究對象。通過遙感解譯分析礦區(qū)的生態(tài)環(huán)境變化情況;采用原位監(jiān)測手段分析植被對土壤水、地下水位的影響規(guī)律;在此基礎(chǔ)上，利用數(shù)值仿真計(jì)算分析植被覆蓋度增加對降水補(bǔ)給地下水的影響機(jī)制。研究成果對豐富我國西部生態(tài)脆弱礦區(qū)植被恢復(fù)與地下水資源的協(xié)調(diào)發(fā)展具有重要意義。

1 研究區(qū)概況

榆神礦區(qū)位于陜西省神木縣榆陽區(qū)境內(nèi)，煤炭資源儲量豐富、地質(zhì)構(gòu)造簡單、開采技術(shù)條件優(yōu)越，其延安組含煤地層自上而下賦存6層可采煤層。主采2-2煤層，平均埋深261.15 m，可采厚度0.80～12.49 m，平均厚度達(dá)到10.04 m[17]。受隔水層性質(zhì)、煤層厚度和埋深等因素影響，采動裂隙發(fā)育至基巖含水層內(nèi)部甚至近地表的松散層含水層，改變了地下水和土壤水的運(yùn)移規(guī)律，表現(xiàn)為地下水位下降和土壤含水率虧缺，繼而加劇了礦區(qū)水資源短缺和生態(tài)環(huán)境脆弱的現(xiàn)狀[18]。近年來，為了落實(shí)綠色礦山建設(shè)要求，各大煤礦開展了廣泛的礦山植被恢復(fù)工程，礦區(qū)植被覆蓋度得到整體提升。為評估礦區(qū)生態(tài)恢復(fù)植被覆蓋度增加對地下水補(bǔ)給的影響，研究團(tuán)隊(duì)在榆神礦區(qū)無定河流域(距榆林市40 km)搭建了原位試驗(yàn)平臺，試驗(yàn)區(qū)多年平均降水量為400 mm，平均潛在蒸散發(fā)量為1 153 mm，地貌以風(fēng)積沙為主，沙柳、沙蒿為主要植被類型。

2 研究方法

2.1 原位試驗(yàn)

本次試驗(yàn)在原位試驗(yàn)平臺的植被耗水區(qū)進(jìn)行，該區(qū)設(shè)置玻璃鋼材質(zhì)的圓柱形監(jiān)測井2個，分別為裸土區(qū)和植被覆蓋區(qū)2個原位試驗(yàn)監(jiān)測井。監(jiān)測井高度為4 m，直徑為2 m，側(cè)向和底部為隔水邊界，上部接受降水補(bǔ)給和蒸散發(fā)排泄。為合理控制和觀測地下水位，監(jiān)測井中間設(shè)置地下水注(抽)水孔和地下水位觀測孔。植被選取研究區(qū)的優(yōu)勢植被沙柳作為研究對象，初始地下水位埋深設(shè)置為3 m(極限蒸發(fā)深度)，試驗(yàn)介質(zhì)采用礦區(qū)典型的風(fēng)積沙。含水率測定采用Decagon公司生產(chǎn)的ECH2O土壤水分傳感器，地下水位測定采用斯倫貝謝水務(wù)公司生產(chǎn)的Mini-Diver傳感器測定，測定頻率均為1 h/次。裸土區(qū)含水率傳感器在垂直方向上的布置為:3,10,20,30,50,80,150,250,350 cm，由于填土均勻，裸土區(qū)水平方向不設(shè)置對照傳感器。植被覆蓋區(qū)含水率傳感器在垂直方向的布置和裸土區(qū)一致，考慮植被耗水對水平方向土壤水的影響，水平方向設(shè)置對照傳感器，分別距離植被中心點(diǎn)0,20,50,80 cm(圖1)。

圖1 原位試驗(yàn)示意(裸土區(qū)+植被覆蓋區(qū))

2.2 數(shù)值模擬

在原位試驗(yàn)的基礎(chǔ)上，筆者采用美國國家鹽改中心開發(fā)的Hydrus軟件[19]模擬裸土和不同植被覆蓋條件下的降水—土壤水—地下水相互轉(zhuǎn)化過程，并求解植被耗水量和降水入滲系數(shù)。植被覆蓋條件下，一維非飽和土壤水運(yùn)動的控制方程為

(1)

式中，C(h)為容水度，cm-1;h為土壤水壓力水頭，cm;t為時間，d;z為垂向坐標(biāo)，cm，向上為正;K(h)為非飽和滲透系數(shù)，cm/d;S(z,t)為植被根系吸水速率,裸土條件下，S(z,t)為0。

文中采用van-Genuchten-Mualem模型[20]描述土壤水分特征曲線和滲透系數(shù)曲線:

(2)

(3)

(4)

式中，θs為飽和含水量，cm3/cm3;θr為殘余含水量，cm3/cm3;Ks為飽和滲透系數(shù)，cm/d;Se為有效飽和度；α,m,n為相關(guān)土壤參數(shù)，其中m=1-1/n，n>1;l為彎曲度參數(shù)。

筆者采用Feddes模型[21]描述植被的根系吸水過程:

S(z,t)=α(h)b(z)Tp

(5)

式中，α(h)為水分脅迫函數(shù);b(z)為根系分布函數(shù);Tp為植被潛在蒸騰量，cm/d,可由下式計(jì)算:

Tp=ETp(1-e-kLAI)

(6)

式中，ETp為潛在蒸散發(fā)量，cm/d;k為冠層的消光系數(shù);LAI為葉面積指數(shù)。

模型構(gòu)建過程中的初始條件(土壤水?dāng)?shù)據(jù))、上邊界條件(氣象數(shù)據(jù))、下邊界條件(地下水位數(shù)據(jù))來源于本次原位試驗(yàn)的實(shí)時監(jiān)測，土壤水分特征參數(shù)和植被根系密度函數(shù)見文獻(xiàn)[3]。在模擬計(jì)算之前，本文根據(jù)原位試驗(yàn)的觀測結(jié)果，對相關(guān)參數(shù)進(jìn)行了識別和驗(yàn)證，30,50和80 cm處的土壤水的觀測值和模擬值動態(tài)變化趨勢一致(圖2)。利用式(7)，(8)求得相對誤差和均方根誤差，見表1。

(7)

(8)

式中，RE為相對誤差，%;RMSE為均方根誤差，cm3/cm3;n為模擬天數(shù);Si和Oi分別為第i天的模擬值和實(shí)測值，cm3/cm3。

圖2 模型識別和驗(yàn)證過程土壤水的觀測值和模擬值(2016年)

表1 模型識別和驗(yàn)證過程中的誤差分析

由表1可知模型識別和驗(yàn)證過程中土壤水模擬值與實(shí)測值的相對誤差<10%、均方根誤差<0.01 cm3/cm3，反映了基于野外觀測數(shù)據(jù)構(gòu)建的數(shù)值模型可以較好的分析研究區(qū)裸土和植被覆蓋條件下的降水—土壤水—地下水相互轉(zhuǎn)化過程。

3 結(jié)果分析

3.1 礦區(qū)生態(tài)環(huán)境演變特征

文中選用數(shù)據(jù)質(zhì)量高、信息量大的MODIS數(shù)據(jù)產(chǎn)品研究榆神礦區(qū)生態(tài)環(huán)境演變特征，其時間分辨率為16 d，空間分辨率為250 m。為增強(qiáng)數(shù)據(jù)的可靠程度，分別選取研究區(qū)8月份植被指數(shù)(NDVI)和植被凈初級生產(chǎn)力(NPP)兩個指標(biāo)聯(lián)合反映礦區(qū)生態(tài)環(huán)境演變特征，其中NPP值采用光能利用率模型(CASA)估算[22]。研究結(jié)果顯示，礦區(qū)尺度的NDVI值和NPP值都呈持續(xù)增加現(xiàn)象，NDVI值從2000年的0.247增長到2015年的0.358，NPP值從2000年的12.67 g/m2增長到2015年的20.78 g/m2(圖3)，升高幅度分別為45%和64%，表現(xiàn)為低覆蓋區(qū)連續(xù)減小，中、高(較高)覆蓋區(qū)持續(xù)增大(表2)，反映礦區(qū)生態(tài)環(huán)境持續(xù)好轉(zhuǎn)。

圖3 榆神礦區(qū)平均NPP和NDVI變化特征

3.2 植被影響下的土壤水變化規(guī)律

根據(jù)原位試驗(yàn)觀測結(jié)果，分別統(tǒng)計(jì)相同地下水埋深和氣象條件下裸土區(qū)和植被覆蓋區(qū)0～80 cm內(nèi)的土壤水監(jiān)測數(shù)據(jù)，并繪制裸土區(qū)(一維)和植被覆蓋區(qū)(二維)月平均土壤含水率分布特征圖(圖4)。裸土區(qū)平均土壤含水率為0.101 cm3/cm3，受氣象條件和地下水位的雙重影響，含水率隨埋深的增大呈增加趨勢;受植被根系耗水影響，植被覆蓋區(qū)0～80 cm內(nèi)的平均土壤含水率僅為0.062 cm3/cm3，較裸土區(qū)減少38.6%。

表2 不同植被覆蓋度條件下的面積分布比例

圖4 裸土區(qū)和植被覆蓋區(qū)平均土壤含水率分布特征

為更好的分析不同條件下土壤水分的運(yùn)移規(guī)律，利用風(fēng)積沙的土壤水分特征曲線，將裸土區(qū)(一維)和植被覆蓋區(qū)(二維)0～80 cm內(nèi)土壤含水率換算成土壤基質(zhì)勢，利用式(9)計(jì)算土壤總水勢。

H=h(θ)+Z

(9)

式中，H為土壤總水勢，cm;h(θ)為土壤基質(zhì)勢，cm;θ為土壤含水率，cm3/cm3;Z為重力勢，cm。

裸土區(qū)(一維)和植被覆蓋區(qū)(二維)土壤總水勢分布特征如圖5所示(箭頭為水分運(yùn)移方向)。裸土區(qū)的平均土壤總水勢為-92.64 cm，20 cm處總水勢最大(-75.54 cm)，表現(xiàn)為土壤水自20 cm處向兩側(cè)運(yùn)移，較易接受大氣降水的補(bǔ)給;植被覆蓋區(qū)平均土壤水勢為-146.88 cm，遠(yuǎn)低于裸土區(qū)的土壤水勢，在根系密集區(qū)(20 cm，50 cm)出現(xiàn)水勢低點(diǎn)(-170 cm)，表現(xiàn)為上部、下部和側(cè)向的土壤水向植被根系處運(yùn)移，不易接受大氣降水的補(bǔ)給。

圖5 裸土區(qū)和植被覆蓋區(qū)平均土壤總水勢分布特征

3.3 植被影響下的地下水變化規(guī)律

圖6反映了裸土區(qū)和植被覆蓋區(qū)的地下水位埋深動態(tài)變化特征。在無降水期，裸土區(qū)地下水位呈現(xiàn)出緩慢的上升趨勢，表明盡管受土壤蒸發(fā)作用，深部土壤水仍維持向下運(yùn)動的趨勢;受植被蒸騰耗水影響，植被覆蓋區(qū)的地下水位平均以1.87 cm/d的速度呈現(xiàn)持續(xù)下降的趨勢。在累計(jì)124 mm降水量的情景下，裸土區(qū)的地下水位表現(xiàn)出明顯的上升趨勢，上升幅度高達(dá)27 cm;相反，植被覆蓋區(qū)的地下水位沒有回升，呈現(xiàn)出較低速度(0.29 cm/d)的地下水位下降現(xiàn)象，表明植被蒸騰耗水縮減了降水入滲補(bǔ)給量。

圖6 裸土區(qū)和植被覆蓋區(qū)地下水埋深變化

基于地下水位動態(tài)觀測數(shù)據(jù)，本次試驗(yàn)地下水蒸發(fā)(補(bǔ)給)強(qiáng)度可由式(10)估算:

Eg(Rg)=±μΔh

(10)

式中，Eg(Rg)為地下水蒸發(fā)(補(bǔ)給)強(qiáng)度，cm/d;μ為風(fēng)積沙給水度，取值為0.26[23];Δh為地下水位降(升)幅，cm。

圖7反映了裸土區(qū)和植被覆蓋區(qū)的蒸發(fā)和補(bǔ)給強(qiáng)度動態(tài)變化特征。通常，裸土區(qū)的地下水以低強(qiáng)度的入滲補(bǔ)給為主，平均補(bǔ)給強(qiáng)度為0.035 cm/d，降水時補(bǔ)給強(qiáng)度大幅增加;植被覆蓋區(qū)地下水表現(xiàn)為高強(qiáng)度的蒸發(fā)現(xiàn)象，平均蒸發(fā)強(qiáng)度為0.48 cm/d，降水時蒸發(fā)強(qiáng)度大幅減小。受植被蒸騰耗水影響，地下水面水分垂向交換量由-0.035 cm/d(入滲)增加至0.480 cm/d(排泄)，水分運(yùn)移方向和水分垂向交換量都發(fā)生明顯改變。由此可見，植被的存在一定程度改變了地下水的補(bǔ)給和排泄條件，從而引發(fā)地下水位的下降。

圖7 裸土區(qū)和植被覆蓋區(qū)地下水蒸發(fā)(補(bǔ)給)強(qiáng)度變化

3.4 植被覆蓋度增加對地下水補(bǔ)給的影響

高強(qiáng)度植被恢復(fù)的外在表現(xiàn)是植被覆蓋度的增加，在數(shù)值模型里面概化為植被葉面積指數(shù)的升高。杜春雨等[24]運(yùn)用歸一化植被指數(shù)(NDVI)推求植被葉面積指數(shù)(LAI)，并構(gòu)建了線性模型LAI=7.67NDVI-4.01，其相關(guān)系數(shù)為0.66。結(jié)合遙感數(shù)據(jù)，礦區(qū)最大的NDVI值為0.98，基于上述NDVI和LAI的經(jīng)驗(yàn)公式，推測最大LAI為3.50。筆者定義裸土條件下植被的葉面積指數(shù)(LAI)為0，模擬過程中以0.5/次的幅度逐步增加植被的葉面積指數(shù)，直至葉面積指數(shù)為3.50?？紤]到當(dāng)?shù)叵滤裆畲笥? m時，榆神礦區(qū)降水入滲系數(shù)隨埋深的增大不發(fā)生變化[25]，故本文采用3 m處的土壤水流通量估算降水入滲系數(shù)，定義為

α′=Fs/P

(11)

式中，α′為降水入滲系數(shù);Fs為平均土壤水通量，mm;P為平均降水量，mm。

圖8反映了不同葉面積指數(shù)條件下的植被耗水量和降水入滲系數(shù)變化規(guī)律。隨著葉面積指數(shù)的增加，植被耗水量呈增加趨勢，由水均衡原理可知土壤中向下運(yùn)移的水流通量逐步減小，從而導(dǎo)致降水入滲補(bǔ)給量減小，表現(xiàn)出隨著葉面積指數(shù)的增加，降水入滲系數(shù)逐步減小。數(shù)值模擬結(jié)果表明，裸土條件下(LAI=0)風(fēng)積沙的降水入滲系數(shù)為0.54，隨著葉面積指數(shù)由0增大至3.5，入滲系數(shù)減小至0.198，減小幅度高達(dá)63.3%。

圖8 植被耗水量和降水入滲系數(shù)隨葉面積指數(shù)的變化

4 討 論

近年來，隨著生態(tài)脆弱礦區(qū)植被恢復(fù)工程的實(shí)施，區(qū)域生態(tài)環(huán)境質(zhì)量得到改善，植被覆蓋度持續(xù)增加，這種現(xiàn)象在我國西部干旱缺水礦區(qū)尤為顯著(圖3)。榆神礦區(qū)地表由大面積的風(fēng)積沙組成，本文研究結(jié)果顯示，無植被覆蓋條件下風(fēng)積沙的降水入滲系數(shù)高達(dá)0.54，降水入滲直接補(bǔ)給地下水，使得該區(qū)域地下水埋深普遍較小，表現(xiàn)為地下水與植被、地表水系等關(guān)系密切[26]。在礦區(qū)植被覆蓋度大幅提升的背景下，植被蒸騰耗水量的增加導(dǎo)致土壤水分大幅度降低[27]，從而減少降水對地下水的垂向補(bǔ)給[28]。結(jié)合原位試驗(yàn)觀測和數(shù)值模擬，尹立河等認(rèn)為當(dāng)植被覆蓋率增加至100%時，降水對地下水的補(bǔ)給量基本為0[9];KEESE K E通過搭建HYDRUS數(shù)值模型研究了植被對地下水補(bǔ)給的影響，認(rèn)為當(dāng)干旱區(qū)植被覆蓋度增加后，地下水補(bǔ)給量降低幅度最高可達(dá)30倍[29]。結(jié)合本文研究成果，原位試驗(yàn)顯示植被根系區(qū)的土壤含水率出現(xiàn)較大范圍的低值區(qū)(圖4)，造成土壤總水勢的分布狀況和大小出現(xiàn)變異(圖5)，使得水分運(yùn)移方向和垂向交換量發(fā)生改變，表現(xiàn)出裸土區(qū)和植被覆蓋區(qū)截然不同的地下水位動態(tài)特征(圖6);數(shù)值模擬結(jié)果表明，隨著植被葉面積指數(shù)增加植被耗水量呈增加趨勢，降水入滲系數(shù)呈減小趨勢(圖8)。因此，榆神礦區(qū)在植被恢復(fù)的工程實(shí)踐中，應(yīng)優(yōu)先考慮植被類型和植被覆蓋度兩個基本要素，選擇耗水量較小的植被類型和適宜的植被覆蓋度，避免出現(xiàn)由于區(qū)域生態(tài)耗水量過大引發(fā)的降水入滲補(bǔ)給量大幅減小的現(xiàn)象發(fā)生。值得注意的是，筆者在研究中沒有考慮植被降水截留對地下水補(bǔ)給的影響，徐先英等研究發(fā)現(xiàn)灌木對降水的截留量可達(dá)30%以上[30]，這同樣會減少降水對地下水的補(bǔ)給量，因此本文的數(shù)值計(jì)算結(jié)果相對保守。

5 結(jié) 論

(1)礦區(qū)尺度的NDVI和NPP都呈持續(xù)增加現(xiàn)象，NDVI值從2000年的0.247增長到2015年的0.358，NPP值從2000年的12.67 g/(m2·a)增長到2015年的20.78 g/(m2·a)，表現(xiàn)為低覆蓋區(qū)連續(xù)減小，中、高覆蓋區(qū)持續(xù)增大。

(2)裸土區(qū)20 cm處土壤總水勢最大，表現(xiàn)為土壤水自20 cm處向上、下兩側(cè)運(yùn)移，較易接受大氣降水的補(bǔ)給;植被覆蓋區(qū)土壤總水勢遠(yuǎn)低于裸土區(qū)的土壤水勢，在根系密集區(qū)出現(xiàn)水勢低點(diǎn)，表現(xiàn)為上部、下部和側(cè)向的土壤水向植被根系處運(yùn)移，不易接受大氣降水的補(bǔ)給。

(3)裸土區(qū)地下水位呈現(xiàn)出緩慢上升的趨勢，以低強(qiáng)度的入滲補(bǔ)給為主;植被覆蓋區(qū)地下水位呈現(xiàn)出明顯的持續(xù)下降趨勢，以高強(qiáng)度的蒸發(fā)排泄為主。受植被蒸騰耗水影響，地下水面水分垂向交換量由-0.035 cm/d(入滲)增加至0.480 cm/d(排泄)，水分運(yùn)移方向和水分垂向交換量都發(fā)生改變。

(4)隨著葉面積指數(shù)的增加植被耗水量呈增加趨勢，土壤中向下運(yùn)移的水流通量逐步減小，導(dǎo)地下水接受降水的補(bǔ)給量減小，表現(xiàn)出隨著葉面積指數(shù)的增加，降水入滲系數(shù)逐步減小。數(shù)值模擬結(jié)果表明，裸土條件下風(fēng)積沙的入滲系數(shù)為0.54，隨著葉面積指數(shù)由0增大至3.5，入滲系數(shù)減小至0.198，減小幅度高達(dá)63.3%。

(5)榆神礦區(qū)在植被恢復(fù)的工程實(shí)踐中，應(yīng)優(yōu)先考慮植被類型和植被覆蓋度兩個基本要素，選擇耗水量較小的植被類型和適宜的植被覆蓋度，避免出現(xiàn)由于區(qū)域生態(tài)耗水量過大引發(fā)的降水入滲補(bǔ)給量大幅減小的現(xiàn)象發(fā)生，以實(shí)現(xiàn)生態(tài)脆弱礦區(qū)植被恢復(fù)與地下水資源保護(hù)的協(xié)調(diào)發(fā)展。