高寒內陸河流域植被覆蓋增加對地下水補給的影響

鄭 麗,金 鑫,3,*,金彥香,3,傅 笛,翟婧雅

1 青海師范大學地理科學學院,西寧 810016 2 青海省自然地理與環境過程重點實驗室,西寧 810016 3 高原科學與可持續發展研究院,西寧 810016

地下水在干旱區內陸河流域水循環中占有主導地位,是自然系統中最重要的基礎性資源之一[1—3]。植被是土壤和大氣間水量交換的關鍵驅動因素,其通過降雨截留和消耗土壤水分等方式間接影響地下水的補給,在地下水循環中具有重要作用[4—5]。相關研究表明,在冰島,當草地及荒地被人工林代替后,植被覆蓋增加,區域地下水補給量降低了10%[6]。在荷蘭,當農田被林地代替后區域地下水補給量減少約485 mm/a[7]。澳大利亞昆士蘭州東南部部分流域的稀疏草地被稀疏林地替代后,地下水補給量將減少25%；當稀疏草地被茂密林地替代后,地下水補給量將減少 48%[8]。我國黃土高原涇河流域的研究表明,平均森林覆蓋率每增加10%,將導致流域年地下水補給量低11.1 mm[9]。這些區域尺度的研究已經證實,植被覆蓋增加對地下水補給的影響非常顯著。

一般來說,干旱/半干旱區的地下水補給比濕潤地區更容易受到地表覆蓋條件的影響[10—11]。由于降水少且季節性強、補給量通常很小等原因,干旱/半干旱區地下水補給量的估算具有挑戰性,一般的土壤水資源平衡法[3]、水文資料分析法[8]等難以適用。目前,水文模型因其可對自然界中復雜水循環過程近似描述和再現,可作為區域植被變化之水文效應研究的有效手段[12—14]。如SWAT(Soil and Water Assessment Tools)模型,在我國干旱/半干旱區內陸河流域的水循環及其影響因素研究中具有廣泛應用[15—17]。區域植被覆蓋的增加或減少可能由人類活動導致的土地利用類型變化引起,也可能因氣候變化導致的植被生長狀態改變等因素引起。SWAT模型中的土地利用更新模塊可在一定程度上表達土地利用類型變化情況[10],但其缺乏對于其它因素導致的植被覆蓋度時空變異性的有效考慮：在SWAT模型中,葉面積指數(LAI)作為表征地表植被生長狀況、植被覆蓋度等的重要指標,是連接植被與水循環的重要橋梁[18]。但是,SWAT僅基于溫度計算LAI,忽略了降水、地形等因素。這在濕潤的熱帶、亞熱帶地區水文過程模擬中被更多被關注[19]。但在降水時空異質性強、地形復雜的干旱區內陸河流域,這一問題被忽視。

基于衛星遙感的長時間序列LAI數據產品因能較好展示LAI的空間異質性、季節性變化等,在評估區域植被生長狀況、植被覆蓋條件等方面具有重要應用。如,中分辨率成像光譜儀葉面積指數MODIS(MODerate-resolution Imaging Spectroradiometer) LAI[20]、全球地表衛星葉面積指數GLASS (Global Land Surface Satellite) LAI[21]、全球長時間序列葉面積指數GLOBMAP (Long-term Global Mapping) LAI[22]等。這些產品在不同區域的時空精度具有差異,且部分產品在中國西北部分區域存在空值區[21—23]。相對來說,GLASS LAI 產品在中國范圍內具有相對較高的準確性[21,23],并且具有較好的時空連續性和完整性[21]。基于該數據集進行區域生態-水文過程模擬更為可靠。

SWAT模型僅采用簡單的線性方程來表達地下徑流變化與潛水面埋深變化速率之間的關系,無法真實表達流域地下水過程[24]。利用成熟的地下水模型MODFLOW取代SWAT模型的地下水模塊,二者通過同參變量的傳輸和反饋進行耦合,既發揮SWAT模型在作物生長、地表水模擬中的優勢,又結合了MODFLOW對于地下水模擬的長處。土壤和水評價工具與模塊化有限拆分地下水流耦合模型(SWAT-MODFLOW)耦合模型已在美國科羅拉多州南普拉特河流域[25]、非洲林波波河流域[26]、巴音河流域[10,27]的地下水-地表水轉換關系、地下水補給、地下水位預測等方面得到了應用且效果優于原始SWAT模型。

巴音河位于柴達木盆地東北緣,是盆地第四大內陸河。近年來,在人工植被恢復以及氣候暖濕化背景下,巴音河流域植被覆蓋條件明顯變好[10]。本文為探究干旱區內陸河流域植被覆蓋變化對地下水補給的影響,以巴音河中下游為例,改進SWAT模型,采用GLASS LAI數據代替其LAI計算模塊,再結合流域土地利用/覆被類型變化,準確刻畫區域植被覆蓋度變化。將改進后的SWAT模型與MODFLOW模型耦合。在準確模擬植被覆蓋變化下流域地下水循環過程的基礎上,基于控制變量法分析巴音河中下游地區植被覆蓋變化對地下水補給的影響。本研究可為干旱區內陸河流域生態環境保護與可持續發展提供科學依據。

1 材料與方法

1.1 研究區概況

巴音河位于青海省德令哈市,是柴達木盆地東北部最大的內陸河,也是德令哈市境內最大的河流。該河流中下游段為德令哈市區北部黑石山水庫以下到可魯克湖、尕海湖入湖口以上的流域[27],該流域中游段由北向南穿過德令哈市市區,下游段由東向經郭里木、戈壁鄉最終匯入克魯克湖-托索湖湖區(圖1)。巴音河出山后河水大量下滲,河道徑流量較少,部分河段河水全部潛入地下。氣候寒冷、干旱,水資源短缺,日照時間長,屬于高原荒漠半荒漠氣候,主要的土地覆蓋類型為荒漠和草地[28]。

圖1 研究區概況圖Fig.1 The study area

1.2 GLASS LAI數據集

本文采用時空分辨率為 8 d/250 m 的最新 GLASS LAI v5 數據,時間跨度為 2001—2019年。該數據集由北京師范大學全球變化處理與分析中心(http://www.glass.umd.edu/Download.html)發布。該數據集采用了加權線性組合的方式,首先將 MODIS LAI和基于生物地球物理參數的葉面積指數(CYCLOPES LAI) 融合,產生二者的融合 LAI[20],再與經過預處理的中分辨率成像光譜儀/高級甚高分辨率輻射計(MODIS/AVHRR) 地表反射率數據,建立廣義回歸神經網絡訓練樣本,訓練出地表反射率與 LAI 值的關系模型,最后將 MODIS/AVHRR 地表反射率數據作為關系模型輸入數據,得到全球長時間序列的GLASS LAI 產品[29]。使用該數據之前進行了特異值去除并將有效值歸為 0—10。

1.3 土地利用/覆被數據的獲取

本研究基于攜帶專題制圖儀的美國陸地衛星系列第五顆衛生Landsat5 TM、攜帶陸地成像儀的美國陸地衛星計劃第8顆衛星Landsat8 OLI、全球環境與安全監測計劃的第二顆衛星(Sentinel-2A) 影像數據獲取巴音河中下游2001—2019年土地利用/覆被數據。首先在研究區針對耕地、林地、草地、水體、城鎮用地、裸地六種地類開展野外調查。利用 Landsat5 TM(2011 年 11 月之前)、Landsat8 OLI(2013 年 4 月之后)、以及 Sentinel-2A(2015 年 8 月之后)影像數據,創建歸一化植被指數(NDVI)、增強型植被指數(EVI)和歸一化水體指數(NDWI)時間序列數據集,結合野外調查結果,實現各地類光譜特征、植被指數特征、時空變化特征等特征提取,進而利用隨機森林算法對研究區五種地類進行分類。運用混淆矩陣對獲得的土地利用/覆被數據進行精度評價,其總體精度達到94%,kappa系數為0.96。圖2顯示了本研究獲取的2001年及2019年研究區土地利用/覆被數據。

圖2 2001年及2019年巴音河中下游土地利用/覆被數據Fig.2 Land use/ land cover maps of the Middle and Lower Reaches of the Bayin River in year 2001 and 2019

1.4 SWAT模型改進

本研究利用GLASS LAI代替原始SWAT模型基于理想葉面積發育模型計算的LAI。過程如下：將基于GLASS數據集的250 m 分辨率、8 d LAI 數據集,映射至SWAT模型每個基本計算單元(HRU) 上。由于 LAI數據集與HRU的空間尺度不同,首先對二者進行空間疊加。疊加完成后,對于 1 個 HRU 包含多個網格的情況,對所有網格的 LAI 值進行平均后,賦予相應的 HRU。由于SWAT 模型模擬水文過程的時間步長為 1 d,故在 LAI 數據集與動態 HRU 完成對應后,對每個 HRU上的8 d LAI 數據進行插值,使其時間分辨率變為 1 d,從而實現GLASS LAI與SWAT 模型的耦合(如圖3)。生物量及流域產流、產沙等后續過程均基于改進后的作物生長模塊計算。由于改進后的SWAT模型更能有效刻畫區域植被動態變化,故本研究將其命名為DVSWAT。

圖3 GLASS LAI與SWAT -MODFLOW模型的耦合Fig.3 Coupling of the GLASS LAI and SWAT -MODFLOWHRU：水文響應單元 Hydrologic response unit；GALSS：全球地表衛星 Global land surface satellite；LAI：葉面積指數 Leaf area index；SWAT：土壤和水評價工具 Soil and water assessment tools；MODFLOW：模塊化有限差分地下水流耦合模型 Modular finite difference groundwater flow model

1.5 DVSWAT -MODFLOW模型及構建

通過GIS平臺建立SWAT模型的基本計算單元與MODFLOW模型的基本計算單元之間的映射關系,然后將每個HRUs分配到MODFLOW模型的相應網格單元上,將DVSWAT模型計算出的地下水補給量、蒸發量等數值用來作為地下水模擬的邊界條[27]。同時,利用映射關系把MODFLOW模型所計算出的地下水排泄量添加到DVSWAT模型的相應地表水計算單元中,以此來實現DVSWAT模型與MODFLOW模型之間的耦合(圖3)[30]。其中,由DVSWAT模型負責計算降水對地下水的直接補給以及灌溉對地下水的補給。由DVSWAT模型負責計算河水水位并結合MODFLOW模型的河流子模塊計算河水滲漏對地下水的補給。

本文選取了德令哈站的逐日降水、溫度、相對濕度、風速及太陽輻射數據作為DVSWAT模型輸入數據,結合30 m分辨率的DEM數據(http://gdex.cr.usgs.gov/gdex/)、全國1∶400萬土壤類型數據、2001、2010、2015、2019年土地利用類型數據以及2001—2019年逐日GLASS LAI數據建立DVSWAT模型。土壤屬性數據查閱自《青海土壤》[31]及《德令哈市志》[32]。

采用研究區鉆孔數據、地下水位觀測數據、河網數據、數字高程模型(DEM)數據等建立MODFLOW模型,建模方法主要參考已有研究[27,33]。以巴音河中下游流域邊界以及淺層含水層底板作為隔水邊界,流域入口及出口作為定水頭邊界,河網作為河流邊界。將研究區概化為250 m×250 m的規則網格,全區共251行、272列。本研究根據需要僅模擬淺層含水層地下水循環過程,且將淺層含水層概化為單層非均質各向同性,將地下水運動概化為二維非穩定流[27,33]。其中,給水度、滲透系數、河流滲透系數等參數主要參考已有研究[33],再結合研究區地下水位觀測數據進行微調。

采用了基于遙感的全球陸地蒸散發阿姆斯特丹模型GLEAMv3(最新)蒸散發數據(https://www.gleam.eu/)、地下水位觀測數據對DVSWAT-MODFLOW模型進行校準。其中,GLEAM數據集空間分辨率為0.25°,時間分辨率為1天。該數據集在中國區域內精度較高[34],被廣泛應用于水文模型校準及驗證[35—36],具有較高可靠性。此外,利用決定系數(R2)、納什系數(NSE)、誤差百分數(PBIAS)、絕對誤差(AE)評價模型的適用性[10,37]。

2 結果與討論

2.1 巴音河中下游土地覆被類型變化

基于2001年和2019年土地利用/覆被數據(圖2)提取了研究區土地利用轉移矩陣(表1)。從2001年至2019 年,土地覆被類型發生了明顯變化。其中,草地、建設用地、林地和水體面積增加,2019年建設用地、林地和水體地面積約是2001年的5.41倍、2倍和2.35倍,草地面積增加了98.96%；裸地和耕地面積減少,分別減少了30.75%和52.88%。本研究主要關注植被覆蓋增加對流域地下水補給的影響。故將2001及2019年土地利用/覆被數據進行空間疊加,提取研究區植被恢復斑塊,結果如圖4所示。其中,草地恢復主要在流域北部及東南部山區,林地恢復主要在流域中部農業區。

表1 2001—2019 年土地利用轉移矩陣/km2Table 1 Land use transfer matrix from year 2001 to 2019

圖4 研究區植被恢復情況Fig.4 Revegetation of the study area

2.2 巴音河中下游植被覆蓋度時空變化

圖5顯示了基于GLASS LAI數據集的2001年及2019年研究區7月平均LAI空間分布情況。從空間上看,研究區北部高海拔山區、中部農業灌區、中南部地下水位高值區植被覆蓋狀況相對較好,這些區域2019年7月植被覆蓋度明顯高于2001年同期。2001年研究區7月(植物生長旺季)平均LAI值為4.37,2019年7月平均LAI值為5.63,比2001年增加了28.83%。這與巴音河流域人工植被恢復及氣候暖濕化有關[27,33]。

圖5 2001 年及2019 年巴音河中下游年平均LAI空間分布Fig.5 Distribution of the annual average LAI in the Middle and Lower Reaches of the Bayin River in year 2001 and 2019

圖6 基于植被動態變化的土壤和水評價工具與模塊化有限拆分地下水流耦合模型(DVSWAT-MODFLOW)模型月蒸散發模擬效果Fig.6 The monthly ET simulated by the DVSWAT-MODFLOWNSE：納什系數 Nash-Sutcliffe efficiency；PBIAS：誤差百分數 Percent bias；R2：擬合優度 Absolute error；ET：蒸散發 Evapotranspiration

2.3 DVSWAT-MODFLOW蒸散發模擬效果評價

巴音河流至中下游后大量滲漏,補給地下水。故河道徑流量較小且無觀測數據。基于河道徑流數據的傳統參數校準方法難以適用。本研究基于GLEAM v3蒸散發數據集校準DVSWAT-MODFLOW模型蒸散發相關參數。GLEAM v3蒸散發數據集在資料缺乏區水文模型校準中具有重要應用[33]。圖6顯示了DVSWAT-MODFLOW模型的月蒸散發模擬效果。各子流域R2值在0.83以上,NSE值在0.68以上,PBIAS絕對值在22%以內。可見模型對于月蒸散發的模擬效果較好[38]。圖6亦顯示了DVSWAT-MODFLOW模型模擬的各子流域多年平均蒸散發分布情況。流域下游子流域地下水位較高,年蒸散發量相對較大,超過200 mm。個別子流域年蒸發量超過300 mm,這是相應子流域僅包含耕地且被大水漫灌所致。

2.4 DVSWAT-MODFLOW地下水位模擬效果評價

圖7 巴音河中下游月地下水位實測及模擬值對比Fig.7 The simulated and observed groundwater level in the Middle and Lower Reaches of the Bayin River1—10號觀測井觀測數據涵蓋2019年各月；11號觀測井觀測數據涵蓋2009—2011年各月；12—14號觀測井觀測數據涵蓋2013—2015年各月；15號觀測井觀測數據涵蓋2004—2005年各月；16—20號觀測井觀測數據涵蓋2001—2014年各月

圖7顯示了DVSWAT-MODFLOW模型月地下水位模擬效果。20個地下水位觀測井的空間分布較為均勻(圖1)。流域南部觀測井較為密集,月地下水位模擬效果相對優于其它區域,R2達到0.95以上,絕對誤差(AE)小于0.35 m。1號及17號觀測井處絕對誤差較大,接近1 m。總體上來說,各觀測井月地下水位模擬效果較好。

圖8 不同地表覆蓋類型對應的年平均地下水補給量 Fig.8 The annual average groundwater recharge corresponding to the different land cover types

2.5 植被覆蓋增加對地下水補給的影響

圖8顯示了研究區不同地表植被覆蓋類型對應的年平均地下水補給量。為剔除氣象因素、坡度、土壤類型等對地下水補給的影響,本研究對具有相同氣象條件、土壤類型及坡度的HRUs所對應的土地覆被類型進行了統計。在年平均降水量為244.23 mm的情況下,耕地對應的年平均蒸散發量、地下水補給量最大,分別為393.12 mm及373.46 mm這主要是由大水漫灌所致。裸地對應的年平均蒸發量為178.14 mm,年平均地下水補給量為55.68 mm。林地對應的年平均蒸發量為193.04 mm,年平均地下水補給量為39.53 mm。草地對應的年平均蒸發量為179.79 mm,年平均地下水補給量為50.14 mm。可見研究區地下水補給量耕地>裸地>草地>林地。表明巴音河中下游耕地、裸地的大量減少和草地的大面積增加使地下水補給量降低。

為探究巴音河中下游植被覆蓋度增加對地下水補給量的影響,本研究選取年平均LAI增加較明顯的26號子流域為典型區(圖9),將2001年LAI與2019年LAI分別輸入DVSWAT-MODFLOW模型。結果表明,2001年該典型區平均LAI在0—0.48之間,在年平均降水量為244.23 mm的情況下,對應的年平均蒸散發量為195.53 mm,對應的年平均地下水補給量為38.14 mm。2019年該典型區平均LAI在0—0.79之間,在年平均降水量為244.23 mm的情況下,對應的年平均蒸散發量為208.10 mm,對應的年平均地下水補給量為26.92 mm(圖9)。可見植被覆蓋度增加使蒸散發量增加,地下水補給量減少。

圖9 不同植被覆蓋度對應的地下水補給量Fig.9 The annual average groundwater recharge corresponding to the different vegetation coverage

為在全流域尺度綜合探究植被覆蓋增加對流域地下水補給的影響,首先基于研究區2001年對應的土地利用/覆被數據以及LAI數據建立DVSWAT-MODFLOW1。在保證它地表特征參數及氣象數據等輸入參數不變的前提下,基于圖4所顯示的植被恢復斑塊建立研究區草地、林地面積增加(其它土地利用類型不變)后的土地利用/覆被數據。再結合2019年LAI數據建立DVSWAT-MODFLOW2。圖10顯示了DVSWAT-MODFLOW1及DVSWAT-MODFLOW2模擬的年際尺度河水滲漏補給地下水量。表明植被覆蓋增加將使2001—2019年河水滲漏補給量減少。減少量最小值為2 mm,出現在2006年,對應的年降水量亦最小,為143 mm。最大減少值為10.34 mm,出現在2019年,對應的年降水量為256.5 mm。圖10顯示了模型模擬的地下水其它補給量(包括降水直接補給、灌溉直接補給等),可見植被覆蓋增加后,地下水其它補給量減少4.1—16.18 mm。最小減少量出現在2006年。最大減少量出現在2002年,對應的年降水量為302 mm。可見年降水量可在一定程度上決定植被覆蓋增加對地下水補給量的影響強弱。在月際尺度上,植被覆蓋增加將使河水滲漏補給量減少0—7.18 mm(圖11)。減少量為0的月份一般處于較為寒冷、干旱的植物休眠期或河流枯水期。河水滲漏補給最大減少量出現在2012年7月,對應的月降水量為109 mm。此外,植被覆蓋增加將使地下水其它補給量減少0—7.85 mm(圖11)。減少量為0的月份同樣處于較為寒冷、干旱的植物休眠期或其它降水較少時期。最大減少量出現在2016年8月,對應的月降水量為113 mm。植被覆蓋增加對地下水補給量的影響較強的月份集中于植物生長旺盛且降水較多時期。綜合來看,2019年植被覆蓋情況對應的年際及月際尺度地下水補給量較2001年分別減少了6.1—26.52 mm以及0—15.03 mm。

圖10 植被覆蓋增加對年地下水補給的影響Fig.10 The impacts of the increasing vegetation coverage on annual groundwater rechargeDVSWAT：耦合植被覆蓋變化的SWAT模型 Soil and water assessment tool coupled with dynamic vegetation coverage

圖11 植被覆蓋增加對月地下水補給的影響Fig.11 The impacts of the increasing vegetation coverage on monthly groundwater recharge

3 結論

(1)DVSWAT-MODFLOW模型的月蒸散發模擬效果較好,各子流域R2值在0.83以上,NSE值在0.68以上,PBIAS絕對值在22%以內。月地下水位模擬效果亦較好,各觀測井R2值在0.92以上,絕對誤差在0.28—0.97 m之間。

(2)2019年巴音河中下游林地及草地面積較2001年增加明顯,分別增加了5.41倍及98.96%。耕地及裸地面積減少明顯。此外,2019年植物生長旺季平均LAI值為5.63,比2001年增加了28.83%。

(3)年際尺度上,巴音河中下游植被覆蓋增加導致河水滲漏補給地下水量減少2—10.34 mm,導致地下水其它補給量減少4.1—16.18 mm。植被覆蓋增加對年際尺度地下水補給量的影響強弱在一定程度上取決于年降水量。

(4)月際尺度上,巴音河中下游植被覆蓋增加導致河水滲漏補給地下水量減少0—7.18 mm,地下水其它補給量減少0—7.85 mm。植被覆蓋增加對月際尺度地下水補給量影響較強的月份集中于植物生長旺盛且降水較多時期。