氣候和土地利用變化對渭河流域水資源短缺的影響

岳思羽, 李懷恩, 趙 麗

(1.陜西理工大學 化學與環境科學學院, 陜西 漢中 723001; 2.省部共建西北旱區生態水利國家重點實驗室 西安理工大學, 西安 710048; 3.漢中市生態環境保護綜合執法支隊, 陜西 漢中 723001)

渭河是黃河最大的支流，發源于甘肅省渭源縣鳥鼠山，流經隴東高原、天水盆地、關中平原，至潼關匯入黃河，流域面積13.48萬km2，涇河和北洛河分別是渭河的第一、第二大支流[1]。渭河位于我國干旱半干旱區，水資源總量匱乏[2]；但流經區域人口密集、經濟發達，為滿足生產生活用水需求，渭河沿岸引水量不斷增加，導致水資源供需矛盾突出[3]。根據陜西省、甘肅省和寧夏回族自治區《水資源公報》，2018年渭河水資源量為95億m3，而渭河流域的用水總量高達65億m3，其中地表水用量為36億m3，約占用水總量的55%；農業用水量為31億m3，是主要的用水行業。近年來，隨著氣候變化和人類活動的影響，渭河天然來水量減少，用水量卻呈增加趨勢，河道生態基流難以保障[4]，水資源短缺已經成為制約渭河流域水生態環境保護和經濟社會發展的重要因素。

為評價水資源短缺狀況，Falkenmark指數[5](用水量與可用水量的比值)得到廣泛應用。在此基礎上，Ohlsson[6]將聯合國提出的人類發展指數(Human Development Index)作為Falkenmark指數的加權，得到社會缺水指數(Social Water Stress Index)。McNulty等[7]以需水量和供水量的比值作為供水壓力指數(Water Supply Stress Index)，這些指標主要考慮人類用水需求的滿足程度。隨著河道生態基流研究的深入[8]，河流生態系統對水資源的需求逐漸得到重視，Asheesh[9]建立了由人口增長率、可用水量、生活用水量、工業用水量和生態需水量等組成的水資源短缺指數(Water Scarcity Index)，Smakhtin等[10]將河流的生態環境需水作為缺水指數(Water Stress Indicator)的重要組成部分。與傳統水資源短缺指標相比，這兩項指標同時考慮了人類和生態系統的用水需求，評價結果更為客觀。

針對渭河流域水資源短缺和河道生態基流保障難度大等問題，本文引入考慮河道生態基流需求的流域水資源短缺評價方法，利用生態系統服務與權衡綜合評估模型(Integrated Valuation of Ecosystem Services and Tradeoffs，InVEST模型)產水量模塊對渭河流域進行子流域尺度的水資源短缺評價，探索水資源短缺的時空變化及其影響因素，并揭示水資源短缺的關鍵區域，為渭河流域的可持續管理以及河道生態基流保護提供一個更為全面的視角。

1 試驗數據與方法

1.1 研究方法

1.1.1 產水量模型 InVEST模型產水量模塊[11]根據Budyko曲線建立，以每個單元的降水量及蒸發量估算每個單元的產水量，再進行求和平均得到子流域的產水量。模型不區分地表水、地下水和基流，只假設所有的水量都從各個單元通過這3種途徑之一到達模型的控制點。模型的輸出結果為子流域尺度。

柵格單元x的年產水量WY(x)為：

(1)

式中:AET(x)為柵格單元x的年實際蒸發量(mm);P(x)為柵格單元x的年降水量(mm)。

(2)

式中:PET(x)為柵格單元x的潛在蒸發量(mm);ω(x)為描述柵格單元x自然氣候—土壤特性的非物理參數。

PET(x)的計算公式為：

PET(x)=Kc(lx)ET0(x)

(3)

式中:ET0(x)為柵格單元x的參考作物騰發量,指一種假想的參考作物(如苜蓿)冠層騰發速率,用于表征氣象條件對騰發過程的影響(mm);Kc(lx)為土地利用類型lx的植被蒸散系數,由土地覆蓋的植被類型決定。

ω(x)是經驗參數，計算公式為[14]：

(4)

式中:AWC(x)為柵格單元x的土壤有效含水量(mm);指土壤中可供植物吸收的水量;Z為表示當地降水模式和水文地質特征的經驗常數,取值為1～30。

AWC(x)的計算公式為：

AWC(x)=min(Rest.layer.depth, root.depth).PAWC

(5)

式中:Rest.layer.depth為土壤根系最大埋深,本文用土層深度代替(mm);root.depth為根系深度(mm);PAWC為植被可利用含水量(一般是田間持水量和萎蔫濕度的差值含水量),取值為0～1[11]。

PAWC的計算公式為[15]：

PAWC=54.509-0.132sand-0.003(sand)2-0.055silt-0.006(silt)2-0.738clay+0.007(clay)2-

2.699OM+0.501(OM)2

(6)

式中:sand為土壤砂粒含量(%);silt為土壤粉粒含量(%);clay為土壤黏粒含量(%);OM為土壤有機質含量(%)。

(2) 對于其他土地利用類型(水域、濕地等)，實際蒸散量可直接根據參考作物騰發量計算得出，且不超過降水量：

AET(x)=min[Kc(lc)·ET0(x),P(x)]

(7)

1.1.2 水資源短缺評價方法 InVEST模型產水量模塊可以通過輸入耗水量數據，對流域的水資源短缺進行評價。產水量模塊中所定義的“耗水量”指為滿足人類需求從天然水生態系統中調取并耗竭的水資源量[11]。這部分水在使用過程中被耗竭，不再回到原來的水生態系統，不同于為人類提供生物多樣性保護、休閑娛樂等水生態系統服務功能的水資源。在產水量模塊中，耗水量數據是以土地利用類型為類別輸入的，指各土地利用類型中從流域水平衡中去除的水資源量，一般指農業、工業和生活的耗水量。

Boithias等[16]提出了在InVEST模型產水量模塊下考慮河道生態基流需求的流域水資源短缺評價方法：

(8)

式中:WSi為子流域i的水資源盈缺指數;WYi為子流域i的產水量(m3);EFi為子流域i的河道生態基流量(m3);CWi為子流域i的耗水量(m3)。

該方法在充分考慮河道生態基流需求的基礎上，利用產水量模塊對流域的水資源短缺狀況進行評價。根據現有水資源短缺指標的研究成果[17]，WSi<1代表水資源短缺的子流域，即該子流域內的產水量不能同時滿足河道生態基流和耗水量的需求，河道生態基流可能被占用，而且WSi值越小，說明水資源短缺狀況越嚴重；WSi≥1代表水資源盈余的子流域，即該子流域內的產水量在滿足河道生態基流和耗水量的需求后仍有盈余；若CWi為0，認為WSi≥1，屬于水資源盈余的子流域。

1.2 數據來源與處理

1.2.1 情景設置 根據渭河流域26個氣象站1980—2017年的降水量數據，將歷年各氣象站的降水量數據進行算術平均，得到1980—2017年渭河流域的歷年降水量值，并進行頻率分析。選擇1992年(p=25%)，1989年(p=50%)和2004年(p=75%)作為典型年，并設置6種情景(表1)，研究不同氣候情景和不同土地利用/覆被情景下渭河流域的水資源短缺變化。

表1 不同氣候和不同土地利用/覆被情景設置

1.2.2 數據來源

(1) 模型輸入數據。InVEST模型產水量模塊需要輸入的圖形數據有降水量、參考作物騰發量(ET0)、土地利用/覆被、土層深度、PAWC和子流域定義圖。根據1992年、1989年和2004年各氣象站的年降水量數據，用克里金插值法得到渭河流域不同典型年的降水量數據。采用聯合國糧食及農業組織(FAO)推薦并修訂的彭曼公式[18]計算各氣象站在1992年、1989年和2004年的年均ET0，并用克里金插值法得到渭河流域不同典型年的ET0數據，各氣象站的所有數據均來自中國氣象數據網。以1980年和2015年的Landsat TM/ETM遙感影像為主要數據源，通過人工目視解譯得到渭河流域土地利用數據，在ArcGIS 10.4下將土地利用類型分為6個一級類型，即耕地、林地、草地、水域、建設用地和未利用地，數據來自中國科學院資源環境科學數據中心。渭河流域的PAWC數據用式(6)計算，為0～0.17，所需數據均來自寒區旱區科學數據中心。渭河流域土層深度數據由寒區旱區科學數據中心提供，為300～1 000 mm。以渭河流域的數字高程圖為基礎，設定集水面積閾值為400 km2，并在ArcGIS 10.4下生成211個子流域。

模型輸入的參數植被蒸散系數Kc為0～1.5，用來修正ET0得到潛在蒸發量，來自FAO的推薦值。根系深度指土地利用類型上95%的植被根系能達到的最大深度，本文計算采用Allen等[19]推薦的標準作物根深系數。Z常數表征區域降水分布和水文特征，取值設定為1～30[11]，吳娜等[20]用InVEST模型產水量模塊對渭河甘肅段進行研究時得到Z常數取值為6，本文將Z常數的初始取值范圍設定為5～7。

(2) 河道生態基流量。華縣斷面是渭河流域的控制斷面，現有河道生態基流基準值為20 m3/s[21]，計算可得河道生態基流的對應水量為6.31億m3。假設河道生態基流均勻分布在華縣斷面以上所有子流域的每個柵格單元x內，根據對應水量和華縣斷面以上的柵格單元數量得到每個柵格單元的年河道生態基流量EF(x)為6 000 m3，柵格精度為1 000 m×1 000 m。為簡化計算，本文認為渭河流域除華縣斷面以上子流域外的EF(x)也都取值為6 000 m3，即渭河流域的EF(x)均為6 000 m3。EFi為子流域i內所有EF(x)之和。

(3) 耗水量。本文采用一種相對簡單的方法估算渭河流域不同土地利用類型的耗水量CW。(1) 根據《陜西省水資源公報》、《甘肅省水資源公報》和《寧夏回族自治區水資源公報》，得到不同典型年渭河流域的農業耗水量、工業耗水量和生活耗水量。(2) 根據產水量模塊中對耗水量的定義，將農業耗水量對應土地利用類型中耕地的耗水量，工業耗水量和生活耗水量對應土地利用類型中建設用地的耗水量，其余土地利用類型的耗水量為0 m3。(3) 假設在同一土地利用類型中耗水量在各個柵格單元內均勻分布，且不隨時間變化，這樣可根據年耗水量和柵格單元數量得到每個柵格單元的耗水量CW(x)：耕地CW(x)在豐水年為4萬m3，平水年為5萬m3，枯水年為6萬m3；建設用地CW(x)在不同典型年均為15萬m3；其余土地利用類型CW(x)在不同典型年均為0 m3；柵格精度為1 000 m×1 000 m。(4) CWi為子流域i內所有CW(x)之和。

1.2.3 模型率定與驗證 選取NSE[22]、相對誤差RE和決定系數R2對模型的適用性進行評價，1980—2000年作為率定期，2001—2017年作為驗證期。率定期指標計算結果：NSE=0.79，RE=2.15%，R2=0.86；驗證期指標計算結果：NSE=0.66，RE=-3.91%，R2=0.85，說明InVEST模型產水量模塊可以用于渭河流域的模擬研究，模型最終輸入的Z常數為6.2。

1.3 氣候和土地利用變化的影響分析

根據水資源盈缺指數WSi的計算結果，利用GeoDa模型，分別對降水量、土地利用強度與水資源盈缺指數WSi進行雙變量空間相關性分析，研究氣候和土地利用因素對水資源短缺時空變化的驅動特征?？臻g權重矩陣設置為簡單的二進制鄰接矩陣。

(1) 全局空間自相關計算方法。全局空間自相關用莫蘭指數I計算，用于研究空間單元的總體關聯性，并判斷是否存在空間集聚現象，計算公式為[23]：

(9)

I的取值范圍為-1～1，I>0表示某空間屬性值與周邊屬性值具有空間集聚效應，且I越接近1，空間集聚效應越顯著；I<0表示某空間屬性值與周邊屬性值具有空間差異效應，且I越接近-1，空間差異效應越顯著；I=0表示不存在空間自相關。

(2) 局部空間自相關計算方法。局部空間自相關用LISA指數計算，用于判斷每一個空間單元的屬性與鄰近單元屬性的相關程度，計算公式為[24]：

(10)

Ii>0表示空間單元相關性強，包括高—高型和低—低型；Ii<0表示空間單元差異大，包括高—低型和低—高型。高—高型指降水量或土地利用強度與WSi都高；低—低型指降水量或土地利用強度與WSi都低；高—低型指降水量或土地利用強度高，而WSi低；低—高型指降水量或土地利用強度低，而WSi高。

2 結果與分析

2.1 產水量的時空變化

將6種情景的基礎數據分別輸入已經率定過的InVEST模型產水量模塊，得到渭河流域單位面積產水量的時空分布及其變化特征，輸出結果為子流域尺度(圖1)。

圖1 渭河流域單位面積產水量的時空變化

渭河流域的單位面積產水量由北向南逐漸遞增。涇河流域和北洛河流域的單位面積產水量多低于1 000 m3/hm2，而渭河干流的單位面積產水量可以達到1 500 m3/hm2。對典型年的產水量進行比較，豐水年(情景1和情景2)渭河流域的單位面積產水量平均值為933 m3/hm2，高于平水年(情景3和情景4)約589 m3/hm2和枯水年(情景5和情景6)約289 m3/hm2，在豐水年秦嶺山區子流域的單位面積產水量甚至超過2 000 m3/hm2。渭河流域降水量南多北少，ET0呈現北高南低的趨勢，導致產水量空間分布不均勻；而且，不同典型年渭河流域的降水量和ET0差異明顯，使得產水量的時空變化特征顯著。

2.2 水資源短缺的時空變化

由渭河流域的產水量、耗水量及河道生態基流量，用式(8)計算得到不同情景下渭河流域水資源短缺的時空分布及變化特征，輸出結果為子流域尺度，見圖2。

圖2 渭河流域水資源短缺的時空變化

渭河流域的水資源短缺區域集中在2個地區，但原因各不相同。(1) 渭河干流上游和涇河流域北部，該區域降水量少、蒸發量高，由產水量圖(圖1)可知，該區域的單位面積產水量一般低于500 m3/hm2，產水量少是水資源短缺的原因。(2) 渭河干流中下游，該區域是沖積平原，地勢平坦、土壤肥沃，耕地數量多、面積大，而且聚集著關中平原城市群，是陜西省乃至我國西北地區人口最密集、最發達的地區之一，由耗水量可知，耕地和建設用地的耗水量為40～150 mm，是水資源短缺的原因。

分析氣候因素對水資源短缺的影響：豐水年(情景1和情景2)，除涇河流域北部的少量子流域以外，大部分區域都可以同時滿足河道生態基流和人類耗水量的需求；到平水年(情景3和情景4)，水資源短缺的區域不斷擴大，渭河干流中下游出現了多個水資源短缺的子流域；枯水年(情景5和情景6)渭河流域的水資源短缺程度明顯比豐水年和平水年嚴重，僅在以林地和草地為主的子流域內沒有出現水資源短缺的情況。

分析土地利用因素對水資源短缺的影響：受人類活動干擾較大區域的水資源短缺程度比干擾小的區域加劇，如渭河干流中下游多是以耕地和建設用地為主的子流域，是水資源短缺區域，而秦嶺山區多是以林地和草地為主的子流域，在不同典型年都沒有出現水資源短缺的狀況。2015年土地利用情景下(情景2，情景4和情景6)渭河流域的水資源短缺狀況整體好于1980年土地利用情景(情景1，情景3和情景5)，原因有兩個：(1) 2015年土地利用情景下渭河流域的產水量整體高于1980年情景。(2) 受退耕還林(草)政策的影響，1980—2015年渭河流域的耕地面積下降了1 909 km2，這樣使部分以耕地為主的子流域耗水量下降，因此緩解了渭河流域的水資源短缺狀況。

實際上，根據水平衡原理，水資源短缺是氣候和土地利用綜合作用的結果。(1) 氣候：降水和蒸發是氣候影響水資源短缺程度的主要形式，枯水年降水量減少，蒸發量增高，導致產水量急劇減少，這就是枯水年水資源短缺程度更為嚴重的原因。(2) 土地利用：土地利用變化對水資源短缺程度的影響主要有兩個方面：一是下墊面的變化直接影響水文循環過程，繼而影響產水量；另一是耗水量和土地利用格局密切相關，所以耕地和建設用地更容易出現水資源短缺的情況，這就是受人類活動干擾較大區域的水資源短缺程度更嚴重的原因。

將渭河流域不同情景下的產水量、耗水量和河道生態基流量代入式(8)，對渭河流域整體的水資源短缺狀況進行綜合評價，評價結果見表2。

表2 不同情景下渭河流域的水資源短缺評價

渭河流域在豐水年(情景1和情景2)的水資源盈缺指數WSi為4，平水年(情景3和情景4)WSi為2，屬于水資源盈余流域。而枯水年(情景5和情景6)WSi僅為0.8，屬于水資源短缺流域，產水量無法同時滿足河道生態基流量和人類耗水量的需求，說明枯水期是渭河流域河道生態基流保障的關鍵時期。

2.3 氣候變化對水資源短缺的影響

選擇降水量代表氣候因素，以情景2(豐水年)、情景4(平水年)和情景6(枯水年)的數據為基礎，研究氣候因素對水資源短缺時空變化的驅動特征。用GeoDa模型得到降水量—水資源盈缺指數WSi的雙變量莫蘭指數I為0.206～0.277，莫蘭指數為正值，說明降水量與WSi呈正相關關系，存在空間集聚效應，即降水量高的區域WSi高，水資源盈余，降水量低的地區WSi低，水資源短缺。p值為0.001，置信度達到99%，說明渭河流域降水量變化對水資源短缺驅動的空間集聚性顯著。

用GeoDa模型得到降水量—WSi的雙變量LISA集聚圖(圖3)，降水量與WSi的空間相關性以高—高型和低—低型區域為主，高—高型區域分布在北洛河流域和秦嶺山區，這些區域的植被覆蓋度好，降水豐沛，耗水量低，水資源充盈。低—低型區域分布在涇河流域北部和天水地區，降水量少而蒸發量高導致水資源短缺。

圖3 渭河流域降水量-WSi的空間集聚分布

2.4 土地利用變化對水資源短缺的影響

根據土地利用類型的分級指數[25]，未利用地為1級，林地、草地、水域為2級，耕地為3級，建設用地為4級。本文用分級指數表征土地的開發利用程度，即分級指數越高，土地開發利用程度越高。以情景5(1980年)和情景6(2015年)的土地利用/覆被數據為基礎，研究土地利用因素對水資源短缺時空變化的驅動特征。用GeoDa模型得到土地利用強度—水資源盈缺指數WSi的雙變量莫蘭指數I為-0.080～-0.045，莫蘭指數為負值，說明土地利用強度與WSi是負相關關系，存在空間差異效應，即土地利用強度低的區域WSi高，水資源盈余；土地開發利用強度增加，WSi降低，水資源短缺。p值為0.007～0.075，置信度為90%，說明渭河流域土地利用變化對水資源短缺的驅動顯著性一般。

用GeoDa模型得到土地利用強度-WSi的雙變量LISA集聚圖(圖4)，土地利用與WSi的空間相關性高—低型區域集中在天水地區，人口密集，土地利用開發程度高，水資源短缺；低—高型區域分布在北洛河流域和秦嶺山區，這些區域的植被覆蓋度好，耗水量低，水資源充盈；低—低型區域集聚在涇河流域北部，土地開發程度低，但仍存在水資源短缺的情況。

圖4 渭河流域土地利用-WSi的空間集聚分布

3 結 論

(1) 渭河流域的產水量由北向南逐漸遞增，涇河流域和北洛河流域的單位面積產水量多低于1 000 m3/hm2，渭河干流可以達到1 500 m3/hm2。

(2) 渭河流域的水資源短缺區域集中在渭河干流和涇河流域北部；枯水年渭河流域的水資源盈缺指數WSi僅為0.8，是河道生態基流保障的關鍵時期。

(3) 以降水量為代表的氣候變化與WSi存在顯著的空間集聚效應，空間相關性以高—高型和低—低型區域為主，高—高型區域分布在北洛河流域和秦嶺山區，低—低型區域分布在涇河流域北部和天水地區；土地利用變化與WSi存在的空間差異效應顯著性一般。